@@ -64,6 +64,40 @@ def __init__(self, base_variable, solution, known_evals=None, warn=True):
6464 self ._sensitivity = None
6565 self .solution_sensitivity = solution .sensitivity
6666
67+ # Special case: symbolic solution, with casadi
68+ if isinstance (solution .y , casadi .Function ):
69+ # Evaluate solution at specific inputs value
70+ inputs_stacked = casadi .vertcat (* solution .inputs .values ())
71+ self .u_sol = solution .y (inputs_stacked ).full ()
72+ # Convert variable to casadi
73+ t_MX = casadi .MX .sym ("t" )
74+ y_MX = casadi .MX .sym ("y" , self .u_sol .shape [0 ])
75+ # Make all inputs symbolic first for converting to casadi
76+ symbolic_inputs_dict = {
77+ name : casadi .MX .sym (name , value .shape [0 ])
78+ for name , value in solution .inputs .items ()
79+ }
80+
81+ # The symbolic_inputs will be used for sensitivity
82+ symbolic_inputs = casadi .vertcat (* symbolic_inputs_dict .values ())
83+ try :
84+ var_casadi = base_variable .to_casadi (
85+ t_MX , y_MX , inputs = symbolic_inputs_dict
86+ )
87+ except :
88+ n = 1
89+ self .base_variable_sym = casadi .Function (
90+ "variable" , [t_MX , y_MX , symbolic_inputs ], [var_casadi ]
91+ )
92+ # Store symbolic inputs for sensitivity
93+ self .symbolic_inputs = symbolic_inputs
94+ self .y_sym = solution .y (symbolic_inputs )
95+ else :
96+ self .u_sol = solution .y
97+ self .base_variable_sym = None
98+ self .symbolic_inputs = None
99+ self .y_sym = None
100+
67101 # Set timescale
68102 self .timescale = solution .model .timescale .evaluate ()
69103 self .t_pts = self .t_sol * self .timescale
@@ -78,8 +112,8 @@ def __init__(self, base_variable, solution, known_evals=None, warn=True):
78112 # Evaluate base variable at initial time
79113 if self .known_evals :
80114 self .base_eval , self .known_evals [solution .t [0 ]] = base_variable .evaluate (
81- solution . t [0 ],
82- solution . y [:, 0 ],
115+ self . t_sol [0 ],
116+ self . u_sol [:, 0 ],
83117 inputs = {name : inp [:, 0 ] for name , inp in solution .inputs .items ()},
84118 known_evals = self .known_evals [solution .t [0 ]],
85119 )
@@ -571,10 +605,20 @@ def sensitivity(self):
571605 return {}
572606 # Otherwise initialise and return sensitivity
573607 if self ._sensitivity is None :
574- self .initialise_sensitivity ()
608+ # Check that we can compute sensitivities
609+ if self .base_variable_sym is None and self .solution_sensitivity == {}:
610+ raise ValueError (
611+ "Cannot compute sensitivities. The 'sensitivity' argument of the "
612+ "solver should be changed from 'None' to allow sensitivity "
613+ "calculations. Check solver documentation for details."
614+ )
615+ if self .base_variable_sym is None :
616+ self .initialise_sensitivity_explicit_forward ()
617+ else :
618+ self .initialise_sensitivity_casadi ()
575619 return self ._sensitivity
576620
577- def initialise_sensitivity (self ):
621+ def initialise_sensitivity_explicit_forward (self ):
578622 "Set up the sensitivity dictionary"
579623 inputs_stacked = casadi .vertcat (* [p for p in self .inputs .values ()])
580624
@@ -628,6 +672,79 @@ def initialise_sensitivity(self):
628672 # Save attribute
629673 self ._sensitivity = sensitivity
630674
675+ def initialise_sensitivity_casadi (self ):
676+ def initialise_0D_symbolic ():
677+ "Create a 0D symbolic variable"
678+ # Evaluate the base_variable index-by-index
679+ for idx in range (len (self .t_sol )):
680+ t = self .t_sol [idx ]
681+ u = self .y_sym [:, idx ]
682+ next_entries = self .base_variable_sym (t , u , self .symbolic_inputs )
683+ if idx == 0 :
684+ entries = next_entries
685+ else :
686+ entries = casadi .horzcat (entries , next_entries )
687+
688+ return entries
689+
690+ def initialise_1D_symbolic ():
691+ "Create a 1D symbolic variable"
692+ # Evaluate the base_variable index-by-index
693+ for idx in range (len (self .t_sol )):
694+ t = self .t_sol [idx ]
695+ u = self .y_sym [:, idx ]
696+ next_entries = self .base_variable_sym (t , u , self .symbolic_inputs )
697+ if idx == 0 :
698+ entries = next_entries
699+ else :
700+ entries = casadi .vertcat (entries , next_entries )
701+
702+ return entries
703+
704+ inputs_stacked = casadi .vertcat (* self .inputs .values ())
705+ self .base_eval = self .base_variable_sym (
706+ self .t_sol [0 ], self .u_sol [:, 0 ], inputs_stacked
707+ )
708+ if (
709+ isinstance (self .base_eval , numbers .Number )
710+ or len (self .base_eval .shape ) == 0
711+ or self .base_eval .shape [0 ] == 1
712+ ):
713+ entries_MX = initialise_0D_symbolic ()
714+ else :
715+ n = self .mesh .npts
716+ base_shape = self .base_eval .shape [0 ]
717+ # Try shape that could make the variable a 1D variable
718+ if base_shape == n :
719+ entries_MX = initialise_1D_symbolic ()
720+ else :
721+ # Raise error for 2D variable
722+ raise NotImplementedError (
723+ "Shape not recognized for {} " .format (self .base_variable )
724+ + "(note processing of 2D and 3D variables is not yet "
725+ + "implemented)"
726+ )
727+
728+ # Make entries a function and compute jacobian
729+ casadi_entries_fn = casadi .Function (
730+ "variable" , [self .symbolic_inputs ], [entries_MX ]
731+ )
732+
733+ sens_MX = casadi .jacobian (entries_MX , self .symbolic_inputs )
734+ casadi_sens_fn = casadi .Function ("variable" , [self .symbolic_inputs ], [sens_MX ])
735+
736+ sens_eval = casadi_sens_fn (inputs_stacked )
737+ sensitivity = {"all" : sens_eval }
738+
739+ # Add the individual sensitivity
740+ start = 0
741+ for name , inp in self .inputs .items ():
742+ end = start + inp .shape [0 ]
743+ sensitivity [name ] = sens_eval [:, start :end ]
744+ start = end
745+
746+ self ._sensitivity = sensitivity
747+
631748
632749def eval_dimension_name (name , x , r , y , z ):
633750 if name == "x" :
0 commit comments