Changeset 16246 for branches/2925_AutoDiffForDynamicalModels/HeuristicLab.Problems.DynamicalSystemsModelling/3.3/Problem.cs

Timestamp:

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

Author:

gkronber

Message:

#2925: working example with forward sensitivity analysis

File:

• branches/2925_AutoDiffForDynamicalModels/HeuristicLab.Problems.DynamicalSystemsModelling/3.3/Problem.cs (modified) (6 diffs)

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* !
+
 
     /*