Changeset 16246 for branches/2925_AutoDiffForDynamicalModels

Timestamp:

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

Author:

gkronber

Message:

#2925: working example with forward sensitivity analysis

Location:

branches/2925_AutoDiffForDynamicalModels

Files:

• AutoDiffForDynamicalModelsTest/AutoDiffForDynamicalModelsTest.csproj (modified) (1 diff)
• AutoDiffForDynamicalModelsTest/TestCvodes.cs (modified) (4 diffs)
• HeuristicLab.Problems.DynamicalSystemsModelling/3.3/Problem.cs (modified) (6 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

branches/2925_AutoDiffForDynamicalModels/HeuristicLab.Problems.DynamicalSystemsModelling/3.3/Problem.cs

@@ -223 +223 @@
     public const int CV_REIFWD_FAIL = -105;
     public const int CV_FWD_FAIL = -106;
+
+
     public const int CV_GETY_BADT = -107;
 
@@ -284 +286 @@
       );
 
+
+    // to calculate sensitivities RHS for all equations at once
+    // must compute (∂f/∂y)s_i(t) + ∂f/∂p_i
+    [UnmanagedFunctionPointer(CallingConvention.Cdecl)]
+    public delegate int CVSensRhsFn(
+      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
+    );
+
     [DllImport("sundials_cvodes.dll", EntryPoint = "CVodeCreate", ExactSpelling = true, CallingConvention = CallingConvention.Cdecl)]
     // returns a void* to the cvodes memory block if successful otherwise NULL
@@ -289 +307 @@
 
     [DllImport("sundials_cvodes.dll", EntryPoint = "CVodeInit", ExactSpelling = true, CallingConvention = CallingConvention.Cdecl)]
-
     public static extern int CVodeInit(
       IntPtr cvode_mem, // pointer returned by CVodeCreate
@@ -316 +333 @@
       IntPtr cvode_mem,
       CVDlsJacFunc jacFunc
+      );
+
+    [DllImport("sundials_cvodes.dll", EntryPoint = "CVodeSensInit", ExactSpelling = true, CallingConvention = CallingConvention.Cdecl)]
+    public static extern int CVodeSensInit(
+      IntPtr cvode_mem,
+      int Ns, // number of parameters
+      int ism, // sensitivity solution method, CV_SIMULTANEOUS or CV_STAGGERED
+      CVSensRhsFn fS, // right hand side function which computes all sensitivity RHS at the same time
+      IntPtr yS0 // N_Vector
+      );
+
+    [DllImport("sundials_cvodes.dll", EntryPoint = "CVodeSensSStolerances", ExactSpelling = true, CallingConvention = CallingConvention.Cdecl)]
+    public static extern int CVodeSensSStolerances(
+      IntPtr cvode_mem,
+      double reltol,
+      double[] abstol
+    );
+
+    /* Call CVodeSensEEtolerances to estimate tolerances for sensitivity 
+   variables based on the rolerances supplied for states variables and 
+   the scaling factor pbar */
+    [DllImport("sundials_cvodes.dll", EntryPoint = "CVodeSensEEtolerances", ExactSpelling = true, CallingConvention = CallingConvention.Cdecl)]
+    public static extern int CVodeSensEEtolerances(
+      IntPtr cvode_mem
+      );
+
+    [DllImport("sundials_cvodes.dll", EntryPoint = "CVodeGetSens", ExactSpelling = true, CallingConvention = CallingConvention.Cdecl)]
+    public static extern int CVodeGetSens(
+      IntPtr cvode_mem,
+      ref double tret,
+      IntPtr yS //N_Vector*, one vector for each parameter
       );
 
@@ -452 +500 @@
     [DllImport("sundials_cvodes.dll", EntryPoint = "SUNMatDestroy", ExactSpelling = true, CallingConvention = CallingConvention.Cdecl)]
     public static extern void SUNMatDestroy(IntPtr A);
+
+    [DllImport("sundials_cvodes.dll", EntryPoint = "SUNDenseMatrix_Data", ExactSpelling = true, CallingConvention = CallingConvention.Cdecl)]
+    public unsafe static extern double* SUNDenseMatrix_Data(IntPtr matrix);
+
+    [DllImport("sundials_cvodes.dll", EntryPoint = "SUNDenseMatrix_Cols", ExactSpelling = true, CallingConvention = CallingConvention.Cdecl)]
+    public static extern long SUNDenseMatrix_Cols(IntPtr matrix);
+
+    [DllImport("sundials_cvodes.dll", EntryPoint = "SUNDenseMatrix_Rows", ExactSpelling = true, CallingConvention = CallingConvention.Cdecl)]
+    public static extern long SUNDenseMatrix_Rows(IntPtr matrix);
+
+
+    public unsafe static double SUNDenseMatrix_Get(IntPtr mat, long i, long j) {
+      long M = SUNDenseMatrix_Rows(mat);
+      // the(i, j)th element of A(with 0 <= i<M and 0 <= j<N) is given by(A->data)[j*M+i]
+      return SUNDenseMatrix_Data(mat)[j*M+i];
+    }
+
+    public unsafe static void SUNDenseMatrix_Set(IntPtr mat, long i, long j, double val) {
+      long M = SUNDenseMatrix_Rows(mat);
+      // the(i, j)th element of A(with 0 <= i<M and 0 <= j<N) is given by(A->data)[j*M+i]
+      SUNDenseMatrix_Data(mat)[j * M + i] = val;
+    }
+
     #endregion
 
@@ -515 +586 @@
 
     [DllImport("sundials_cvodes.dll", EntryPoint = "N_VCloneVectorArray_Serial", ExactSpelling = true, CallingConvention = CallingConvention.Cdecl)]
-    public static extern IntPtr N_VCloneVectorArray_Serial(int count, IntPtr vec);
+    public static extern IntPtr N_VCloneVectorArray_Serial(int count, IntPtr vec); // returns N_Vector* !
+
 
     /*