Can anyone tell me what's wrong in my code below ?
I am a casual user of pymc2, generally for solving physical equations. I have troubles to adapt a fit to pymc3 and the documentation seems to me unclear. Also I did not recognize my problem on forums, probably because I don’t know what is my problem…
I use the find_MAP method to get a first guess of fitted values but this first guess is completly wrong (not even inside the physical limits) and a warning tells me there are discrete variables (which is wrong), implying the gradient is not available.
The aim is to fit some parameters in a diffusion equation: here alpha0, alpha1 and epsilon, that are continuous and a priori uniformly distributed. During a long debugging time I anti-optimized the code, so I don’t think the code is interesting by itself. Just know that the pymc2 version is ok and works fine. As I don’t know where the problem(s) is(are), I also give the inside of the ‘simul_DifferentialEq’ function, but the pymc3 stuff is below the corresponding comment.
import numpy as np
from scipy.interpolate import interp1d
import pymc3 as pm3
import theano.tensor as tt
import theano.compile
import config
@theano.compile.ops.as_op(itypes=[tt.dscalar,tt.dscalar,tt.dscalar],
otypes=[tt.dvector])
def simul_DifferentialEq(alpha0,alpha1,epsilon):
observed_depth = np.array([0.5,1.5,3,5,7,9,13,17])# in cm
observed_values = np.array([6.25,2.75,1.25,1.25,1.5,1.75,1.5,1])# mmol.l-1
#useful ?
observed_values = observed_values[np.argsort(-observed_depth)]
observed_depth = -np.sort(-observed_depth)
depth = config.depth
matA = config.matA
matB = config.matB
matC = config.matC
concentration = config.concentration
alpha = alpha0 * np.exp(-alpha1*depth) # in day^-1
# simplification for a Constant phi
phi = np.empty(len(depth))
phi[:]=0.6
#########################
beta = config.beta * phi * epsilon # dimensionless
delta = beta * phi / config.Deltax
eta = alpha * config.Deltat
f1 = f2 = delta
f3 = 1 - 2*delta + eta/2
matB[0,0] = f1[0] - eta[0]/2 +1
matB[0,1] = matA[0,1] = -2*f1[0]
matB[0,2] = matA[0,2] = delta[0]
matB[-1,-1] = f2[-1] + eta[-1]/2 +1
matB[-1,-2] = matA[-1,-2] = -2*f2[-1]
matB[-1,-3] = matA[-1,-3] = delta[-1]
matB[range(1,concentration.sizex-1),
range(1,concentration.sizex-1)] = \
f3[1:concentration.sizex-1]
matA[range(1,concentration.sizex-1),
range(concentration.sizex-2)] = \
matB[range(1,concentration.sizex-1),
range(concentration.sizex-2)] = \
f1[1:concentration.sizex-1]
matA[range(1,concentration.sizex-1),
range(2,concentration.sizex)] = \
matB[range(1,concentration.sizex-1),
range(2,concentration.sizex)] = \
f2[1:concentration.sizex-1]
matB[range(1,concentration.sizex),0] = -eta[1:]
matA[range(concentration.sizex),
range(concentration.sizex)] = \
matB[range(concentration.sizex),
range(concentration.sizex)] -2
matA[0,0] += eta[0]
matC = np.dot(np.linalg.matrix_power(-matA,-1),matB)
for tcount in range(concentration.sizet-1):
#the variable 'temp' has no interest (just convenient for debugging)
temp = np.dot(matC,concentration.values[:,tcount])
# condition limit
temp[0] = config.C0
# a priori useless (but convenient for debugging))
temp[np.where(temp>config.C0)] = config.C0
# everything for that...
concentration.values[:,tcount+1] = temp
interpolated_concentration = interp1d(depth,concentration.values[:,-1])
return interpolated_concentration(observed_depth)
# the pymc3 stuff is below
model = pm3.Model()
with model:
alpha0 = pm3.Uniform("alpha0",-2,0)
alpha1 = pm3.Uniform("alpha1",-1,2)
epsilon = pm3.Uniform("epsilon",0.1,15)
DifferentialEq = simul_DifferentialEq(alpha0,alpha1,epsilon)
# it is awkward to repeat observed values
#some previous tries made me think it could solve the problem but it didn't
observed_depth = np.array([0.5,1.5,3,5,7,9,13,17])# in cm
observed_values = np.array([6.25,2.75,1.25,1.25,1.5,1.75,1.5,1])# mmol.l-1
# useful ?
observed_values = observed_values[np.argsort(-observed_depth)]
observed_depth = -np.sort(-observed_depth)
obs = pm3.Normal('obs', mu=DifferentialEq, sd=0.1, observed=observed_values)
print('running test170127, find_MAP...')
testfindmap = pm3.find_MAP()
Thank you for your attention, the content of the config.py is:
C0=Cowl0 = 10 # in mmol/l: concentration at the surface (at t=0), sometimes noted C0
Dsw = 1.6 # in cm^2.d-1
Cdefault = 1e-10 # concentration at t=0, depth>0
# maximum depth and time in the simulation for solving the ED (assuming it begins at x=t=0)
maxdepth = 17 # in cm
maxtime = 1 # in day
#steps in depth and time in the simulation for solving the ED
Deltax = 0.05# in cm
Deltat = 0.02# in day
##############################################
# internal cooking
from numpy import arange, empty, zeros
from solve_ED_crank import sph_2Dfunct
depth = (arange(maxdepth/Deltax +1))*Deltax # in cm
time = (arange(maxtime/Deltat +1))*Deltat # in day
beta = Dsw * Deltat / (2 * Deltax)
matA = zeros([len(depth),len(depth)])
matB = zeros([len(depth),len(depth)])
matC = empty([len(depth),len(depth)])
concentration_t0 = empty(len(depth))
concentration_t0[1:] = Cdefault
concentration_t0[0] = Cowl0
concentration = sph_2Dfunct(sizex=len(depth),
sizet=len(time),
firstline=concentration_t0)
comment tuesday the 07/02 at ~12:30.
I replaced the last line (ie. the find_MAP stuff) with:
pm3.sample(500)
when I run the main code, I get:
Auto-assigning NUTS sampler...
INFO:pymc3:Auto-assigning NUTS sampler...
Initializing NUTS using advi...
INFO:pymc3:Initializing NUTS using advi...
Traceback (most recent call last):
File "<ipython-input-1-8395e07601b2>", line 1, in <module>
runfile('/Users/steph/work/profiles/profiles-pymc/test170127.py', wdir='/Users/steph/work/profiles/profiles-pymc')
File "/Users/steph/anaconda/lib/python3.5/site-packages/spyder/utils/site/sitecustomize.py", line 866, in runfile
execfile(filename, namespace)
File "/Users/steph/anaconda/lib/python3.5/site-packages/spyder/utils/site/sitecustomize.py", line 102, in execfile
exec(compile(f.read(), filename, 'exec'), namespace)
File "/Users/steph/work/profiles/profiles-pymc/test170127.py", line 84, in <module>
pm3.sample(500)
File "/Users/steph/anaconda/lib/python3.5/site-packages/pymc3/sampling.py", line 149, in sample
start_, step = init_nuts(init=init, n_init=n_init, model=model)
File "/Users/steph/anaconda/lib/python3.5/site-packages/pymc3/sampling.py", line 434, in init_nuts
v_params = pm.variational.advi(n=n_init)
File "/Users/steph/anaconda/lib/python3.5/site-packages/pymc3/variational/advi.py", line 139, in advi
updates = optimizer(loss=-1 * elbo, param=[uw_shared])
File "/Users/steph/anaconda/lib/python3.5/site-packages/pymc3/variational/advi.py", line 259, in optimizer
grad = tt.grad(loss, param_)
File "/Users/steph/anaconda/lib/python3.5/site-packages/theano/gradient.py", line 561, in grad
grad_dict, wrt, cost_name)
File "/Users/steph/anaconda/lib/python3.5/site-packages/theano/gradient.py", line 1324, in _populate_grad_dict
rval = [access_grad_cache(elem) for elem in wrt]
File "/Users/steph/anaconda/lib/python3.5/site-packages/theano/gradient.py", line 1324, in <listcomp>
rval = [access_grad_cache(elem) for elem in wrt]
File "/Users/steph/anaconda/lib/python3.5/site-packages/theano/gradient.py", line 1279, in access_grad_cache
term = access_term_cache(node)[idx]
File "/Users/steph/anaconda/lib/python3.5/site-packages/theano/gradient.py", line 973, in access_term_cache
output_grads = [access_grad_cache(var) for var in node.outputs]
File "/Users/steph/anaconda/lib/python3.5/site-packages/theano/gradient.py", line 973, in <listcomp>
output_grads = [access_grad_cache(var) for var in node.outputs]
File "/Users/steph/anaconda/lib/python3.5/site-packages/theano/gradient.py", line 1279, in access_grad_cache
term = access_term_cache(node)[idx]
File "/Users/steph/anaconda/lib/python3.5/site-packages/theano/gradient.py", line 973, in access_term_cache
output_grads = [access_grad_cache(var) for var in node.outputs]
File "/Users/steph/anaconda/lib/python3.5/site-packages/theano/gradient.py", line 973, in <listcomp>
output_grads = [access_grad_cache(var) for var in node.outputs]
File "/Users/steph/anaconda/lib/python3.5/site-packages/theano/gradient.py", line 1279, in access_grad_cache
term = access_term_cache(node)[idx]
File "/Users/steph/anaconda/lib/python3.5/site-packages/theano/gradient.py", line 973, in access_term_cache
output_grads = [access_grad_cache(var) for var in node.outputs]
File "/Users/steph/anaconda/lib/python3.5/site-packages/theano/gradient.py", line 973, in <listcomp>
output_grads = [access_grad_cache(var) for var in node.outputs]
File "/Users/steph/anaconda/lib/python3.5/site-packages/theano/gradient.py", line 1279, in access_grad_cache
term = access_term_cache(node)[idx]
File "/Users/steph/anaconda/lib/python3.5/site-packages/theano/gradient.py", line 973, in access_term_cache
output_grads = [access_grad_cache(var) for var in node.outputs]
File "/Users/steph/anaconda/lib/python3.5/site-packages/theano/gradient.py", line 973, in <listcomp>
output_grads = [access_grad_cache(var) for var in node.outputs]
File "/Users/steph/anaconda/lib/python3.5/site-packages/theano/gradient.py", line 1279, in access_grad_cache
term = access_term_cache(node)[idx]
File "/Users/steph/anaconda/lib/python3.5/site-packages/theano/gradient.py", line 973, in access_term_cache
output_grads = [access_grad_cache(var) for var in node.outputs]
File "/Users/steph/anaconda/lib/python3.5/site-packages/theano/gradient.py", line 973, in <listcomp>
output_grads = [access_grad_cache(var) for var in node.outputs]
File "/Users/steph/anaconda/lib/python3.5/site-packages/theano/gradient.py", line 1279, in access_grad_cache
term = access_term_cache(node)[idx]
File "/Users/steph/anaconda/lib/python3.5/site-packages/theano/gradient.py", line 973, in access_term_cache
output_grads = [access_grad_cache(var) for var in node.outputs]
File "/Users/steph/anaconda/lib/python3.5/site-packages/theano/gradient.py", line 973, in <listcomp>
output_grads = [access_grad_cache(var) for var in node.outputs]
File "/Users/steph/anaconda/lib/python3.5/site-packages/theano/gradient.py", line 1279, in access_grad_cache
term = access_term_cache(node)[idx]
File "/Users/steph/anaconda/lib/python3.5/site-packages/theano/gradient.py", line 973, in access_term_cache
output_grads = [access_grad_cache(var) for var in node.outputs]
File "/Users/steph/anaconda/lib/python3.5/site-packages/theano/gradient.py", line 973, in <listcomp>
output_grads = [access_grad_cache(var) for var in node.outputs]
File "/Users/steph/anaconda/lib/python3.5/site-packages/theano/gradient.py", line 1279, in access_grad_cache
term = access_term_cache(node)[idx]
File "/Users/steph/anaconda/lib/python3.5/site-packages/theano/gradient.py", line 973, in access_term_cache
output_grads = [access_grad_cache(var) for var in node.outputs]
File "/Users/steph/anaconda/lib/python3.5/site-packages/theano/gradient.py", line 973, in <listcomp>
output_grads = [access_grad_cache(var) for var in node.outputs]
File "/Users/steph/anaconda/lib/python3.5/site-packages/theano/gradient.py", line 1279, in access_grad_cache
term = access_term_cache(node)[idx]
File "/Users/steph/anaconda/lib/python3.5/site-packages/theano/gradient.py", line 973, in access_term_cache
output_grads = [access_grad_cache(var) for var in node.outputs]
File "/Users/steph/anaconda/lib/python3.5/site-packages/theano/gradient.py", line 973, in <listcomp>
output_grads = [access_grad_cache(var) for var in node.outputs]
File "/Users/steph/anaconda/lib/python3.5/site-packages/theano/gradient.py", line 1279, in access_grad_cache
term = access_term_cache(node)[idx]
File "/Users/steph/anaconda/lib/python3.5/site-packages/theano/gradient.py", line 973, in access_term_cache
output_grads = [access_grad_cache(var) for var in node.outputs]
File "/Users/steph/anaconda/lib/python3.5/site-packages/theano/gradient.py", line 973, in <listcomp>
output_grads = [access_grad_cache(var) for var in node.outputs]
File "/Users/steph/anaconda/lib/python3.5/site-packages/theano/gradient.py", line 1279, in access_grad_cache
term = access_term_cache(node)[idx]
File "/Users/steph/anaconda/lib/python3.5/site-packages/theano/gradient.py", line 973, in access_term_cache
output_grads = [access_grad_cache(var) for var in node.outputs]
File "/Users/steph/anaconda/lib/python3.5/site-packages/theano/gradient.py", line 973, in <listcomp>
output_grads = [access_grad_cache(var) for var in node.outputs]
File "/Users/steph/anaconda/lib/python3.5/site-packages/theano/gradient.py", line 1279, in access_grad_cache
term = access_term_cache(node)[idx]
File "/Users/steph/anaconda/lib/python3.5/site-packages/theano/gradient.py", line 973, in access_term_cache
output_grads = [access_grad_cache(var) for var in node.outputs]
File "/Users/steph/anaconda/lib/python3.5/site-packages/theano/gradient.py", line 973, in <listcomp>
output_grads = [access_grad_cache(var) for var in node.outputs]
File "/Users/steph/anaconda/lib/python3.5/site-packages/theano/gradient.py", line 1279, in access_grad_cache
term = access_term_cache(node)[idx]
File "/Users/steph/anaconda/lib/python3.5/site-packages/theano/gradient.py", line 1113, in access_term_cache
input_grads = node.op.grad(inputs, new_output_grads)
AttributeError: 'FromFunctionOp' object has no attribute 'grad'
I would add: if I keep the call to find_MAP, the code runs without any error but the resulting values look absurd and I get this double-warning:
Warning: gradient not available.(E.g. vars contains discrete variables). MAP estimates may not be accurate for the default parameters. Defaulting to non-gradient minimization fmin_powell.WARNING:pymc3:Warning: gradient not available.(E.g. vars contains discrete variables). MAP estimates may not be accurate for the default parameters. Defaulting to non-gradient minimization fmin_powell.
Optimization terminated successfully.
Current function value: 36.569283
Iterations: 10
Function evaluations: 415
