Free cookie consent management tool by TermsFeed Policy Generator

Changeset 16246


Ignore:
Timestamp:
10/19/18 16:05:20 (4 years ago)
Author:
gkronber
Message:

#2925: working example with forward sensitivity analysis

Location:
branches/2925_AutoDiffForDynamicalModels
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • branches/2925_AutoDiffForDynamicalModels/AutoDiffForDynamicalModelsTest/AutoDiffForDynamicalModelsTest.csproj

    r16225 r16246  
    6767    </ProjectReference>
    6868    <ProjectReference Include="..\HeuristicLab.Problems.DynamicalSystemsModelling\3.3\HeuristicLab.Problems.DynamicalSystemsModelling-3.3.csproj">
    69       <Project>{EDCC8202-4463-4122-B01F-21B664343DFB}</Project>
     69      <Project>{edcc8202-4463-4122-b01f-21b664343dfb}</Project>
    7070      <Name>HeuristicLab.Problems.DynamicalSystemsModelling-3.3</Name>
    7171    </ProjectReference>
  • branches/2925_AutoDiffForDynamicalModels/AutoDiffForDynamicalModelsTest/TestCvodes.cs

    r16245 r16246  
    77  public class TestCvodes {
    88    [TestMethod]
    9     public unsafe void TestCvodesMethod() {
     9    public unsafe void TestCvodesASAMethod() {
    1010      // test vectors
    1111      var arr = new double[] { 3.14, 2.71 };
     
    177177    }
    178178
    179 
    180     private int JacF(double t, IntPtr y, IntPtr fy, IntPtr Jac, IntPtr user_data, IntPtr tmp1, IntPtr tmp2, IntPtr tmp3) {
    181       throw new NotImplementedException();
    182       return 0;
    183     }
    184 
    185     private int JacFB(double t, IntPtr y, IntPtr yB, IntPtr fyB, IntPtr Jac, IntPtr user_data, IntPtr tmp1, IntPtr tmp2, IntPtr tmp3) {
    186       throw new NotImplementedException();
    187       return 0;
    188     }
     179    [TestMethod]
     180    public void TestCvodesFSAMethod() {
     181
     182      // linear oscillator
     183      int numberOfEquations = 2;
     184      var y = Problem.N_VNew_Serial(numberOfEquations);
     185      // y must be initialized before calling CVodeInit
     186      // Problem.N_VConst_Serial(100.0, y);
     187      Problem.NV_Set_Ith_S(y, 0, 0.5); // x
     188      Problem.NV_Set_Ith_S(y, 1, 1);  // v
     189
     190      var cvode_mem = Problem.CVodeCreate(Problem.MultistepMethod.CV_ADAMS, Problem.NonlinearSolverIteration.CV_FUNCTIONAL);
     191
     192      var flag = Problem.CVodeInit(cvode_mem, F, 0.0, y);
     193      Assert.AreEqual(Problem.CV_SUCCESS, flag);
     194
     195      flag = Problem.CVodeSStolerances(cvode_mem, 1E-4, 1.0);
     196      Assert.AreEqual(Problem.CV_SUCCESS, flag);
     197
     198      var A = Problem.SUNDenseMatrix(numberOfEquations, numberOfEquations);
     199      Assert.AreNotSame(A, IntPtr.Zero);
     200      Console.WriteLine(Problem.SUNDenseMatrix_Get(A, 0, 0));
     201
     202
     203      var linearSolver = Problem.SUNDenseLinearSolver(y, A);
     204      Assert.AreNotSame(linearSolver, IntPtr.Zero);
     205
     206      flag = Problem.CVDlsSetLinearSolver(cvode_mem, linearSolver, A);
     207      Assert.AreEqual(Problem.CV_SUCCESS, flag);
     208
     209      flag = Problem.CVDlsSetJacFn(cvode_mem, JacF);
     210      Assert.AreEqual(Problem.CV_SUCCESS, flag);
     211
     212      var ns = 1; // number of parameters
     213      var p = new double[1]; // set as user-data
     214      var yS0 = Problem.N_VCloneVectorArray_Serial(ns, y); // clone the output vector for each parameter, TODO: free
     215      unsafe {
     216        for (int i = 0; i < ns; i++) {
     217          var yS0_i = *((IntPtr*)yS0.ToPointer() + i);
     218          Problem.N_VConst_Serial(0.0, yS0_i);
     219        }
     220      }
     221
     222
     223      flag = Problem.CVodeSensInit(cvode_mem, ns, Problem.CV_SIMULTANEOUS, FS, yS0);
     224      Assert.AreEqual(Problem.CV_SUCCESS, flag);
     225
     226      double reltolS = 1E-5;
     227      var abstolS = new double[] { 1E-6 };
     228      // flag = Problem.CVodeSensSStolerances(cvode_mem, reltolS, abstolS);
     229      // Assert.AreEqual(Problem.CV_SUCCESS, flag);
     230      flag = Problem.CVodeSensEEtolerances(cvode_mem);
     231      Assert.AreEqual(Problem.CV_SUCCESS, flag);
     232
     233
     234      // var q = Problem.N_VNew_Serial(1);
     235      // Problem.NV_Set_Ith_S(q, 0, 0.0);
     236      // flag = Problem.CVodeQuadInit(cvode_mem, FQ, q);
     237      // Assert.AreEqual(Problem.CV_SUCCESS, flag);
     238
     239
     240      // int steps = 10;
     241      // flag = Problem.CVodeAdjInit(cvode_mem, steps, Problem.CV_HERMITE);
     242      // Assert.AreEqual(Problem.CV_SUCCESS, flag);
     243
     244      // step by step forward integration
     245      double t = 0.0;
     246      double tout = 0.1; // first output time
     247      int nout = 1000; // number of output times
     248      for (int iout = 0; iout < nout; iout++) {
     249        flag = Problem.CVode(cvode_mem, tout, y, ref t, Problem.CV_NORMAL);
     250        Assert.AreEqual(Problem.CV_SUCCESS, flag);
     251
     252        // get sensitivities
     253        flag = Problem.CVodeGetSens(cvode_mem, ref t, yS0);
     254        Assert.AreEqual(Problem.CV_SUCCESS, flag);
     255        unsafe {
     256          var ySP0 = *(IntPtr*)yS0.ToPointer(); // parameter of sensitivities for first component
     257          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));
     258        }
     259
     260        tout += 0.1;
     261      }
     262
     263
     264
     265      Problem.N_VDestroy_Serial(y);
     266      Problem.CVodeFree(cvode_mem);
     267      Problem.SUNLinSolFree(linearSolver);
     268      Problem.SUNMatDestroy(A);
     269    }
     270
     271
    189272
    190273    public static int F(
     
    194277      IntPtr user_data) {
    195278
     279      // the system
     280      // ∂y1/∂t = y2
     281      // ∂y2/∂t = -0.3 y1
    196282      Problem.NV_Set_Ith_S(ydot, 0, Problem.NV_Get_Ith_S(y, 1));
    197283      Problem.NV_Set_Ith_S(ydot, 1, -0.3 * Problem.NV_Get_Ith_S(y, 0));
     
    199285      ;
    200286    }
     287
     288    // computes ∂f/∂y or an approximation
     289    private int JacF(double t, IntPtr y, IntPtr fy, IntPtr Jac, IntPtr user_data, IntPtr tmp1, IntPtr tmp2, IntPtr tmp3) {
     290      // t is current time
     291      // y is current y
     292      // fy is current f(t,y)
     293      // Jac is output
     294
     295      // Matrix:
     296      //  ∂f1/∂y1  ∂f1/∂y2
     297      //  ∂f2/∂y1  ∂f2/∂y2
     298
     299      Problem.SUNDenseMatrix_Set(Jac, 0, 0, 0.0);
     300      Problem.SUNDenseMatrix_Set(Jac, 0, 1, 1.0);
     301      Problem.SUNDenseMatrix_Set(Jac, 1, 0, -0.3);
     302      Problem.SUNDenseMatrix_Set(Jac, 1, 1, 0.0);
     303      return 0;
     304    }
     305
     306    private int JacFB(double t, IntPtr y, IntPtr yB, IntPtr fyB, IntPtr Jac, IntPtr user_data, IntPtr tmp1, IntPtr tmp2, IntPtr tmp3) {
     307      throw new NotImplementedException();
     308      return 0;
     309    }
     310
     311    // sensitivity function must compute the vector (∂f/∂y)s_i(t) + ∂f/∂p_i and store in ySdot[i]
     312    public static int FS(
     313      int Ns,
     314      double t,
     315      IntPtr y, // N_Vector
     316      IntPtr ydot, // N_Vector
     317      IntPtr yS, // N_Vector* one vector for each parameter
     318      IntPtr ySdot, // N_Vector* one vector for each parameter
     319      IntPtr user_data,
     320      IntPtr tmp1, // N_Vector
     321      IntPtr tmp2 // N_Vector
     322      ) {
     323
     324      // (∂f /∂y)s_i(t) + ∂f /∂p_i
     325      var y1 = Problem.NV_Get_Ith_S(y, 0); var y2 = Problem.NV_Get_Ith_S(y, 1);
     326
     327      unsafe {
     328        var yS_p0 = *(IntPtr*)yS.ToPointer();
     329        var s1 = Problem.NV_Get_Ith_S(yS_p0, 0); var s2 = Problem.NV_Get_Ith_S(yS_p0, 1);
     330
     331
     332        // (∂f /∂y)s_i(t)  = Jac s_i(t)
     333        var sd1 = s2;
     334        var sd2 = -0.3 * s1;
     335
     336        var s_p0 = *(IntPtr*)ySdot.ToPointer();
     337        Problem.NV_Set_Ith_S(s_p0, 0, sd1 + 0.0);
     338        Problem.NV_Set_Ith_S(s_p0, 1, sd2 - 0.3);
     339      }
     340
     341      return 0;
     342    }
     343
    201344    public static int FB(
    202345        double t, // realtype
  • branches/2925_AutoDiffForDynamicalModels/HeuristicLab.Problems.DynamicalSystemsModelling/3.3/Problem.cs

    r16245 r16246  
    223223    public const int CV_REIFWD_FAIL = -105;
    224224    public const int CV_FWD_FAIL = -106;
     225
     226
    225227    public const int CV_GETY_BADT = -107;
    226228
     
    284286      );
    285287
     288
     289    // to calculate sensitivities RHS for all equations at once
     290    // must compute (∂f/∂y)s_i(t) + ∂f/∂p_i
     291    [UnmanagedFunctionPointer(CallingConvention.Cdecl)]
     292    public delegate int CVSensRhsFn(
     293      int Ns,
     294      double t,
     295      IntPtr y, // N_Vector
     296      IntPtr ydot, // N_Vector
     297      IntPtr yS, // N_Vector* one vector for each parameter
     298      IntPtr ySdot, // N_Vector* one vector for each parameter
     299      IntPtr user_data,
     300      IntPtr tmp1, // N_Vector
     301      IntPtr tmp2 // N_Vector
     302    );
     303
    286304    [DllImport("sundials_cvodes.dll", EntryPoint = "CVodeCreate", ExactSpelling = true, CallingConvention = CallingConvention.Cdecl)]
    287305    // returns a void* to the cvodes memory block if successful otherwise NULL
     
    289307
    290308    [DllImport("sundials_cvodes.dll", EntryPoint = "CVodeInit", ExactSpelling = true, CallingConvention = CallingConvention.Cdecl)]
    291 
    292309    public static extern int CVodeInit(
    293310      IntPtr cvode_mem, // pointer returned by CVodeCreate
     
    316333      IntPtr cvode_mem,
    317334      CVDlsJacFunc jacFunc
     335      );
     336
     337    [DllImport("sundials_cvodes.dll", EntryPoint = "CVodeSensInit", ExactSpelling = true, CallingConvention = CallingConvention.Cdecl)]
     338    public static extern int CVodeSensInit(
     339      IntPtr cvode_mem,
     340      int Ns, // number of parameters
     341      int ism, // sensitivity solution method, CV_SIMULTANEOUS or CV_STAGGERED
     342      CVSensRhsFn fS, // right hand side function which computes all sensitivity RHS at the same time
     343      IntPtr yS0 // N_Vector
     344      );
     345
     346    [DllImport("sundials_cvodes.dll", EntryPoint = "CVodeSensSStolerances", ExactSpelling = true, CallingConvention = CallingConvention.Cdecl)]
     347    public static extern int CVodeSensSStolerances(
     348      IntPtr cvode_mem,
     349      double reltol,
     350      double[] abstol
     351    );
     352
     353    /* Call CVodeSensEEtolerances to estimate tolerances for sensitivity
     354   variables based on the rolerances supplied for states variables and
     355   the scaling factor pbar */
     356    [DllImport("sundials_cvodes.dll", EntryPoint = "CVodeSensEEtolerances", ExactSpelling = true, CallingConvention = CallingConvention.Cdecl)]
     357    public static extern int CVodeSensEEtolerances(
     358      IntPtr cvode_mem
     359      );
     360
     361    [DllImport("sundials_cvodes.dll", EntryPoint = "CVodeGetSens", ExactSpelling = true, CallingConvention = CallingConvention.Cdecl)]
     362    public static extern int CVodeGetSens(
     363      IntPtr cvode_mem,
     364      ref double tret,
     365      IntPtr yS //N_Vector*, one vector for each parameter
    318366      );
    319367
     
    452500    [DllImport("sundials_cvodes.dll", EntryPoint = "SUNMatDestroy", ExactSpelling = true, CallingConvention = CallingConvention.Cdecl)]
    453501    public static extern void SUNMatDestroy(IntPtr A);
     502
     503    [DllImport("sundials_cvodes.dll", EntryPoint = "SUNDenseMatrix_Data", ExactSpelling = true, CallingConvention = CallingConvention.Cdecl)]
     504    public unsafe static extern double* SUNDenseMatrix_Data(IntPtr matrix);
     505
     506    [DllImport("sundials_cvodes.dll", EntryPoint = "SUNDenseMatrix_Cols", ExactSpelling = true, CallingConvention = CallingConvention.Cdecl)]
     507    public static extern long SUNDenseMatrix_Cols(IntPtr matrix);
     508
     509    [DllImport("sundials_cvodes.dll", EntryPoint = "SUNDenseMatrix_Rows", ExactSpelling = true, CallingConvention = CallingConvention.Cdecl)]
     510    public static extern long SUNDenseMatrix_Rows(IntPtr matrix);
     511
     512
     513    public unsafe static double SUNDenseMatrix_Get(IntPtr mat, long i, long j) {
     514      long M = SUNDenseMatrix_Rows(mat);
     515      // the(i, j)th element of A(with 0 <= i<M and 0 <= j<N) is given by(A->data)[j*M+i]
     516      return SUNDenseMatrix_Data(mat)[j*M+i];
     517    }
     518
     519    public unsafe static void SUNDenseMatrix_Set(IntPtr mat, long i, long j, double val) {
     520      long M = SUNDenseMatrix_Rows(mat);
     521      // the(i, j)th element of A(with 0 <= i<M and 0 <= j<N) is given by(A->data)[j*M+i]
     522      SUNDenseMatrix_Data(mat)[j * M + i] = val;
     523    }
     524
    454525    #endregion
    455526
     
    515586
    516587    [DllImport("sundials_cvodes.dll", EntryPoint = "N_VCloneVectorArray_Serial", ExactSpelling = true, CallingConvention = CallingConvention.Cdecl)]
    517     public static extern IntPtr N_VCloneVectorArray_Serial(int count, IntPtr vec);
     588    public static extern IntPtr N_VCloneVectorArray_Serial(int count, IntPtr vec); // returns N_Vector* !
     589
    518590
    519591    /*
Note: See TracChangeset for help on using the changeset viewer.