Free cookie consent management tool by TermsFeed Policy Generator

source: branches/2925_AutoDiffForDynamicalModels/AutoDiffForDynamicalModelsTest/TestCvodes.cs

Last change on this file was 16254, checked in by gkronber, 6 years ago

#2925: made usage of x86 native library explicit.
File size: 12.8 KB
Line 
1using System;
2using HeuristicLab.Problems.DynamicalSystemsModelling;
3using Microsoft.VisualStudio.TestTools.UnitTesting;
4
5namespace AutoDiffForDynamicalModelsTest {
6  [TestClass]
7  public class TestCvodes {
8    [TestMethod]
9    public unsafe void TestCvodesASAMethod() {
10      // test vectors
11      var arr = new double[] { 3.14, 2.71 };
12      var vec = CVODES.N_VMake_Serial(2, arr);
13      CVODES.N_VDestroy_Serial(vec);
14
15      vec = CVODES.N_VNew_Serial(10);
16
17      unsafe {
18        int* content = (int*)*(int*)vec.ToPointer();
19        long length = *content;
20        Console.WriteLine(*content);
21        int own_data = *(content + 2);
22        Console.WriteLine(own_data);
23        double* data = (double*)*(content + 3);
24        double data0 = *data;
25      }
26      CVODES.N_VConst_Serial(2.0, vec);
27      CVODES.N_VPrint_Serial(vec);
28      Assert.AreEqual(20, CVODES.N_VL1Norm_Serial(vec));
29      CVODES.N_VDestroy_Serial(vec);
30
31
32      // linear oscillator
33      int numberOfEquations = 2;
34      var y = CVODES.N_VNew_Serial(numberOfEquations);
35      // y must be initialized before calling CVodeInit
36      // CVODES.N_VConst_Serial(100.0, y);
37      CVODES.NV_Set_Ith_S(y, 0, 0.5); // x
38      CVODES.NV_Set_Ith_S(y, 1, 1);  // v
39
40      var cvode_mem = CVODES.CVodeCreate(CVODES.MultistepMethod.CV_ADAMS, CVODES.NonlinearSolverIteration.CV_FUNCTIONAL);
41
42      var flag = CVODES.CVodeInit(cvode_mem, F, 0.0, y);
43      Assert.AreEqual(CVODES.CV_SUCCESS, flag);
44
45      flag = CVODES.CVodeSStolerances(cvode_mem, 1E-4, 1.0);
46      Assert.AreEqual(CVODES.CV_SUCCESS, flag);
47
48      var A = CVODES.SUNDenseMatrix(numberOfEquations, numberOfEquations);
49      Assert.AreNotSame(A, IntPtr.Zero);
50
51
52      var linearSolver = CVODES.SUNDenseLinearSolver(y, A);
53      Assert.AreNotSame(linearSolver, IntPtr.Zero);
54
55      flag = CVODES.CVDlsSetLinearSolver(cvode_mem, linearSolver, A);
56      Assert.AreEqual(CVODES.CV_SUCCESS, flag);
57
58      // flag = CVODES.CVDlsSetJacFn(cvode_mem, JacF);
59      // Assert.AreEqual(CVODES.CV_SUCCESS, flag);
60
61      // var ns = 1; // number of parameters
62      // var p = new double[1]; // set as user-data
63      // var yS0 = CVODES.N_VCloneVectorArray_Serial(ns, y); // clone the output vector for each parameter, TODO: free
64      // for(int i=0;i<ns;i++) {
65      //
66      // }
67
68
69      var q = CVODES.N_VNew_Serial(1);
70      CVODES.NV_Set_Ith_S(q, 0, 0.0);
71      flag = CVODES.CVodeQuadInit(cvode_mem, FQ, q);
72      Assert.AreEqual(CVODES.CV_SUCCESS, flag);
73
74
75      int steps = 10;
76      flag = CVODES.CVodeAdjInit(cvode_mem, steps, CVODES.CV_HERMITE);
77      Assert.AreEqual(CVODES.CV_SUCCESS, flag);
78
79      /* step by step forward integration
80      double t = 0.0;
81      double tout = 0.1; // first output time
82      int nout = 100; // number of output times
83      for (int iout = 0; iout < nout; iout++) {
84        flag = CVODES.CVode(cvode_mem, tout, y, ref t, CVODES.CV_NORMAL);
85        Assert.AreEqual(CVODES.CV_SUCCESS, flag);
86        Console.WriteLine("{0} {1} {2}", t, CVODES.NV_Get_Ith_S(y, 0), CVODES.NV_Get_Ith_S(y, 1));
87        // CVODES.N_VPrint_Serial(y);
88        tout += 0.1;
89      }
90      */
91
92      // complete forward integration
93      double tout = 100.0; // last time
94      double time = 0.0;
95      int ncheck = 0; // number of checkpoints
96      flag = CVODES.CVodeF(cvode_mem, tout, y, ref time, CVODES.CV_NORMAL, ref ncheck);
97      Assert.AreEqual(CVODES.CV_SUCCESS, flag);
98
99      long numSteps = 0;
100      flag = CVODES.CVodeGetNumSteps(cvode_mem, ref numSteps);
101      Assert.AreEqual(CVODES.CV_SUCCESS, flag);
102
103      var yB = CVODES.N_VNew_Serial(numberOfEquations);
104      CVODES.N_VConst_Serial(0.0, yB);
105
106      int numberOfParameters = 2;
107      var qB = CVODES.N_VNew_Serial(numberOfParameters);
108      CVODES.N_VConst_Serial(0.0, qB);
109
110      int indexB = 0;
111      flag = CVODES.CVodeCreateB(cvode_mem, CVODES.MultistepMethod.CV_BDF, CVODES.NonlinearSolverIteration.CV_NEWTON, ref indexB);
112      Assert.AreEqual(CVODES.CV_SUCCESS, flag);
113
114      var TB1 = tout;
115      flag = CVODES.CVodeInitB(cvode_mem, indexB, FB, TB1, yB);
116      Assert.AreEqual(CVODES.CV_SUCCESS, flag);
117
118      var relTolB = 1E-6;
119      var absTolB = 1E-8;
120      flag = CVODES.CVodeSStolerancesB(cvode_mem, indexB, relTolB, absTolB);
121      Assert.AreEqual(CVODES.CV_SUCCESS, flag);
122
123      var AB = CVODES.SUNDenseMatrix(numberOfEquations, numberOfEquations);
124      Assert.AreNotSame(linearSolver, IntPtr.Zero);
125
126      var lsB = CVODES.SUNDenseLinearSolver(yB, AB);
127      Assert.AreNotSame(linearSolver, IntPtr.Zero);
128
129      flag = CVODES.CVDlsSetLinearSolverB(cvode_mem, indexB, lsB, AB);
130      Assert.AreEqual(CVODES.CV_SUCCESS, flag);
131
132      flag = CVODES.CVDlsSetJacFnB(cvode_mem, indexB, JacFB);
133      Assert.AreEqual(CVODES.CV_SUCCESS, flag);
134
135      flag = CVODES.CVodeQuadInitB(cvode_mem, indexB, FQB, qB);
136      Assert.AreEqual(CVODES.CV_SUCCESS, flag);
137
138      /* First get results at t = TBout1 */
139
140      /* Call CVodeB to integrate the backward ODE CVODES. */
141      var tBackOut = 50.0;
142      flag = CVODES.CVodeB(cvode_mem, tBackOut, CVODES.CV_NORMAL);
143      Assert.AreEqual(CVODES.CV_SUCCESS, flag);
144
145      /* Call CVodeGetB to get yB of the backward ODE CVODES. */
146      flag = CVODES.CVodeGetB(cvode_mem, indexB, ref time, yB);
147      Assert.AreEqual(CVODES.CV_SUCCESS, flag);
148
149      /* Call CVodeGetAdjY to get the interpolated value of the forward solution
150         y during a backward integration. */
151      flag = CVODES.CVodeGetAdjY(cvode_mem, tBackOut, y);
152      Assert.AreEqual(CVODES.CV_SUCCESS, flag);
153
154      /* Then at t = T0 */
155
156      double t0 = 0.0;
157      flag = CVODES.CVodeB(cvode_mem, t0, CVODES.CV_NORMAL);
158      Assert.AreEqual(CVODES.CV_SUCCESS, flag);
159
160      long nstB = 0;
161      CVODES.CVodeGetNumSteps(CVODES.CVodeGetAdjCVodeBmem(cvode_mem, indexB), ref nstB);
162
163      flag = CVODES.CVodeGetB(cvode_mem, indexB, ref time, yB);
164      Assert.AreEqual(CVODES.CV_SUCCESS, flag);
165
166
167      flag = CVODES.CVodeGetQuadB(cvode_mem, indexB, ref time, qB);
168      Assert.AreEqual(CVODES.CV_SUCCESS, flag);
169
170      flag = CVODES.CVodeGetAdjY(cvode_mem, t0, y);
171      Assert.AreEqual(CVODES.CV_SUCCESS, flag);
172
173      CVODES.N_VDestroy_Serial(y);
174      CVODES.CVodeFree(ref cvode_mem);
175      CVODES.SUNLinSolFree(linearSolver);
176      CVODES.SUNMatDestroy(A);
177    }
178
179    [TestMethod]
180    public void TestCvodesFSAMethod() {
181
182      // linear oscillator
183      int numberOfEquations = 2;
184      var y = CVODES.N_VNew_Serial(numberOfEquations);
185      // y must be initialized before calling CVodeInit
186      // CVODES.N_VConst_Serial(100.0, y);
187      CVODES.NV_Set_Ith_S(y, 0, 0.5); // x
188      CVODES.NV_Set_Ith_S(y, 1, 1);  // v
189
190      var cvode_mem = CVODES.CVodeCreate(CVODES.MultistepMethod.CV_BDF, CVODES.NonlinearSolverIteration.CV_NEWTON);
191
192      var flag = CVODES.CVodeInit(cvode_mem, F, 0.0, y);
193      Assert.AreEqual(CVODES.CV_SUCCESS, flag);
194
195      flag = CVODES.CVodeSStolerances(cvode_mem, 1E-4, 1.0);
196      Assert.AreEqual(CVODES.CV_SUCCESS, flag);
197
198      var A = CVODES.SUNDenseMatrix(numberOfEquations, numberOfEquations);
199      Assert.AreNotSame(A, IntPtr.Zero);
200      Console.WriteLine(CVODES.SUNDenseMatrix_Get(A, 0, 0));
201
202
203      var linearSolver = CVODES.SUNDenseLinearSolver(y, A);
204      Assert.AreNotSame(linearSolver, IntPtr.Zero);
205
206      flag = CVODES.CVDlsSetLinearSolver(cvode_mem, linearSolver, A);
207      Assert.AreEqual(CVODES.CV_SUCCESS, flag);
208
209      flag = CVODES.CVDlsSetJacFn(cvode_mem, JacF);
210      Assert.AreEqual(CVODES.CV_SUCCESS, flag);
211
212      var ns = 1; // number of parameters
213      var p = new double[1]; // set as user-data
214      var yS0 = CVODES.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          CVODES.N_VConst_Serial(0.0, yS0_i);
219        }
220      }
221
222
223      flag = CVODES.CVodeSensInit(cvode_mem, ns, CVODES.CV_SIMULTANEOUS, FS, yS0);
224      Assert.AreEqual(CVODES.CV_SUCCESS, flag);
225
226      // flag = CVODES.CVodeSensSStolerances(cvode_mem, reltolS, abstolS);
227      // Assert.AreEqual(CVODES.CV_SUCCESS, flag);
228      flag = CVODES.CVodeSensEEtolerances(cvode_mem);
229      Assert.AreEqual(CVODES.CV_SUCCESS, flag);
230
231
232      // var q = CVODES.N_VNew_Serial(1);
233      // CVODES.NV_Set_Ith_S(q, 0, 0.0);
234      // flag = CVODES.CVodeQuadInit(cvode_mem, FQ, q);
235      // Assert.AreEqual(CVODES.CV_SUCCESS, flag);
236
237
238      // int steps = 10;
239      // flag = CVODES.CVodeAdjInit(cvode_mem, steps, CVODES.CV_HERMITE);
240      // Assert.AreEqual(CVODES.CV_SUCCESS, flag);
241
242      // step by step forward integration
243      double t = 0.0;
244      double tout = 0.1; // first output time
245      int nout = 1000; // number of output times
246      for (int iout = 0; iout < nout; iout++) {
247        flag = CVODES.CVode(cvode_mem, tout, y, ref t, CVODES.CV_NORMAL);
248        Assert.AreEqual(CVODES.CV_SUCCESS, flag);
249
250        // get sensitivities
251        flag = CVODES.CVodeGetSens(cvode_mem, ref t, yS0);
252        Assert.AreEqual(CVODES.CV_SUCCESS, flag);
253        unsafe {
254          var ySP0 = *(IntPtr*)yS0.ToPointer(); // parameter of sensitivities for first component
255          Console.WriteLine("{0} {1} {2} {3} {4}", t, CVODES.NV_Get_Ith_S(y, 0), CVODES.NV_Get_Ith_S(y, 1), CVODES.NV_Get_Ith_S(ySP0, 0), CVODES.NV_Get_Ith_S(ySP0, 1));
256        }
257
258        tout += 0.1;
259      }
260
261
262
263      CVODES.N_VDestroy_Serial(y);
264      CVODES.CVodeFree(ref cvode_mem);
265      CVODES.SUNLinSolFree(linearSolver);
266      CVODES.SUNMatDestroy(A);
267    }
268
269
270
271    public static int F(
272      double t, // realtype
273      IntPtr y, // N_Vector
274      IntPtr ydot, // N_Vector
275      IntPtr user_data) {
276
277      // the system
278      // ∂y1/∂t = y2
279      // ∂y2/∂t = -0.3 y1
280      CVODES.NV_Set_Ith_S(ydot, 0, CVODES.NV_Get_Ith_S(y, 1));
281      CVODES.NV_Set_Ith_S(ydot, 1, -0.3 * CVODES.NV_Get_Ith_S(y, 0));
282      return 0;
283      ;
284    }
285
286    // computes ∂f/∂y or an approximation
287    private int JacF(double t, IntPtr y, IntPtr fy, IntPtr Jac, IntPtr user_data, IntPtr tmp1, IntPtr tmp2, IntPtr tmp3) {
288      // t is current time
289      // y is current y
290      // fy is current f(t,y)
291      // Jac is output
292
293      // Matrix:
294      //  ∂f1/∂y1  ∂f1/∂y2
295      //  ∂f2/∂y1  ∂f2/∂y2
296
297      CVODES.SUNDenseMatrix_Set(Jac, 0, 0, 0.0);
298      CVODES.SUNDenseMatrix_Set(Jac, 0, 1, 1.0);
299      CVODES.SUNDenseMatrix_Set(Jac, 1, 0, -0.3);
300      CVODES.SUNDenseMatrix_Set(Jac, 1, 1, 0.0);
301      return 0;
302    }
303
304    private int JacFB(double t, IntPtr y, IntPtr yB, IntPtr fyB, IntPtr Jac, IntPtr user_data, IntPtr tmp1, IntPtr tmp2, IntPtr tmp3) {
305      throw new NotImplementedException();
306      return 0;
307    }
308
309    // sensitivity function must compute the vector (∂f/∂y)s_i(t) + ∂f/∂p_i and store in ySdot[i]
310    public static int FS(
311      int Ns,
312      double t,
313      IntPtr y, // N_Vector
314      IntPtr ydot, // N_Vector
315      IntPtr yS, // N_Vector* one vector for each parameter
316      IntPtr ySdot, // N_Vector* one vector for each parameter
317      IntPtr user_data,
318      IntPtr tmp1, // N_Vector
319      IntPtr tmp2 // N_Vector
320      ) {
321
322      // (∂f /∂y)s_i(t) + ∂f /∂p_i
323      var y1 = CVODES.NV_Get_Ith_S(y, 0); var y2 = CVODES.NV_Get_Ith_S(y, 1);
324
325      unsafe {
326        var yS_p0 = *(IntPtr*)yS.ToPointer();
327        var s1 = CVODES.NV_Get_Ith_S(yS_p0, 0); var s2 = CVODES.NV_Get_Ith_S(yS_p0, 1);
328
329
330        // (∂f /∂y)s_i(t)  == Jac s_i(t)
331        var sd1 = s2;
332        var sd2 = -0.3 * s1;
333
334        var s_p0 = *(IntPtr*)ySdot.ToPointer();
335        CVODES.NV_Set_Ith_S(s_p0, 0, sd1 + 0.0);
336        CVODES.NV_Set_Ith_S(s_p0, 1, sd2 - 0.3);
337      }
338
339      return 0;
340    }
341
342    public static int FB(
343        double t, // realtype
344        IntPtr y, // N_Vector
345        IntPtr yB, // N_Vector
346        IntPtr yBdot, // N_Vector
347        IntPtr user_data) {
348
349      // y' = f(y,p)
350
351      // yB = λ
352      // λ' = -(∂f / ∂y)^T λ - (∂g / ∂y)
353      // the goal is to find dG/dp where G = integrate( g(y,t,p), t=0, t=T) )
354
355      // for F above:
356      // ∂f / ∂ =
357      //  0.0 1.0
358      // -0.3 0.0
359
360      // we have no g!?
361      // therefore λ' =
362      // λ1' = 0.3 λ2 - 0.0
363      // λ2' = 1.0 λ1 - 0.0
364
365      CVODES.NV_Set_Ith_S(yBdot, 0, 0.3 * CVODES.NV_Get_Ith_S(yB, 1));
366      CVODES.NV_Set_Ith_S(yBdot, 1, CVODES.NV_Get_Ith_S(yB, 0));
367
368      return 0;
369      ;
370    }
371
372    private int FQ(double t, IntPtr y, IntPtr yQdot, IntPtr user_data) {
373      // TODO: squared error
374      CVODES.NV_Set_Ith_S(yQdot, 0, CVODES.NV_Get_Ith_S(y, 2));
375      return 0;
376    }
377
378
379    private int FQB(double t, IntPtr y, IntPtr yB, IntPtr qBdot, IntPtr user_data) {
380      throw new NotImplementedException();
381    }
382  }
383}
Note: See TracBrowser for help on using the repository browser.