I have a model in which rates of a random (Poisson) process are distributed according to a function p_nu(gamma,delta,nu_max). With pymc3/theano I want to obtain the posterior distribution of gamma,delta,nu_max when observing N_AP events (N_AP being a vector, each entry being a measurement of one unit, whose rate is drawn from p_nu).
I'm running into issues when evaluating the jacobian of the scan function, required to iterate over multiple measurements given in N_AP - details are described below.
This is my first post on stackoverflow, as well as part of my first dive into theano and I've been trying to break the model down to a minimum example, providing explanations where needed and including all code snippets necessary to run the example. Thus, I'm also happy to receive feedback on how the issue was posted and if things should be more concise / detailed, etc, or if my way of using theano could be improved.
Used imports for this example are
import numpy as np
from scipy.integrate import quad
import theano
import theano.tensor as T
from theano.graph.op import Op
from theano.graph.basic import Apply
For this, I have built a custom theano Op (or rather slightly adapted the one found in Custom Theano Op to do numerical integration), as it involves integration over (complete code below)
p_nu(nu,gamma,delta,nu_max)*poisson_spikes(nu,N,T_total)
with function definitions
def p_nu(NU,gamma,delta,nu_max):
return gamma / ( nu_max * np.sqrt( -np.pi * np.log( NU / nu_max ) ) ) * \
np.exp( - delta**2/2.) * ( NU / nu_max )**(gamma**2 - 1) * \
np.cosh( gamma * delta * np.sqrt( -2 * np.log( NU / nu_max) ) )
def poisson_spikes(nu,N,T_total):
return (nu*T_total)**N / T.gamma(N+1) * T.exp(-nu*T_total)
and variables
T_total = 600. # measurement time
N_AP = T.lvector('N_AP') # number of events (for vector values)
N = T.dscalar('N') # number of events (for single values)
nu = T.dscalar('nu') # event rate; is integrated out in the process
gamma = T.dscalar('gamma') # model parameter
delta = T.dscalar('delta') # model parameter
nu_max = T.dscalar('nu_max') # model parameter
With the integration Op defined as (using a couple of prints for debugging purposes)
class integrateOut(Op):
"""
Integrate out a variable from an expression, computing
the definite integral w.r.t. the variable specified
!!! Only implemented in this for scalars !!!
"""
def __init__(self,f,t,*args,**kwargs):
super(integrateOut,self).__init__()
self.f = f
self.t = t
def make_node(self,*inputs):
self.fvars=list(inputs)
#print(f'[make node]: inputs: {inputs}, fvars:{self.fvars}')
# This will fail when taking the gradient... don't be concerned
try:
self.gradF = T.grad(self.f,self.fvars)
print(f'[make node (fun)]: inputs: {inputs}, fvars:{self.fvars}')
except:
self.gradF = None
print(f'[make node (grad)]: inputs: {inputs}, fvars:{self.fvars}')
return Apply(self,self.fvars,[T.dscalar().type()])
def perform(self,node, inputs, output_storage):
if self.gradF is None:
print(f'[perform (grad)]: inputs: {inputs}, fvars:{self.fvars}')
else:
print(f'[perform (fun)]: inputs: {inputs}, fvars:{self.fvars}')
# create a function to evaluate the integral
f = theano.function([self.t]+self.fvars,self.f)
# integrate the function from 0 to maximum firing rate nu_max
nu_max = inputs[3]
output_storage[0][0] = np.array(quad(f,0,nu_max,args=tuple(inputs))[0],dtype='float64')
def grad(self,inputs,output_grads):
nu_max = inputs[3]
return [integrateOut(g,self.t,0,nu_max)(*inputs)*output_grads[0] \
for g in self.gradF]
The code works smoothly for calculating the actual likelihood, but as I'd like to use NUTS sampling, I am also required to evaluate the grad / jacobian. And here comes the tricky part, or rather the part that I've been trying to wrap my head around for the past couple of days: Evaluating the jacobian for a scalar, defined as
p_N_AP = integrateOut(p_nu(nu,gamma,delta,nu_max)*poisson_spikes(nu,N,T_total),nu)(N,gamma,delta,nu_max)
pGrad = T.jacobian(p_N_AP,[gamma,delta,nu_max])
funcGrad = theano.function([N,gamma,delta,nu_max],pGrad)
grad_vals = funcGrad(0,1.2,4.8,30.)
with N being a single, scalar valued datapoint, works.
Evaluating it for N_AP, a vector of datapoints, the functions are defined as
p_N_AP_arr, updates = theano.scan(
fn = lambda N, gamma, delta, nu_max : integrateOut(p_nu(nu,gamma,delta,nu_max)*poisson_spikes(nu,N,T_total),nu)(N,gamma,delta,nu_max),
sequences=[N_AP],
non_sequences=[gamma,delta,nu_max]
)
pGrad = T.jacobian(p_N_AP_arr,[gamma,delta,nu_max])
funcGrad = theano.function([N_AP,gamma,delta,nu_max],pGrad) ### somehow in here, copies are made, which dont fit ..
grad_vals = funcGrad([0,3,5],1.2,4.8,30.)
fails with the error:
[make node (fun)]: inputs: (N_AP[t], gamma, delta, nu_max), fvars:[N_AP[t], gamma, delta, nu_max]
[make node (grad)]: inputs: (N_AP[t], gamma_copy, delta_copy, nu_max_copy), fvars:[N_AP[t], gamma_copy, delta_copy, nu_max_copy]
... couple of same outputs here ...
[make node (grad)]: inputs: (N_AP[t], gamma_copy, delta_copy, nu_max_copy), fvars:[N_AP[t], gamma_copy, delta_copy, nu_max_copy]
[perform (fun)]: inputs: [array(0), array(1.2), array(4.8), array(30.)], fvars:[N_AP[t]012, gamma, delta, nu_max]
[perform (fun)]: inputs: [array(3), array(1.2), array(4.8), array(30.)], fvars:[N_AP[t]012, gamma, delta, nu_max]
[perform (fun)]: inputs: [array(5), array(1.2), array(4.8), array(30.)], fvars:[N_AP[t]012, gamma, delta, nu_max]
[perform (grad)]: inputs: [array(5), array(1.2), array(4.8), array(30.)], fvars:[N_AP[t]012, gamma_copy012, delta_copy012, nu_max_copy012]
---------------------------------------------------------------------------
UnusedInputError Traceback (most recent call last)
File scan_perform.pyx:399, in theano.scan.scan_perform.perform()
File ~/Data/Programmes/anaconda3/envs/inference/lib/python3.9/site-packages/theano/graph/op.py:476, in Op.make_py_thunk.<locals>.rval(p, i, o, n)
475 def rval(p=p, i=node_input_storage, o=node_output_storage, n=node):
--> 476 r = p(n, [x[0] for x in i], o)
477 for o in node.outputs:
Input In [29], in integrateOut.perform(self, node, inputs, output_storage)
30 print(f'[perform (fun)]: inputs: {inputs}, fvars:{self.fvars}')
---> 32 f = theano.function([self.t]+self.fvars,self.f)
34 ## integrate the function from 0 to maximum firing rate nu_max
File ~/Data/Programmes/anaconda3/envs/inference/lib/python3.9/site-packages/theano/compile/function/__init__.py:337, in function(inputs, outputs, mode, updates, givens, no_default_updates, accept_inplace, name, rebuild_strict, allow_input_downcast, profile, on_unused_input)
334 else:
335 # note: pfunc will also call orig_function -- orig_function is
336 # a choke point that all compilation must pass through
--> 337 fn = pfunc(
338 params=inputs,
339 outputs=outputs,
340 mode=mode,
341 updates=updates,
342 givens=givens,
343 no_default_updates=no_default_updates,
344 accept_inplace=accept_inplace,
345 name=name,
346 rebuild_strict=rebuild_strict,
347 allow_input_downcast=allow_input_downcast,
348 on_unused_input=on_unused_input,
349 profile=profile,
350 output_keys=output_keys,
351 )
352 return fn
File ~/Data/Programmes/anaconda3/envs/inference/lib/python3.9/site-packages/theano/compile/function/pfunc.py:524, in pfunc(params, outputs, mode, updates, givens, no_default_updates, accept_inplace, name, rebuild_strict, allow_input_downcast, profile, on_unused_input, output_keys)
522 inputs.append(si)
--> 524 return orig_function(
525 inputs,
526 cloned_outputs,
527 mode,
528 accept_inplace=accept_inplace,
529 name=name,
530 profile=profile,
531 on_unused_input=on_unused_input,
532 output_keys=output_keys,
533 )
File ~/Data/Programmes/anaconda3/envs/inference/lib/python3.9/site-packages/theano/compile/function/types.py:1970, in orig_function(inputs, outputs, mode, accept_inplace, name, profile, on_unused_input, output_keys)
1969 Maker = getattr(mode, "function_maker", FunctionMaker)
-> 1970 m = Maker(
1971 inputs,
1972 outputs,
1973 mode,
1974 accept_inplace=accept_inplace,
1975 profile=profile,
1976 on_unused_input=on_unused_input,
1977 output_keys=output_keys,
1978 name=name,
1979 )
1980 with config.change_flags(compute_test_value="off"):
File ~/Data/Programmes/anaconda3/envs/inference/lib/python3.9/site-packages/theano/compile/function/types.py:1573, in FunctionMaker.__init__(self, inputs, outputs, mode, accept_inplace, function_builder, profile, on_unused_input, fgraph, output_keys, name)
1572 # Check if some input variables are unused
-> 1573 self._check_unused_inputs(inputs, outputs, on_unused_input)
1575 # Make a list of (SymbolicInput|SymblicInputKits, indices,
1576 # [SymbolicInput,...]), one tuple for each input. (See
1577 # Function.indices for more details)
File ~/Data/Programmes/anaconda3/envs/inference/lib/python3.9/site-packages/theano/compile/function/types.py:1745, in FunctionMaker._check_unused_inputs(self, inputs, outputs, on_unused_input)
1744 elif on_unused_input == "raise":
-> 1745 raise UnusedInputError(msg % (inputs.index(i), i.variable, err_msg))
1746 else:
UnusedInputError: theano.function was asked to create a function computing outputs given certain inputs, but the provided input variable at index 2 is not part of the computational graph needed to compute the outputs: gamma_copy012.
To make this error into a warning, you can pass the parameter on_unused_input='warn' to theano.function. To disable it completely, use on_unused_input='ignore'.
During handling of the above exception, another exception occurred:
AttributeError Traceback (most recent call last)
File scan_perform.pyx:399, in theano.scan.scan_perform.perform()
File ~/Data/Programmes/anaconda3/envs/inference/lib/python3.9/site-packages/theano/scan/op.py:1086, in Scan.make_thunk.<locals>.rval(p, i, o, n, allow_gc)
1083 def rval(
1084 p=p, i=node_input_storage, o=node_output_storage, n=node, allow_gc=allow_gc
1085 ):
-> 1086 r = p(n, [x[0] for x in i], o)
1087 for o in node.outputs:
File ~/Data/Programmes/anaconda3/envs/inference/lib/python3.9/site-packages/theano/scan/op.py:1045, in Scan.make_thunk.<locals>.p(node, args, outs)
1044 def p(node, args, outs):
-> 1045 return scan_perform_ext.perform(
1046 self.n_shared_outs,
1047 self.n_mit_mot_outs,
1048 self.n_seqs,
1049 self.n_mit_mot,
1050 self.n_mit_sot,
1051 self.n_sit_sot,
1052 self.n_nit_sot,
1053 args[0],
1054 self.as_while,
1055 cython_mintaps,
1056 cython_tap_array,
1057 cython_tap_array_len,
1058 cython_vector_seqs,
1059 cython_vector_outs,
1060 cython_mit_mot_out_slices,
1061 cython_mit_mot_out_nslices,
1062 cython_mitmots_preallocated,
1063 cython_inps_is_tensor,
1064 cython_outs_is_tensor,
1065 self.fn.fn,
1066 self.fn,
1067 cython_destroy_map,
1068 args,
1069 outs,
1070 self,
1071 node,
1072 )
File scan_perform.pyx:407, in theano.scan.scan_perform.perform()
AttributeError: 'CVM' object has no attribute 'maker'
During handling of the above exception, another exception occurred:
AttributeError Traceback (most recent call last)
File ~/Data/Programmes/anaconda3/envs/inference/lib/python3.9/site-packages/theano/compile/function/types.py:974, in Function.__call__(self, *args, **kwargs)
972 try:
973 outputs = (
--> 974 self.fn()
975 if output_subset is None
976 else self.fn(output_subset=output_subset)
977 )
978 except Exception:
File ~/Data/Programmes/anaconda3/envs/inference/lib/python3.9/site-packages/theano/scan/op.py:1086, in Scan.make_thunk.<locals>.rval(p, i, o, n, allow_gc)
1083 def rval(
1084 p=p, i=node_input_storage, o=node_output_storage, n=node, allow_gc=allow_gc
1085 ):
-> 1086 r = p(n, [x[0] for x in i], o)
1087 for o in node.outputs:
File ~/Data/Programmes/anaconda3/envs/inference/lib/python3.9/site-packages/theano/scan/op.py:1045, in Scan.make_thunk.<locals>.p(node, args, outs)
1044 def p(node, args, outs):
-> 1045 return scan_perform_ext.perform(
1046 self.n_shared_outs,
1047 self.n_mit_mot_outs,
1048 self.n_seqs,
1049 self.n_mit_mot,
1050 self.n_mit_sot,
1051 self.n_sit_sot,
1052 self.n_nit_sot,
1053 args[0],
1054 self.as_while,
1055 cython_mintaps,
1056 cython_tap_array,
1057 cython_tap_array_len,
1058 cython_vector_seqs,
1059 cython_vector_outs,
1060 cython_mit_mot_out_slices,
1061 cython_mit_mot_out_nslices,
1062 cython_mitmots_preallocated,
1063 cython_inps_is_tensor,
1064 cython_outs_is_tensor,
1065 self.fn.fn,
1066 self.fn,
1067 cython_destroy_map,
1068 args,
1069 outs,
1070 self,
1071 node,
1072 )
File scan_perform.pyx:407, in theano.scan.scan_perform.perform()
AttributeError: 'CVM' object has no attribute 'maker'
During handling of the above exception, another exception occurred:
AttributeError Traceback (most recent call last)
Input In [32], in <cell line: 11>()
9 pGrad = T.jacobian(p_N_AP_arr,[gamma,delta,nu_max])
10 funcGrad = theano.function([N_AP,gamma,delta,nu_max],pGrad) ### somehow in here, copies are made, which dont fit ..
---> 11 grad_vals = funcGrad([0,3,5],1.2,4.8,30.)
File ~/Data/Programmes/anaconda3/envs/inference/lib/python3.9/site-packages/theano/compile/function/types.py:987, in Function.__call__(self, *args, **kwargs)
985 if hasattr(self.fn, "thunks"):
986 thunk = self.fn.thunks[self.fn.position_of_error]
--> 987 raise_with_op(
988 self.maker.fgraph,
989 node=self.fn.nodes[self.fn.position_of_error],
990 thunk=thunk,
991 storage_map=getattr(self.fn, "storage_map", None),
992 )
993 else:
994 # old-style linkers raise their own exceptions
995 raise
File ~/Data/Programmes/anaconda3/envs/inference/lib/python3.9/site-packages/theano/link/utils.py:508, in raise_with_op(fgraph, node, thunk, exc_info, storage_map)
505 warnings.warn(f"{exc_type} error does not allow us to add extra error message")
506 # Some exception need extra parameter in inputs. So forget the
507 # extra long error message in that case.
--> 508 raise exc_value.with_traceback(exc_trace)
File ~/Data/Programmes/anaconda3/envs/inference/lib/python3.9/site-packages/theano/compile/function/types.py:974, in Function.__call__(self, *args, **kwargs)
971 t0_fn = time.time()
972 try:
973 outputs = (
--> 974 self.fn()
975 if output_subset is None
976 else self.fn(output_subset=output_subset)
977 )
978 except Exception:
979 restore_defaults()
File ~/Data/Programmes/anaconda3/envs/inference/lib/python3.9/site-packages/theano/scan/op.py:1086, in Scan.make_thunk.<locals>.rval(p, i, o, n, allow_gc)
1083 def rval(
1084 p=p, i=node_input_storage, o=node_output_storage, n=node, allow_gc=allow_gc
1085 ):
-> 1086 r = p(n, [x[0] for x in i], o)
1087 for o in node.outputs:
1088 compute_map[o][0] = True
File ~/Data/Programmes/anaconda3/envs/inference/lib/python3.9/site-packages/theano/scan/op.py:1045, in Scan.make_thunk.<locals>.p(node, args, outs)
1044 def p(node, args, outs):
-> 1045 return scan_perform_ext.perform(
1046 self.n_shared_outs,
1047 self.n_mit_mot_outs,
1048 self.n_seqs,
1049 self.n_mit_mot,
1050 self.n_mit_sot,
1051 self.n_sit_sot,
1052 self.n_nit_sot,
1053 args[0],
1054 self.as_while,
1055 cython_mintaps,
1056 cython_tap_array,
1057 cython_tap_array_len,
1058 cython_vector_seqs,
1059 cython_vector_outs,
1060 cython_mit_mot_out_slices,
1061 cython_mit_mot_out_nslices,
1062 cython_mitmots_preallocated,
1063 cython_inps_is_tensor,
1064 cython_outs_is_tensor,
1065 self.fn.fn,
1066 self.fn,
1067 cython_destroy_map,
1068 args,
1069 outs,
1070 self,
1071 node,
1072 )
File scan_perform.pyx:407, in theano.scan.scan_perform.perform()
AttributeError: 'CVM' object has no attribute 'maker'
Apply node that caused the error: for{cpu,scan_fn}(Shape_i{0}.0, Subtensor{int64:int64:int8}.0, Shape_i{0}.0, Shape_i{0}.0, Shape_i{0}.0, for{cpu,scan_fn}.0, gamma, delta, nu_max, Shape_i{0}.0, Alloc.0, Subtensor{int64:int64:int64}.0, ScalarFromTensor.0, Alloc.0, Subtensor{int64:int64:int64}.0)
Toposort index: 38
Inputs types: [TensorType(int64, scalar), TensorType(int64, vector), TensorType(int64, scalar), TensorType(int64, scalar), TensorType(int64, scalar), TensorType(float64, vector), TensorType(float64, scalar), TensorType(float64, scalar), TensorType(float64, scalar), TensorType(int64, scalar), TensorType(float64, vector), TensorType(float64, vector), Scalar(int64), TensorType(float64, vector), TensorType(int64, vector)]
Inputs shapes: [(), (3,), (), (), (), (3,), (), (), (), (), (4,), (3,), (), (3,), (3,)]
Inputs strides: [(), (8,), (), (), (), (8,), (), (), (), (), (8,), (-8,), (), (8,), (-8,)]
Inputs values: [array(3), array([0, 1, 2]), array(3), array(3), array(3), array([1.06234782e-09, 3.62980190e-02, 2.32467904e-02]), array(1.2), array(4.8), array(30.), array(3), array([0., 0., 0., 0.]), array([2.32467904e-02, 3.62980190e-02, 1.06234782e-09]), 3, array([0., 0., 0.]), array([5, 3, 0])]
Inputs type_num: [7, 7, 7, 7, 7, 12, 12, 12, 12, 7, 12, 12, 7, 12, 7]
Outputs clients: [['output'], ['output'], ['output']]
HINT: Re-running with most Theano optimization disabled could give you a back-trace of when this node was created. This can be done with by setting the Theano flag 'optimizer=fast_compile'. If that does not work, Theano optimizations can be disabled with 'optimizer=None'.
HINT: Use the Theano flag 'exception_verbosity=high' for a debugprint and storage map footprint of this apply node.
To my understanding, theanos scan creates named copies of the symbolic variables when calling make_node, and when evaluating the function, passes down the given inputs, which do not agree with the automatically named variables anymore, thus throwing an UnusedInputError. I'm somewhat at a loss as to how to fix that and would love to get some insight into how to fix this (and would take some info on the way into how theano decides to create copies of symbolic variables and how their names are given).
I did search stackoverflow / the internet for similar / same issues, but couldn't find anything that helped me get past the described issue
