bransfor
2012-04-18
It looks like adjoint sensitivity analysis is not available from the Simulator class. Is the following code an efficient way to get the adjoint sensitivity of a state at multiple time points? More specifically, can some of the operations be pulled out of the loop, i.e. performed only once?
The example is an ODE with 2 states and 2 parameters. In this case I just want the sensitivity of the first state wrt all of the parameters.
integrator = casadi.CVodesIntegrator( ... time = numpy.arange(10, dtype=numpy.float) y0 = [1.0, 2.0] param = [3.0, 4.0] adjoint_sensitivities = [] seed = [1, 0] index = 0 for tf in time[1:]: integrator.setOption('t0', 0.0) integrator.setOption('tf', tf) integrator.init() integrator.input(casadi.INTEGRATOR_X0).set(y0) integrator.input(casadi.INTEGRATOR_P).set(param) integrator.setAdjSeed(seed, casadi.INTEGRATOR_XF) integrator.evaluate(0, 1) adjsens = integrator.adjSens(casadi.INTEGRATOR_P).toArray() adjoint_sensitivities.append(adjsens[index])
Thanks!
Joel Andersson
2012-04-18
Hello! Calling integrator.init() in a loop is never efficient since it results in creating new CVodes instances from scratch. Instead of the above, consider scaling time (possibly by introducing a new parameter) so that the time horizon is always . That way you can do away with
integrator.setOption('t0', 0.0) integrator.setOption('tf', tf) integrator.init()
Best,
Joel
bransfor
2012-04-18
This causes a segmentation fault in the windows binary version of CasADi.
Program recieved signal SIGSEGV, Segmentation fault.
0x106f66e in CVAdataStore ()
import numpy import casadi time = casadi.SX('time') y1 = casadi.SX('y1') y2 = casadi.SX('y2') k1 = casadi.SX('k1') k2 = casadi.SX('k2') tf = casadi.SX('tf') right_hand_side = [-k1 * y1, -k2 * y2] scaled_rhs = [eqn * tf for eqn in right_hand_side] rhs_in = casadi.DAE_NUM_IN * [[]] rhs_in[casadi.DAE_T] = [time] rhs_in[casadi.DAE_Y] = casadi.vertcat([y1, y2]) rhs_in[casadi.DAE_P] = casadi.vertcat([k1, k2] + [tf]) rhs = casadi.SXFunction(rhs_in, [scaled_rhs]) rhs.init() # generate integrator instance and set the options (if they exist) integrator = casadi.CVodesIntegrator(rhs) time = numpy.arange(10, dtype=numpy.float) y0 = [1.0, 2.0] param = [3.0, 4.0, 1.0] adjoint_sensitivities = [] seed = [1, 0] index = 0 integrator.setOption('t0', 0.0) integrator.setOption('tf', 1.0) integrator.init() integrator.input(casadi.INTEGRATOR_X0).set(y0) for param[-1] in time[1:]: integrator.input(casadi.INTEGRATOR_P).set(param) integrator.setAdjSeed(seed, casadi.INTEGRATOR_XF) integrator.evaluate(0, 1) adjsens = integrator.adjSens(casadi.INTEGRATOR_P).toArray() adjoint_sensitivities.append(adjsens[index]) print numpy.asfarray(adjoint_sensitivities)
bransfor
2012-04-18
integrator.setOption("steps_per_checkpoint",1000) removes the error and returns the same result as the original method I posted.
Joel Andersson
2012-04-19
Hello,
this is a known issue, see CasADi's FAQ: http://sourceforge.net/apps/trac/casadi/wiki/FAQ. The segfault occurs not in CasADi but in CVodes. I have not yet seen a satisfactory solution to this issue, except for the workaround fix you already mentioned (increasing "steps_per_checkpoint"). There is a new version of Sundials out now, where this issue might or might not have been resolved, though this version is still not compatible with CasADi's Sundials interface.
I'm working on a major change on integrators in CasADi, in which process I will also update the Sundials interface. It should also give another way of addressing this issue. Might take some weeks before it's ready, though.
Joel Andersson
2012-10-08
The segfault issue has now been resolved. Setting "steps_per_checkpoint" should no longer be necessary, although it can still improve performance. Adj. sensitivities of Simulator instances are still not supported.