Changeset 16246 for branches/2925_AutoDiffForDynamicalModels/AutoDiffForDynamicalModelsTest

Timestamp:

10/19/18 16:05:20 (6 years ago)

Author:

gkronber

Message:

#2925: working example with forward sensitivity analysis

Location:

branches/2925_AutoDiffForDynamicalModels/AutoDiffForDynamicalModelsTest

Files:

• AutoDiffForDynamicalModelsTest.csproj (modified) (1 diff)
• TestCvodes.cs (modified) (4 diffs)

branches/2925_AutoDiffForDynamicalModels/AutoDiffForDynamicalModelsTest/AutoDiffForDynamicalModelsTest.csproj

@@ -67 +67 @@
     </ProjectReference>
     <ProjectReference Include="..\HeuristicLab.Problems.DynamicalSystemsModelling\3.3\HeuristicLab.Problems.DynamicalSystemsModelling-3.3.csproj">
-      <Project>{EDCC8202-4463-4122-B01F-21B664343DFB}</Project>
+      <Project>{edcc8202-4463-4122-b01f-21b664343dfb}</Project>
       <Name>HeuristicLab.Problems.DynamicalSystemsModelling-3.3</Name>
     </ProjectReference>

branches/2925_AutoDiffForDynamicalModels/AutoDiffForDynamicalModelsTest/TestCvodes.cs

@@ -7 +7 @@
   public class TestCvodes {
     [TestMethod]
-    public unsafe void TestCvodesMethod() {
+    public unsafe void TestCvodesASAMethod() {
       // test vectors
       var arr = new double[] { 3.14, 2.71 };
@@ -177 +177 @@
     }
 
-
-    private int JacF(double t, IntPtr y, IntPtr fy, IntPtr Jac, IntPtr user_data, IntPtr tmp1, IntPtr tmp2, IntPtr tmp3) {
-      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;
-    }
+    [TestMethod]
+    public void TestCvodesFSAMethod() {
+
+      // linear oscillator
+      int numberOfEquations = 2;
+      var y = Problem.N_VNew_Serial(numberOfEquations);
+      // y must be initialized before calling CVodeInit
+      // Problem.N_VConst_Serial(100.0, y);
+      Problem.NV_Set_Ith_S(y, 0, 0.5); // x
+      Problem.NV_Set_Ith_S(y, 1, 1);  // v
+
+      var cvode_mem = Problem.CVodeCreate(Problem.MultistepMethod.CV_ADAMS, Problem.NonlinearSolverIteration.CV_FUNCTIONAL);
+
+      var flag = Problem.CVodeInit(cvode_mem, F, 0.0, y);
+      Assert.AreEqual(Problem.CV_SUCCESS, flag);
+
+      flag = Problem.CVodeSStolerances(cvode_mem, 1E-4, 1.0);
+      Assert.AreEqual(Problem.CV_SUCCESS, flag);
+
+      var A = Problem.SUNDenseMatrix(numberOfEquations, numberOfEquations);
+      Assert.AreNotSame(A, IntPtr.Zero);
+      Console.WriteLine(Problem.SUNDenseMatrix_Get(A, 0, 0));
+
+
+      var linearSolver = Problem.SUNDenseLinearSolver(y, A);
+      Assert.AreNotSame(linearSolver, IntPtr.Zero);
+
+      flag = Problem.CVDlsSetLinearSolver(cvode_mem, linearSolver, A);
+      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
+      unsafe {
+        for (int i = 0; i < ns; i++) {
+          var yS0_i = *((IntPtr*)yS0.ToPointer() + i);
+          Problem.N_VConst_Serial(0.0, yS0_i);
+        }
+      }
+
+
+      flag = Problem.CVodeSensInit(cvode_mem, ns, Problem.CV_SIMULTANEOUS, FS, yS0);
+      Assert.AreEqual(Problem.CV_SUCCESS, flag);
+
+      double reltolS = 1E-5;
+      var abstolS = new double[] { 1E-6 };
+      // flag = Problem.CVodeSensSStolerances(cvode_mem, reltolS, abstolS);
+      // Assert.AreEqual(Problem.CV_SUCCESS, flag);
+      flag = Problem.CVodeSensEEtolerances(cvode_mem);
+      Assert.AreEqual(Problem.CV_SUCCESS, flag);
+
+
+      // 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
+      int nout = 1000; // number of output times
+      for (int iout = 0; iout < nout; iout++) {
+        flag = Problem.CVode(cvode_mem, tout, y, ref t, Problem.CV_NORMAL);
+        Assert.AreEqual(Problem.CV_SUCCESS, flag);
+
+        // get sensitivities
+        flag = Problem.CVodeGetSens(cvode_mem, ref t, yS0);
+        Assert.AreEqual(Problem.CV_SUCCESS, flag);
+        unsafe {
+          var ySP0 = *(IntPtr*)yS0.ToPointer(); // parameter of sensitivities for first component
+          Console.WriteLine("{0} {1} {2} {3} {4}", t, Problem.NV_Get_Ith_S(y, 0), Problem.NV_Get_Ith_S(y, 1), Problem.NV_Get_Ith_S(ySP0, 0), Problem.NV_Get_Ith_S(ySP0, 1));
+        }
+
+        tout += 0.1;
+      }
+
+
+
+      Problem.N_VDestroy_Serial(y);
+      Problem.CVodeFree(cvode_mem);
+      Problem.SUNLinSolFree(linearSolver);
+      Problem.SUNMatDestroy(A);
+    }
+
+
 
     public static int F(
@@ -194 +277 @@
       IntPtr user_data) {
 
+      // the system
+      // ∂y1/∂t = y2
+      // ∂y2/∂t = -0.3 y1
       Problem.NV_Set_Ith_S(ydot, 0, Problem.NV_Get_Ith_S(y, 1));
       Problem.NV_Set_Ith_S(ydot, 1, -0.3 * Problem.NV_Get_Ith_S(y, 0));
@@ -199 +285 @@
       ;
     }
+
+    // computes ∂f/∂y or an approximation
+    private int JacF(double t, IntPtr y, IntPtr fy, IntPtr Jac, IntPtr user_data, IntPtr tmp1, IntPtr tmp2, IntPtr tmp3) {
+      // t is current time
+      // y is current y
+      // fy is current f(t,y)
+      // Jac is output
+
+      // Matrix:
+      //  ∂f1/∂y1  ∂f1/∂y2
+      //  ∂f2/∂y1  ∂f2/∂y2
+
+      Problem.SUNDenseMatrix_Set(Jac, 0, 0, 0.0);
+      Problem.SUNDenseMatrix_Set(Jac, 0, 1, 1.0);
+      Problem.SUNDenseMatrix_Set(Jac, 1, 0, -0.3);
+      Problem.SUNDenseMatrix_Set(Jac, 1, 1, 0.0);
+      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;
+    }
+
+    // sensitivity function must compute the vector (∂f/∂y)s_i(t) + ∂f/∂p_i and store in ySdot[i]
+    public static int FS(
+      int Ns,
+      double t,
+      IntPtr y, // N_Vector
+      IntPtr ydot, // N_Vector
+      IntPtr yS, // N_Vector* one vector for each parameter
+      IntPtr ySdot, // N_Vector* one vector for each parameter
+      IntPtr user_data,
+      IntPtr tmp1, // N_Vector
+      IntPtr tmp2 // N_Vector
+      ) {
+
+      // (∂f /∂y)s_i(t) + ∂f /∂p_i
+      var y1 = Problem.NV_Get_Ith_S(y, 0); var y2 = Problem.NV_Get_Ith_S(y, 1);
+
+      unsafe {
+        var yS_p0 = *(IntPtr*)yS.ToPointer();
+        var s1 = Problem.NV_Get_Ith_S(yS_p0, 0); var s2 = Problem.NV_Get_Ith_S(yS_p0, 1);
+
+
+        // (∂f /∂y)s_i(t)  = Jac s_i(t)
+        var sd1 = s2;
+        var sd2 = -0.3 * s1;
+
+        var s_p0 = *(IntPtr*)ySdot.ToPointer();
+        Problem.NV_Set_Ith_S(s_p0, 0, sd1 + 0.0);
+        Problem.NV_Set_Ith_S(s_p0, 1, sd2 - 0.3);
+      }
+
+      return 0;
+    }
+
     public static int FB(
         double t, // realtype