Learn how easy it is to sync an existing GitHub or Google Code repo to a SourceForge project! See Demo

Close

Substitute for Adj. sensitivity in Simulator

bransfor
2012-04-18
2013-04-11
  • bransfor
    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
    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
    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
    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
    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
    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.