Changeset 16245 for branches/2925_AutoDiffForDynamicalModels/AutoDiffForDynamicalModelsTest

Timestamp:

10/19/18 14:22:01 (6 years ago)

Author:

gkronber

Message:

#2925 worked on integration of CVODES library

File:

• branches/2925_AutoDiffForDynamicalModels/AutoDiffForDynamicalModelsTest/TestCvodes.cs (modified) (7 diffs)

branches/2925_AutoDiffForDynamicalModels/AutoDiffForDynamicalModelsTest/TestCvodes.cs

@@ -21 +21 @@
         int own_data = *(content + 2);
         Console.WriteLine(own_data);
-        double* data = (double*) *(content + 3);
+        double* data = (double*)*(content + 3);
         double data0 = *data;
       }
@@ -40 +40 @@
       var cvode_mem = Problem.CVodeCreate(Problem.MultistepMethod.CV_ADAMS, Problem.NonlinearSolverIteration.CV_FUNCTIONAL);
 
-     var flag = Problem.CVodeInit(cvode_mem, F, 0.0, y);
+      var flag = Problem.CVodeInit(cvode_mem, F, 0.0, y);
       Assert.AreEqual(Problem.CV_SUCCESS, flag);
 
@@ -50 +50 @@
 
 
-      var linearSolver = Problem.SUNDenseLinearSolver(y, A); 
+      var linearSolver = Problem.SUNDenseLinearSolver(y, A);
       Assert.AreNotSame(linearSolver, IntPtr.Zero);
 
@@ -56 +56 @@
       Assert.AreEqual(Problem.CV_SUCCESS, flag);
 
-      flag = Problem.CVDlsSetJacFn(cvode_mem, JacF);
-      Assert.AreEqual(Problem.CV_SUCCESS, flag);
-
-
+      // flag = Problem.CVDlsSetJacFn(cvode_mem, JacF);
+      // Assert.AreEqual(Problem.CV_SUCCESS, flag);
+
+      // var ns = 1; // number of parameters
+      // var p = new double[1]; // set as user-data
+      // var yS0 = Problem.N_VCloneVectorArray_Serial(ns, y); // clone the output vector for each parameter, TODO: free
+      // for(int i=0;i<ns;i++) {
+      // 
+      // }
+
+
+      var q = Problem.N_VNew_Serial(1);
+      Problem.NV_Set_Ith_S(q, 0, 0.0);
+      flag = Problem.CVodeQuadInit(cvode_mem, FQ, q);
+      Assert.AreEqual(Problem.CV_SUCCESS, flag);
+
+
+      int steps = 10;
+      flag = Problem.CVodeAdjInit(cvode_mem, steps, Problem.CV_HERMITE);
+      Assert.AreEqual(Problem.CV_SUCCESS, flag);
+
+      /* step by step forward integration 
       double t = 0.0;
       double tout = 0.1; // first output time
@@ -70 +88 @@
         tout += 0.1; 
       }
+      */
+
+      // complete forward integration
+      double tout = 100.0; // last time
+      double time = 0.0;
+      int ncheck = 0; // number of checkpoints
+      flag = Problem.CVodeF(cvode_mem, tout, y, ref time, Problem.CV_NORMAL, ref ncheck);
+      Assert.AreEqual(Problem.CV_SUCCESS, flag);
+
+      long numSteps = 0;
+      flag = Problem.CVodeGetNumSteps(cvode_mem, ref numSteps);
+      Assert.AreEqual(Problem.CV_SUCCESS, flag);
+
+      var yB = Problem.N_VNew_Serial(numberOfEquations);
+      Problem.N_VConst_Serial(0.0, yB);
+
+      int numberOfParameters = 2;
+      var qB = Problem.N_VNew_Serial(numberOfParameters);
+      Problem.N_VConst_Serial(0.0, qB);
+
+      int indexB = 0;
+      flag = Problem.CVodeCreateB(cvode_mem, Problem.MultistepMethod.CV_BDF, Problem.NonlinearSolverIteration.CV_NEWTON, ref indexB);
+      Assert.AreEqual(Problem.CV_SUCCESS, flag);
+
+      var TB1 = tout;
+      flag = Problem.CVodeInitB(cvode_mem, indexB, FB, TB1, yB);
+      Assert.AreEqual(Problem.CV_SUCCESS, flag);
+
+      var relTolB = 1E-6;
+      var absTolB = 1E-8;
+      flag = Problem.CVodeSStolerancesB(cvode_mem, indexB, relTolB, absTolB);
+      Assert.AreEqual(Problem.CV_SUCCESS, flag);
+
+      var AB = Problem.SUNDenseMatrix(numberOfEquations, numberOfEquations);
+      Assert.AreNotSame(linearSolver, IntPtr.Zero);
+
+      var lsB = Problem.SUNDenseLinearSolver(yB, AB);
+      Assert.AreNotSame(linearSolver, IntPtr.Zero);
+
+      flag = Problem.CVDlsSetLinearSolverB(cvode_mem, indexB, lsB, AB);
+      Assert.AreEqual(Problem.CV_SUCCESS, flag);
+
+      flag = Problem.CVDlsSetJacFnB(cvode_mem, indexB, JacFB);
+      Assert.AreEqual(Problem.CV_SUCCESS, flag);
+
+      flag = Problem.CVodeQuadInitB(cvode_mem, indexB, FQB, qB);
+      Assert.AreEqual(Problem.CV_SUCCESS, flag);
+
+      /* First get results at t = TBout1 */
+
+      /* Call CVodeB to integrate the backward ODE problem. */
+      var tBackOut = 50.0;
+      flag = Problem.CVodeB(cvode_mem, tBackOut, Problem.CV_NORMAL);
+      Assert.AreEqual(Problem.CV_SUCCESS, flag);
+
+      /* Call CVodeGetB to get yB of the backward ODE problem. */
+      flag = Problem.CVodeGetB(cvode_mem, indexB, ref time, yB);
+      Assert.AreEqual(Problem.CV_SUCCESS, flag);
+
+      /* Call CVodeGetAdjY to get the interpolated value of the forward solution
+         y during a backward integration. */
+      flag = Problem.CVodeGetAdjY(cvode_mem, tBackOut, y);
+      Assert.AreEqual(Problem.CV_SUCCESS, flag);
+
+      /* Then at t = T0 */
+
+      double t0 = 0.0;
+      flag = Problem.CVodeB(cvode_mem, t0, Problem.CV_NORMAL);
+      Assert.AreEqual(Problem.CV_SUCCESS, flag);
+
+      long nstB = 0;
+      Problem.CVodeGetNumSteps(Problem.CVodeGetAdjCVodeBmem(cvode_mem, indexB), ref nstB);
+
+      flag = Problem.CVodeGetB(cvode_mem, indexB, ref time, yB);
+      Assert.AreEqual(Problem.CV_SUCCESS, flag);
+
+
+      flag = Problem.CVodeGetQuadB(cvode_mem, indexB, ref time, qB);
+      Assert.AreEqual(Problem.CV_SUCCESS, flag);
+
+      flag = Problem.CVodeGetAdjY(cvode_mem, t0, y);
+      Assert.AreEqual(Problem.CV_SUCCESS, flag);
 
       Problem.N_VDestroy_Serial(y);
@@ -77 +177 @@
     }
 
+
     private int JacF(double t, IntPtr y, IntPtr fy, IntPtr Jac, IntPtr user_data, IntPtr tmp1, IntPtr tmp2, IntPtr tmp3) {
-      // TODO
+      throw new NotImplementedException();
+      return 0;
+    }
+
+    private int JacFB(double t, IntPtr y, IntPtr yB, IntPtr fyB, IntPtr Jac, IntPtr user_data, IntPtr tmp1, IntPtr tmp2, IntPtr tmp3) {
+      throw new NotImplementedException();
       return 0;
     }
@@ -93 +199 @@
       ;
     }
+    public static int FB(
+        double t, // realtype
+        IntPtr y, // N_Vector
+        IntPtr yB, // N_Vector
+        IntPtr yBdot, // N_Vector
+        IntPtr user_data) {
+
+      // y' = f(y,p)
+
+      // yB = λ
+      // λ' = -(∂f / ∂y)^T λ - (∂g / ∂y)
+      // the goal is to find dG/dp where G = integrate( g(y,t,p), t=0, t=T) )
+
+      // for F above:
+      // ∂f / ∂ = 
+      //  0.0 1.0
+      // -0.3 0.0
+
+      // we have no g!?
+      // therefore λ' = 
+      // λ1' = 0.3 λ2 - 0.0
+      // λ2' = 1.0 λ1 - 0.0
+
+      Problem.NV_Set_Ith_S(yBdot, 0, 0.3 * Problem.NV_Get_Ith_S(yB, 1));
+      Problem.NV_Set_Ith_S(yBdot, 1, Problem.NV_Get_Ith_S(yB, 0));
+
+      return 0;
+      ;
+    }
+
+    private int FQ(double t, IntPtr y, IntPtr yQdot, IntPtr user_data) {
+      // TODO: squared error
+      Problem.NV_Set_Ith_S(yQdot, 0, Problem.NV_Get_Ith_S(y, 2));
+      return 0;
+    }
+
+
+    private int FQB(double t, IntPtr y, IntPtr yB, IntPtr qBdot, IntPtr user_data) {
+      throw new NotImplementedException();
+    }
   }
 }