AssertionError and suggestions required to improve my model

Questions
1

Dear all, I am getting an error of AssertionError. Can You please help me.

I also need your suggestions to improve my model because I am new with PyMC3 and theano

Please find the code below:
Thanking you all

# This is the code used for pymc3 model

import numpy as np
import pymc3 as pm
import theano.tensor as tt


def test_model(p,q,r,stations):
    
    x,y=stations
    Azim=p
    DIP=q
    Wid=r
    U1 = np.zeros((np.size(x),3),dtype=object)
    
    for i in range(1,np.size(x)):
        U1[i,:]=get_Okada(Azim, DIP,Wid,y[0,i],x[0,i],0)
        
            
    return U1
        

def get_Okada(Azim, DIP,Wid,x,y,z):
    
   
    xrot=tt.stacklists([[tt.cos(Azim*np.pi/180), -tt.sin(Azim*np.pi/180), 0],\
                          [tt.sin(Azim*np.pi/180), tt.cos(Azim*np.pi/180), 0], 
                          [0, 0, 1]])
   
    new_x= xrot*(np.asarray([[x],[y],[z]]))
    
    h  = new_x[2][0]
    
    st=tt.sin(DIP*np.pi/180)
   
    w  = Wid/1000
       
    xpa =tt.exp(h)
    
    ypa = h + w*st
    
      
    (tmpX,tmpY,AZ) = FAULT(xpa,ypa,st)    
    (AX,AY,tmpZ) = FAULT(xpa,ypa,st)
    U1 =np.matrix([AX, AY, AZ])
    
    U2 = (xrot.transpose() * U1.transpose()).transpose()

    return U2


    
def FAULT(xpa,ypa,st):
    
    qs1=xpa**2 + ypa/2.0
    #import sys
    #sys.path.insert(0,'C:/Users/singh/Dropbox/Catania/Python_codes')
   
    [ux1,uy1,uz1] = DISLOK(qs1,st)
    
    ax=ux1/np.pi
    ay=uy1/2/np.pi
    az=uz1/2/np.pi
        
    return ax,ay,az
    
def DISLOK(qsi, st):
    
    r = pm.math.sqrt(qsi**2)
    xx = pm.math.sqrt(qsi**2 +78)
    
    atqe = tt.arctan(qsi * xx/r)
        
    fux=fuy=fuz=0.0
    
    ai_4 = (pm.math.log(r ) - st*pm.math.log(r + st))/st
    
    fuz = ai_4*st**2
    fux = fux/ai_4*st
    fuy = st*atqe  
        
    return fux, fuy,fuz 


### Now I call it as basic model
x0=np.random.rand(1,21)
y0=np.random.rand(1,21)
stations=x0,y0

import pymc3 as pm
import test_pymc_model as tpm

basic_model = pm.Model()

with basic_model:
    
    
    Azim=pm.Uniform('Azim',lower=0,upper=360)
    DIP=pm.Uniform('DIP',lower=10,upper=89.9)
    Wid=pm.Uniform('Wid',lower=100,upper=10000)
    sigma = pm.HalfNormal('sigma', sd=1)
    mu1 = tpm.test_model(Azim,DIP,Wid,stations)
    Y_obs = pm.Normal('Y_obs', mu=mu1, sd=sigma, observed=U)

I got the following error

Traceback (most recent call last):
  File "test_pymc.py", line 40, in <module>
    mu1 = tpm.test_model(Azim,DIP,Wid,stations)
  File "F:\Python_codes\dummy_codes\test_pymc_model.py", line 16, in test_model
    U1[i,:]=get_Okada(Azim, DIP,Wid,y[0,i],x[0,i],0)
  File "F:\Python_codes\dummy_codes\test_pymc_model.py", line 46, in get_Okada
    U2 = (xrot.transpose() * U1.transpose()).transpose()
  File "C:\Users\singh.CT\AppData\Local\Continuum\anaconda3\lib\site-packages\theano\tensor\var.py", line 155, in __mul__
    return theano.tensor.mul(self, other)
  File "C:\Users\singh.CT\AppData\Local\Continuum\anaconda3\lib\site-packages\theano\gof\op.py", line 615, in __call__
    node = self.make_node(*inputs, **kwargs)
  File "C:\Users\singh.CT\AppData\Local\Continuum\anaconda3\lib\site-packages\theano\tensor\elemwise.py", line 480, in make_node
    inputs = list(map(as_tensor_variable, inputs))
  File "C:\Users\singh.CT\AppData\Local\Continuum\anaconda3\lib\site-packages\theano\tensor\basic.py", line 194, in as_tensor_variable
    return constant(x, name=name, ndim=ndim)
  File "C:\Users\singh.CT\AppData\Local\Continuum\anaconda3\lib\site-packages\theano\tensor\basic.py", line 232, in constant
    x_ = scal.convert(x, dtype=dtype)
  File "C:\Users\singh.CT\AppData\Local\Continuum\anaconda3\lib\site-packages\theano\scalar\basic.py", line 284, in convert
    assert type(x_) in [np.ndarray, np.memmap]
AssertionError
2

In case you need it , this is my actual model with that I calculate my data (U). I want to invert U for estimating the (p,q,r OR Azim, DIP, Wid)

import numpy as np
def test_model(p,q,r,stations):
x,y=stations
Azim=p
DIP=q
Wid=r
U1 = np.zeros((np.size(x),3),dtype=object)

for i in range(1,np.size(x)):
    U1[i,:]=getOkada(Azim, DIP,Wid,y[0,i],x[0,i],0)
    
 return U1

def getOkada(Azim, DIP,Wid,x,y,z):

xrot=np.matrix([[np.cos(np.radians(Azim)),\
      -np.sin(np.radians(Azim)), 0],\
      [np.sin(np.radians(Azim)),\
      np.cos(np.radians(Azim)),\
      0], [0, 0, 1]])

new_x= xrot*(np.asarray([[x],[y],[z]]))

h  = new_x[2][0]

st=np.sin(DIP*np.pi/180)

w  = Wid/1000
   
xpa =np.exp(h)

ypa = h + w*st   
  
(tmpX,tmpY,AZ) = FAULT(xpa,ypa,st)    
(AX,AY,tmpZ) = FAULT(xpa,ypa,st)
AX=float(AX)
AY=float(AY)
AZ=float(AZ)
U1 =np.matrix([AX, AY, AZ])

U2 = (xrot.transpose() * U1.transpose()).transpose()


return U2

def FAULT(xpa,ypa,st):
qs1=xpa**2 + ypa/2.0

[ux1,uy1,uz1] = DISLOK(qs1,st)

ax=ux1/np.pi
ay=uy1/2/np.pi
az=uz1/2/np.pi
    
return ax,ay,az

def DISLOK(qsi, st):
r = np.sqrt(qsi2)
xx = np.sqrt(qsi2 +78)
atqe = np.arctan(qsi * xx/r)

fux=fuy=fuz=0.0

ai_4 = (np.log(r ) - st*np.log(r + st))/st

fuz = ai_4*st**2
fux = fux/ai_4*st
fuy = st*atqe  
    
return fux, fuy,fuz

I call it like this

import numpy as np
import matplotlib.pyplot as plt

p=89.9
q=75
r=1000

x0=np.random.rand(1,21)
y0=np.random.rand(1,21)
stations=x0,y0

import test_model as TM

U=TM.test_model(p,q,r,stations)

3

numpy is not aware of theano tensors so the quoted line makes a matrix of dtype=object. You need to replace that line with an appropriate theano.stack.

Furthermore U1 should be created as a tensor of floats using theano.tensor.zeros, and finally you cannot set its elements like you did in the for loop, you need to use inc_subtensor or one of its variants from theano.tensor.

4

Thank You very much for your reply. I am trying to adopt your suggestions. However, is there any alternative way to get rid of the loop.

Actually I have the N number of data points which are calculated by calling the model N times. So, when writing the following line in pymc3 model:

Y_obs = pm.Normal(‘Y_obs’, mu=mu1, sd=sigma, observed=U)

the length/size of mu1 and U should be the same. For that reason, I have to use loop so mu1 return with the same length/size as of U. My precise question, can I some how calculate mu1 iteratively within pymc3 model ? if its possible then I don’t need to use loop inside my customized forward model. I would appriciate if you could forward some particular example of vectorized model or where observed data has to be calculated iteratively.
Thank You!

5

I think that what you are trying to do is called vectorize the loop. The answer can be yes, but you’ll have to work hard to make get_Okada to work with

  • Azim, DIP and Wid that are scalar floats or scalar random variables.
  • x and y that are 1 dimensional arrays.

Your function should then output an (N, 3) array that would set all the elements of U1 without using a for loop. One thing I hadn’t noticed before was that mu1 ends up being a 2 dimensional tensor, so your observations must also be 2 dimensional, and in particular have (N, 3) shape!

1 Like