MUQ  0.4.3
RootfindingIVP.cpp
Go to the documentation of this file.
2 
3 // Sundials includes
4 #include <cvodes/cvodes.h>
5 #include <nvector/nvector_serial.h>
6 #include <cvodes/cvodes_dense.h> /* prototype for CVDense */
7 
8 namespace pt = boost::property_tree;
9 using namespace muq::Utilities;
10 using namespace muq::Modeling;
11 
12 RootfindingIVP::RootfindingIVP(std::shared_ptr<WorkPiece> rhs, std::shared_ptr<WorkPiece> root, pt::ptree const& pt, std::shared_ptr<AnyAlgebra> algebra) : ODEBase(rhs, pt, algebra), root(root), maxSteps(pt.get<unsigned int>("Rootfinder.MaxSteps", (int)1e10)), maxTime(pt.get<double>("Rootfinder.MaxTime", 1.0e3)), maxErrorTests(pt.get<unsigned int>("Rootfinder.MaxErrorTests", 100)) {
13  // we must know the number of inputs for both the rhs and the root and they must have at least one (the state)
14  assert(rhs->numInputs>0);
15  assert(root->numInputs>0);
16 
17  // we must know the number of outputs for the root
18  assert(root->numOutputs>0);
19 
20  // set the input and output types
22 }
23 
25 
26 Eigen::VectorXi RootfindingIVP::FindRoot(ref_vector<boost::any> const& inputs, int const wrtIn, int const wrtOut, DerivativeMode const& mode) {
27  // the number of inputs must be at least the the number of inputs required by the rhs and the root
28  assert(inputs.size()>=rhs->numInputs+root->numInputs-1-(autonomous? 0 : 1));
29 
30  // clear the results
31  ClearResults();
32 
33  // create the state vector (have to do a hard copy --- N_Vector is a pointer to the data, the pointer has been declared const, not the data)
34  N_Vector state;
35  DeepCopy(state, boost::any_cast<const N_Vector&>(inputs[0]));
36 
37  // create a data structure to pass around in Sundials
38  auto data = std::make_shared<ODEData>(rhs, root, inputs, autonomous, wrtIn, wrtOut);
39  if( !autonomous ) {
40  const boost::any time = 0.0;
41  data->inputs.insert(data->inputs.begin(), time);
42  }
43 
44  // set solver to null
45  void* cvode_mem = nullptr;
46 
47  // create the solver memory
48  cvode_mem = CVodeCreate(multiStep, solveMethod);
49  assert(CheckFlag((void*)cvode_mem, "CVodeCreate", 0));
50 
51  // initialize the solver
52  CreateSolverMemory(cvode_mem, state, data);
53 
54  // sensState will hold several N_Vectors (because it's a Jacobian)
55  N_Vector *sensState = nullptr;
56  int paramSize = -1;
57 
58  if( wrtIn>=0 && wrtOut==0 && wrtIn<rhs->numInputs ) { // we are computing the derivative wrt one of the rhs parameters
59  paramSize = algebra->Size(inputs[wrtIn]); // the dimension of the parameter
60 
61  // set up sensitivity vector
62  sensState = N_VCloneVectorArray_Serial(paramSize, state);
63  assert(CheckFlag((void *)sensState, "N_VCloneVectorArray_Serial", 0));
64 
65  // initialize the sensitivies to zero
66  for( int is=0; is<paramSize; ++is ) {
67  N_VConst(0.0, sensState[is]);
68  }
69 
70  // set up solver for sensitivity
71  SetUpSensitivity(cvode_mem, paramSize, sensState);
72 
73  // initalize the derivative information
74  InitializeDerivative(1, NV_LENGTH_S(state), paramSize, mode);
75  }
76 
77  // tell the solver how evaluate the root
78  int flag = CVodeRootInit(cvode_mem, root->numOutputs, EvaluateRoot);
79  assert(CheckFlag(&flag, "CVodeRootInit", 1));
80 
81  // set the maximum number of steps
82  flag = CVodeSetMaxNumSteps(cvode_mem, maxSteps);
83  assert(CheckFlag(&flag, "CVodeSetMaxNumSteps", 1));
84 
85  // set the maximum number of error test failures
86  flag = CVodeSetMaxErrTestFails(cvode_mem, maxErrorTests);
87  assert(CheckFlag(&flag, "CVodeSetMaxErrorTestFails", 1));
88 
89  // set the intial time
90  double t = 0.0;
91 
92  // each element corresponds to a vector of desired times, first: the current index of that vector, second: the size of that vector
93  if( inputs.size()>rhs->numInputs+root->numInputs-1 ) {
94  ref_vector<boost::any> outputTimes(inputs.begin()+rhs->numInputs+root->numInputs-1, inputs.end());
95  std::vector<std::pair<unsigned int, unsigned int> > timeIndices = TimeIndices(outputTimes);
96 
97  // first: the next time to integrate to, second: the output index
98  std::pair<double, int> nextTime;
99 
100  while( NextTime(nextTime, timeIndices, outputTimes) ) {
101  // if we have to move forward --- i.e., not at the initial time or another output does not need the current time
102  if( std::fabs(nextTime.first-t)>1.0e-14 ) { // we have to move forward in time
103  flag = CVode(cvode_mem, nextTime.first, state, &t, CV_NORMAL);
104  assert(CheckFlag(&flag, "CVode", 1));
105  }
106 
107  // if we are at the root, but the root time is not a time we asked for do not record!
108  if( std::fabs(nextTime.first-t)<1.0e-14 ) { // we do not have to move forward in time
109  // save the result at this timestep
110  if( timeIndices[nextTime.second].second>1 ) { // the output has more than one compnent ...
111  // .. save the current state as an element in the vector
112  DeepCopy(boost::any_cast<std::vector<N_Vector>&>(outputs[nextTime.second]) [timeIndices[nextTime.second].first-1], state);
113  } else { // the out has one component ...
114  // ... save the current state (not inside a vector)
115  DeepCopy(boost::any_cast<N_Vector&>(outputs[nextTime.second]), state);
116  }
117  }
118 
119  if( flag==CV_ROOT_RETURN ) { break; }
120  }
121  }
122 
123  if( flag!=CV_ROOT_RETURN ) {
124  // integrate forward in time
125  flag = CVode(cvode_mem, maxTime, state, &t, CV_NORMAL);
126  assert(CheckFlag(&flag, "CVode", 1));
127  }
128 
129  // make sure we found a root
130  assert(flag==CV_ROOT_RETURN);
131 
132  // set the output
133  outputs.insert(outputs.begin(), t);
134  outputs.insert(outputs.begin(), state);
135 
136  if( sensState ) {
137  // get the sensitivity
138  flag = CVodeGetSens(cvode_mem, &t, sensState);
139  assert(CheckFlag(&flag, "CVodeGetSens", 1));
140 
141  // create a new jacobian --- shallow copy
142  DlsMat& jac = boost::any_cast<DlsMat&>(*jacobian);
143  for( unsigned int col=0; col<paramSize; ++col ) {
144  DENSE_COL(jac, col) = NV_DATA_S(sensState[col]);
145  NV_DATA_S(sensState[col]) = nullptr;
146  }
147 
148  // delete sensState
149  N_VDestroyVectorArray_Serial(sensState, paramSize);
150  }
151 
152  // retrieve information about the root
153  Eigen::VectorXi rootsfound(root->numOutputs);
154  flag = CVodeGetRootInfo(cvode_mem, rootsfound.data());
155  assert(CheckFlag(&flag, "CVodeGetRootInfo", 1));
156 
157  // free integrator memory (don't destory the state --- outputs is a pointer to it)
158  CVodeFree(&cvode_mem);
159 
160  return rootsfound;
161 }
162 
164  // find the first root to one of the root outputs
165  FindRoot(inputs);
166 }
167 
168 void RootfindingIVP::JacobianImpl(unsigned int const wrtIn, unsigned int const wrtOut, ref_vector<boost::any> const& inputs) {
169  if( wrtOut>0 ) {
170  std::cerr << std::endl << "ERROR: Cannot compute derivative infomration of muq::Modeling::RootfindingIVP for outputs other than the state at the root (the first output)" << std::endl << std::endl;
171 
172  assert(wrtOut==0);
173  }
174 
175  // find the first root to one of the root outputs
176  const Eigen::VectorXi& rootsfound = FindRoot(inputs, wrtIn, wrtOut);
177  unsigned int ind = 0;
178  for( ; ind<rootsfound.size(); ++ind ) {
179  if( rootsfound(ind)!=0 ) { break; }
180  }
181 
182  // the inputs to the rhs function
183  ref_vector<boost::any> rhsIns;
184  rhsIns.push_back(outputs[0]);
185  rhsIns.insert(rhsIns.end(), inputs.begin()+1, inputs.begin()+rhs->numInputs);
186 
187  // the inputs to the root function
188  ref_vector<boost::any> rootIns;
189  rootIns.push_back(outputs[0]);
190  rootIns.insert(rootIns.end(), inputs.begin()+rhs->numInputs, inputs.begin()+rhs->numInputs+root->numInputs-1);
191 
192  // the derivative of the root function at the final time -- derivative of the root wrt the state
193  const boost::any& dgdy = root->Jacobian(0, ind, rootIns);
194  const DlsMat& dgdyref = boost::any_cast<const DlsMat&>(dgdy); // 1 x stateSize matrix
195 
196  // evaluate the right hand side
197  const std::vector<boost::any>& f = rhs->Evaluate(rhsIns); // stateSize x 1 vector
198  const N_Vector& fref = boost::any_cast<N_Vector>(f[0]);
199  assert(NV_LENGTH_S(fref)==dgdyref->N);
200 
201  double dum = 0.0;
202  for( unsigned int i=0; i<dgdyref->N; ++i ) {
203  dum += DENSE_ELEM(dgdyref, 0, i)*NV_Ith_S(fref, i);
204  }
205 
206  if( wrtIn>=rhs->numInputs ) { // derivative wrt root parameter
207  boost::any dgdpara = root->Jacobian(wrtIn-rhs->numInputs+1, ind, rootIns);
208  const DlsMat& dgdpararef = boost::any_cast<const DlsMat&>(dgdpara); // 1 x paramSize matrix
209  DenseScale(-1.0/dum, dgdpararef);
210 
211  jacobian = NewDenseMat(NV_LENGTH_S(fref), dgdpararef->N);
212  DlsMat& jac = boost::any_cast<DlsMat&>(*jacobian); // stateSize x paramSize matrix
213  SetToZero(jac);
214 
215  // update the jacobian jac += f.transpose()*dgdpara (outer product)
216  for( unsigned int i=0; i<jac->M; ++i ) {
217  for( unsigned int j=0; j<jac->N; ++j ) {
218  DENSE_ELEM(jac, i, j) = NV_Ith_S(fref, i)*DENSE_ELEM(dgdpararef, 0, j);
219  }
220  }
221 
222  return;
223  }
224 
225  // get a reference to the jacobian
226  DlsMat& jac = boost::any_cast<DlsMat&>(*jacobian); // stateSize x paramSize matrix
227  assert(jac->M==dgdyref->N);
228  assert(NV_LENGTH_S(fref)==jac->M);
229 
230  if( wrtIn==0 ) { // derivative wrt the initial conditions
231  // add the identity
232  AddIdentity(jac);
233  }
234 
235  DenseScale(-1.0/dum, dgdyref);
236 
237  // the derivative of the final time wrt to the parameter; dtfdpara = -jac_{state}(root)*jac/dot(jac_{state}(root), f)
238  N_Vector dtfdic = N_VNew_Serial(jac->N);
239  for( unsigned int i=0; i<jac->N; ++i ) {
240  NV_Ith_S(dtfdic, i) = 0.0;
241  for( unsigned int j=0; j<jac->M; ++j ) {
242  NV_Ith_S(dtfdic, i) += DENSE_ELEM(dgdyref, 0, j)*DENSE_ELEM(jac, j, i);
243  }
244  }
245 
246  // update the jacobian jac += f.transpose()*dtfdpara (outer product)
247  for( unsigned int i=0; i<jac->M; ++i ) {
248  for( unsigned int j=0; j<jac->N; ++j ) {
249  DENSE_ELEM(jac, i, j) += NV_Ith_S(fref, i)*NV_Ith_S(dtfdic, j);
250  }
251  }
252 
253  // destroy temp vectors/matrices
254  N_VDestroy(dtfdic);
255 }
256 
258  if( autonomous ) {
259  // the type of the first input (the state) for the rhs and the root
260  assert(rhs->InputType(0, false).compare(root->InputType(0, false))==0);
261  } else { // non-autonomous
262  // the input type of the second input of the rhs is the input of the first input of the root
263  assert(rhs->InputType(1, false).compare(root->InputType(0, false))==0);
264  }
265 
266  // the second set of input parameters are the parameters for the root
267  for( auto intype : root->InputTypes() ) {
268  if( intype.first==0 ) { // we've already set the state type
269  continue;
270  }
271 
272  inputTypes[rhs->numInputs+intype.first-1] = intype.second;
273  }
274 
275  // the first output type is the state
276  outputTypes[0] = root->InputType(0, false);
277 
278  // the second output is the time where the root is reached
279  outputTypes[1] = typeid(double).name();
280 }
281 
282 int RootfindingIVP::EvaluateRoot(realtype t, N_Vector state, realtype *root, void *user_data) {
283  // get the data type
284  ODEData* data = (ODEData*)user_data;
285  assert(data);
286  assert(data->root);
287 
288  // set the state input
289  const boost::any& anyref = state;
290 
291  // the inputs the root function
292  ref_vector<boost::any> rootins(data->inputs.begin()+data->rhs->numInputs, data->inputs.begin()+data->rhs->numInputs+data->root->numInputs-1);
293  rootins.insert(rootins.begin(), anyref);
294 
295  // evaluate the root
296  const std::vector<boost::any>& result = data->root->Evaluate(rootins);
297  assert(result.size()==data->root->numOutputs);
298  for( unsigned int i=0; i<result.size(); ++i ) {
299  root[i] = boost::any_cast<double>(result[i]);
300  }
301 
302  return 0;
303 }
A helper class for Sundial's time integration.
Definition: ODEData.h:18
std::shared_ptr< ModPiece > root
A function we are trying to find the root of along an orbit of the ODE (nullptr if we are not doing a...
Definition: ODEData.h:59
ref_vector< Eigen::VectorXd > inputs
Definition: ODEData.h:64
std::shared_ptr< ModPiece > rhs
The right hand side of the ODE.
Definition: ODEData.h:56
const unsigned int maxSteps
The maximum number of steps the timesteper can take.
std::shared_ptr< WorkPiece > root
The root function.
virtual void JacobianImpl(unsigned int const wrtIn, unsigned int const wrtOut, ref_vector< boost::any > const &inputs) override
Compute the Jacobian of the state at the root.
const unsigned int maxErrorTests
The maximum number of error test failures.
Eigen::VectorXi FindRoot(ref_vector< boost::any > const &inputs, int const wrtIn=-1, int const wrtOut=-1, DerivativeMode const &mode=DerivativeMode::Jac)
Run the CVODES integrator.
const double maxTime
The maximum amount of time to integrate.
virtual void EvaluateImpl(ref_vector< boost::any > const &inputs) override
Integrate the ODE until we find a root.
static int EvaluateRoot(realtype t, N_Vector state, realtype *root, void *user_data)
Evaluate the root function.
void UpdateInputOutputTypes()
Update the input and output types based on the rhs and root muq::Modeling::WorkPiece's.
std::vector< std::reference_wrapper< const T > > ref_vector
A vector of references to something ...
Definition: WorkPiece.h:37
auto get(const nlohmann::detail::iteration_proxy_value< IteratorType > &i) -> decltype(i.key())
Definition: json.h:3956