Changeset 16245

Timestamp:

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

Author:

gkronber

Message:

#2925 worked on integration of CVODES library

Location:

branches/2925_AutoDiffForDynamicalModels

Files:

• AutoDiffForDynamicalModelsTest/TestCvodes.cs (modified) (7 diffs)
• HeuristicLab.Problems.DynamicalSystemsModelling/3.3/Problem.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();
+    }
   }
 }

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

@@ -124 +124 @@
   }
 
+
+  // Eine weitere Möglichkeit ist spline-smoothing der Daten (über Zeit) um damit für jede Zielvariable
+  // einen bereinigten (Rauschen) Wert und die Ableitung dy/dt für alle Beobachtungen zu bekommen
+  // danach kann man die Regression direkt für dy/dt anwenden (mit bereinigten y als inputs)
   [Item("Dynamical Systems Modelling Problem", "TODO")]
   [Creatable(CreatableAttribute.Categories.GeneticProgrammingProblems, Priority = 900)]
@@ -226 +230 @@
         IntPtr y, // N_Vector
         IntPtr ydot, // N_Vector
+        IntPtr user_data
+      );
+
+    [UnmanagedFunctionPointer(CallingConvention.Cdecl)]
+    public delegate int CVRhsFuncB(
+        double t, // realtype
+        IntPtr y, // N_Vector
+        IntPtr yB, // N_Vector
+        IntPtr yBdot, // N_Vector
         IntPtr user_data
       );
@@ -241 +254 @@
       );
 
+    [UnmanagedFunctionPointer(CallingConvention.Cdecl)]
+    public delegate int CVDlsJacFuncB(
+      double t,
+      IntPtr y, // N_Vector
+      IntPtr yB, // N_Vector
+      IntPtr fyB, // N_Vector
+      IntPtr Jac, // SUNMatrix
+      IntPtr user_data,
+      IntPtr tmp1, // N_Vector
+      IntPtr tmp2, // N_Vector
+      IntPtr tmp3 // N_Vector
+    );
+
+    [UnmanagedFunctionPointer(CallingConvention.Cdecl)]
+    public delegate int CVQuadRhsFn(
+      double t,
+      IntPtr y, // N_Vector
+      IntPtr yQdot, // N_Vector
+      IntPtr user_data
+      );
+
+    [UnmanagedFunctionPointer(CallingConvention.Cdecl)]
+    public delegate int CVQuadRhsFnB(
+      double t,
+      IntPtr y, // N_Vector
+      IntPtr yB, // N_Vector
+      IntPtr qBdot, // N_Vector
+      IntPtr user_data
+      );
 
     [DllImport("sundials_cvodes.dll", EntryPoint = "CVodeCreate", ExactSpelling = true, CallingConvention = CallingConvention.Cdecl)]
@@ -274 +316 @@
       IntPtr cvode_mem,
       CVDlsJacFunc jacFunc
+      );
+
+    [DllImport("sundials_cvodes.dll", EntryPoint = "CVodeQuadInit", ExactSpelling = true, CallingConvention = CallingConvention.Cdecl)]
+    public static extern int CVodeQuadInit(
+      IntPtr cvode_mem,
+      CVQuadRhsFn qF,
+      IntPtr yQ0 // N_Vector, initial value of yQ
+      );
+
+    [DllImport("sundials_cvodes.dll", EntryPoint = "CVodeQuadInitB", ExactSpelling = true, CallingConvention = CallingConvention.Cdecl)]
+    public static extern int CVodeQuadInitB(
+      IntPtr cvode_mem,
+      int indexB,
+      CVQuadRhsFnB rhsQB,
+      IntPtr yQB0 // N_Vector, initial value of yQB
+    );
+
+
+    [DllImport("sundials_cvodes.dll", EntryPoint = "CVodeAdjInit", ExactSpelling = true, CallingConvention = CallingConvention.Cdecl)]
+    public static extern int CVodeAdjInit(
+      IntPtr cvode_mem,
+      int nd,  // integration steps between checkpoints
+      int interpType // either CV_POLYNOMIAL or CV_HERMITE
+      );
+
+    [DllImport("sundials_cvodes.dll", EntryPoint = "CVodeF", ExactSpelling = true, CallingConvention = CallingConvention.Cdecl)]
+    public static extern int CVodeF(
+        IntPtr cvode_mem,
+        double t_out, // the next time at which a computed solution is desired
+        IntPtr y, // N_Vector, the solution vector y
+        ref double t_ret, // the time reached by the solver (output)
+        int itask, // CV_NORMAL or CV_ONE_STEP
+        ref int ncheck  // the number of internal checkpoints stored so far.
+      );
+
+    [DllImport("sundials_cvodes.dll", EntryPoint = "CVodeGetNumSteps", ExactSpelling = true, CallingConvention = CallingConvention.Cdecl)]
+    public static extern int CVodeGetNumSteps(
+      IntPtr cvode_mem,
+      ref long ncheck  // the number of steps taken
+    );
+
+    [DllImport("sundials_cvodes.dll", EntryPoint = "CVodeCreateB", ExactSpelling = true, CallingConvention = CallingConvention.Cdecl)]
+    public static extern int CVodeCreateB(
+        IntPtr cvode_mem,
+        MultistepMethod lmmB,
+        NonlinearSolverIteration iterB,
+        ref int which
+      );
+
+    [DllImport("sundials_cvodes.dll", EntryPoint = "CVodeInitB", ExactSpelling = true, CallingConvention = CallingConvention.Cdecl)]
+    public static extern int CVodeInitB(
+        IntPtr cvode_mem,
+        int which,
+        CVRhsFuncB rightHandSide,
+        double tB0, // endpoint T where final conditions are provided (equal to endpoint of forward integration)
+        IntPtr yB0 // N_Vector inital value at t=tb0 of backward solution
       );
 
@@ -285 +383 @@
       );
 
+
+    [DllImport("sundials_cvodes.dll", EntryPoint = "CVodeSStolerancesB", ExactSpelling = true, CallingConvention = CallingConvention.Cdecl)]
+    public static extern int CVodeSStolerancesB(
+      IntPtr cvode_mem,
+      int which,
+      double reltol,
+      double abstol
+    );
+
+    [DllImport("sundials_cvodes.dll", EntryPoint = "CVDlsSetLinearSolverB", ExactSpelling = true, CallingConvention = CallingConvention.Cdecl)]
+    public static extern int CVDlsSetLinearSolverB(
+      IntPtr cvode_mem,
+      int which,
+      IntPtr linearSolver, // SUNLinearSolver
+      IntPtr j // SUNMatrix
+    );
+
+    [DllImport("sundials_cvodes.dll", EntryPoint = "CVDlsSetJacFnB", ExactSpelling = true, CallingConvention = CallingConvention.Cdecl)]
+    public static extern int CVDlsSetJacFnB(
+      IntPtr cvode_mem,
+      int which,
+      CVDlsJacFuncB jacFunc
+    );
+
+    [DllImport("sundials_cvodes.dll", EntryPoint = "CVodeB", ExactSpelling = true, CallingConvention = CallingConvention.Cdecl)]
+    public static extern int CVodeB(
+      IntPtr cvode_mem,
+      double tout, // next time at which a solution is desired
+      int itask // flag indicating the job of the solver for the next step. 
+    );
+
+    [DllImport("sundials_cvodes.dll", EntryPoint = "CVodeGetB", ExactSpelling = true, CallingConvention = CallingConvention.Cdecl)]
+    public static extern int CVodeGetB(
+          IntPtr cvode_mem,
+          int which,
+          ref double tret,
+          IntPtr yB // the backward solution at time tret
+        );
+    [DllImport("sundials_cvodes.dll", EntryPoint = "CVodeGetAdjY", ExactSpelling = true, CallingConvention = CallingConvention.Cdecl)]
+    public static extern int CVodeGetAdjY(
+          IntPtr cvode_mem,
+          double t,
+          IntPtr y // the forward solution y(t)
+        );
+
+    [DllImport("sundials_cvodes.dll", EntryPoint = "CVodeGetAdjCVodeBmem", ExactSpelling = true, CallingConvention = CallingConvention.Cdecl)]
+    public static extern IntPtr CVodeGetAdjCVodeBmem(
+              IntPtr cvode_mem,
+              int which
+            );
+
+    [DllImport("sundials_cvodes.dll", EntryPoint = "CVodeGetQuadB", ExactSpelling = true, CallingConvention = CallingConvention.Cdecl)]
+    public static extern int CVodeGetQuadB(
+        IntPtr cvode_mem,
+        int which,
+        ref double tret,
+        IntPtr yQB // N_Vector
+      );
+
     [DllImport("sundials_cvodes.dll", EntryPoint = "CVodeFree", ExactSpelling = true, CallingConvention = CallingConvention.Cdecl)]
 
@@ -357 +514 @@
     public unsafe static extern double* N_VGetArrayPointer_Serial(IntPtr vec);
 
+    [DllImport("sundials_cvodes.dll", EntryPoint = "N_VCloneVectorArray_Serial", ExactSpelling = true, CallingConvention = CallingConvention.Cdecl)]
+    public static extern IntPtr N_VCloneVectorArray_Serial(int count, IntPtr vec);
+
     /*
 #define NV_CONTENT_S(v)  ( (N_VectorContent_Serial)(v->content) )
@@ -368 +528 @@
 #define NV_Ith_S(v,i)    ( NV_DATA_S(v)[i] )
 */
-// methods for macros
+    // methods for macros
     public unsafe static int* NV_CONTENT_S(IntPtr v) {
       int* content = (int*)*(int*)v.ToPointer();