Free cookie consent management tool by TermsFeed Policy Generator

Changeset 16245


Ignore:
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:
2 edited

Legend:

Unmodified
Added
Removed
  • branches/2925_AutoDiffForDynamicalModels/AutoDiffForDynamicalModelsTest/TestCvodes.cs

    r16227 r16245  
    2121        int own_data = *(content + 2);
    2222        Console.WriteLine(own_data);
    23         double* data = (double*) *(content + 3);
     23        double* data = (double*)*(content + 3);
    2424        double data0 = *data;
    2525      }
     
    4040      var cvode_mem = Problem.CVodeCreate(Problem.MultistepMethod.CV_ADAMS, Problem.NonlinearSolverIteration.CV_FUNCTIONAL);
    4141
    42      var flag = Problem.CVodeInit(cvode_mem, F, 0.0, y);
     42      var flag = Problem.CVodeInit(cvode_mem, F, 0.0, y);
    4343      Assert.AreEqual(Problem.CV_SUCCESS, flag);
    4444
     
    5050
    5151
    52       var linearSolver = Problem.SUNDenseLinearSolver(y, A); 
     52      var linearSolver = Problem.SUNDenseLinearSolver(y, A);
    5353      Assert.AreNotSame(linearSolver, IntPtr.Zero);
    5454
     
    5656      Assert.AreEqual(Problem.CV_SUCCESS, flag);
    5757
    58       flag = Problem.CVDlsSetJacFn(cvode_mem, JacF);
    59       Assert.AreEqual(Problem.CV_SUCCESS, flag);
    60 
    61 
     58      // flag = Problem.CVDlsSetJacFn(cvode_mem, JacF);
     59      // Assert.AreEqual(Problem.CV_SUCCESS, flag);
     60
     61      // var ns = 1; // number of parameters
     62      // var p = new double[1]; // set as user-data
     63      // var yS0 = Problem.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 = Problem.N_VNew_Serial(1);
     70      Problem.NV_Set_Ith_S(q, 0, 0.0);
     71      flag = Problem.CVodeQuadInit(cvode_mem, FQ, q);
     72      Assert.AreEqual(Problem.CV_SUCCESS, flag);
     73
     74
     75      int steps = 10;
     76      flag = Problem.CVodeAdjInit(cvode_mem, steps, Problem.CV_HERMITE);
     77      Assert.AreEqual(Problem.CV_SUCCESS, flag);
     78
     79      /* step by step forward integration
    6280      double t = 0.0;
    6381      double tout = 0.1; // first output time
     
    7088        tout += 0.1;
    7189      }
     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 = Problem.CVodeF(cvode_mem, tout, y, ref time, Problem.CV_NORMAL, ref ncheck);
     97      Assert.AreEqual(Problem.CV_SUCCESS, flag);
     98
     99      long numSteps = 0;
     100      flag = Problem.CVodeGetNumSteps(cvode_mem, ref numSteps);
     101      Assert.AreEqual(Problem.CV_SUCCESS, flag);
     102
     103      var yB = Problem.N_VNew_Serial(numberOfEquations);
     104      Problem.N_VConst_Serial(0.0, yB);
     105
     106      int numberOfParameters = 2;
     107      var qB = Problem.N_VNew_Serial(numberOfParameters);
     108      Problem.N_VConst_Serial(0.0, qB);
     109
     110      int indexB = 0;
     111      flag = Problem.CVodeCreateB(cvode_mem, Problem.MultistepMethod.CV_BDF, Problem.NonlinearSolverIteration.CV_NEWTON, ref indexB);
     112      Assert.AreEqual(Problem.CV_SUCCESS, flag);
     113
     114      var TB1 = tout;
     115      flag = Problem.CVodeInitB(cvode_mem, indexB, FB, TB1, yB);
     116      Assert.AreEqual(Problem.CV_SUCCESS, flag);
     117
     118      var relTolB = 1E-6;
     119      var absTolB = 1E-8;
     120      flag = Problem.CVodeSStolerancesB(cvode_mem, indexB, relTolB, absTolB);
     121      Assert.AreEqual(Problem.CV_SUCCESS, flag);
     122
     123      var AB = Problem.SUNDenseMatrix(numberOfEquations, numberOfEquations);
     124      Assert.AreNotSame(linearSolver, IntPtr.Zero);
     125
     126      var lsB = Problem.SUNDenseLinearSolver(yB, AB);
     127      Assert.AreNotSame(linearSolver, IntPtr.Zero);
     128
     129      flag = Problem.CVDlsSetLinearSolverB(cvode_mem, indexB, lsB, AB);
     130      Assert.AreEqual(Problem.CV_SUCCESS, flag);
     131
     132      flag = Problem.CVDlsSetJacFnB(cvode_mem, indexB, JacFB);
     133      Assert.AreEqual(Problem.CV_SUCCESS, flag);
     134
     135      flag = Problem.CVodeQuadInitB(cvode_mem, indexB, FQB, qB);
     136      Assert.AreEqual(Problem.CV_SUCCESS, flag);
     137
     138      /* First get results at t = TBout1 */
     139
     140      /* Call CVodeB to integrate the backward ODE problem. */
     141      var tBackOut = 50.0;
     142      flag = Problem.CVodeB(cvode_mem, tBackOut, Problem.CV_NORMAL);
     143      Assert.AreEqual(Problem.CV_SUCCESS, flag);
     144
     145      /* Call CVodeGetB to get yB of the backward ODE problem. */
     146      flag = Problem.CVodeGetB(cvode_mem, indexB, ref time, yB);
     147      Assert.AreEqual(Problem.CV_SUCCESS, flag);
     148
     149      /* Call CVodeGetAdjY to get the interpolated value of the forward solution
     150         y during a backward integration. */
     151      flag = Problem.CVodeGetAdjY(cvode_mem, tBackOut, y);
     152      Assert.AreEqual(Problem.CV_SUCCESS, flag);
     153
     154      /* Then at t = T0 */
     155
     156      double t0 = 0.0;
     157      flag = Problem.CVodeB(cvode_mem, t0, Problem.CV_NORMAL);
     158      Assert.AreEqual(Problem.CV_SUCCESS, flag);
     159
     160      long nstB = 0;
     161      Problem.CVodeGetNumSteps(Problem.CVodeGetAdjCVodeBmem(cvode_mem, indexB), ref nstB);
     162
     163      flag = Problem.CVodeGetB(cvode_mem, indexB, ref time, yB);
     164      Assert.AreEqual(Problem.CV_SUCCESS, flag);
     165
     166
     167      flag = Problem.CVodeGetQuadB(cvode_mem, indexB, ref time, qB);
     168      Assert.AreEqual(Problem.CV_SUCCESS, flag);
     169
     170      flag = Problem.CVodeGetAdjY(cvode_mem, t0, y);
     171      Assert.AreEqual(Problem.CV_SUCCESS, flag);
    72172
    73173      Problem.N_VDestroy_Serial(y);
     
    77177    }
    78178
     179
    79180    private int JacF(double t, IntPtr y, IntPtr fy, IntPtr Jac, IntPtr user_data, IntPtr tmp1, IntPtr tmp2, IntPtr tmp3) {
    80       // TODO
     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();
    81187      return 0;
    82188    }
     
    93199      ;
    94200    }
     201    public static int FB(
     202        double t, // realtype
     203        IntPtr y, // N_Vector
     204        IntPtr yB, // N_Vector
     205        IntPtr yBdot, // N_Vector
     206        IntPtr user_data) {
     207
     208      // y' = f(y,p)
     209
     210      // yB = λ
     211      // λ' = -(∂f / ∂y)^T λ - (∂g / ∂y)
     212      // the goal is to find dG/dp where G = integrate( g(y,t,p), t=0, t=T) )
     213
     214      // for F above:
     215      // ∂f / ∂ =
     216      //  0.0 1.0
     217      // -0.3 0.0
     218
     219      // we have no g!?
     220      // therefore λ' =
     221      // λ1' = 0.3 λ2 - 0.0
     222      // λ2' = 1.0 λ1 - 0.0
     223
     224      Problem.NV_Set_Ith_S(yBdot, 0, 0.3 * Problem.NV_Get_Ith_S(yB, 1));
     225      Problem.NV_Set_Ith_S(yBdot, 1, Problem.NV_Get_Ith_S(yB, 0));
     226
     227      return 0;
     228      ;
     229    }
     230
     231    private int FQ(double t, IntPtr y, IntPtr yQdot, IntPtr user_data) {
     232      // TODO: squared error
     233      Problem.NV_Set_Ith_S(yQdot, 0, Problem.NV_Get_Ith_S(y, 2));
     234      return 0;
     235    }
     236
     237
     238    private int FQB(double t, IntPtr y, IntPtr yB, IntPtr qBdot, IntPtr user_data) {
     239      throw new NotImplementedException();
     240    }
    95241  }
    96242}
  • branches/2925_AutoDiffForDynamicalModels/HeuristicLab.Problems.DynamicalSystemsModelling/3.3/Problem.cs

    r16226 r16245  
    124124  }
    125125
     126
     127  // Eine weitere Möglichkeit ist spline-smoothing der Daten (über Zeit) um damit für jede Zielvariable
     128  // einen bereinigten (Rauschen) Wert und die Ableitung dy/dt für alle Beobachtungen zu bekommen
     129  // danach kann man die Regression direkt für dy/dt anwenden (mit bereinigten y als inputs)
    126130  [Item("Dynamical Systems Modelling Problem", "TODO")]
    127131  [Creatable(CreatableAttribute.Categories.GeneticProgrammingProblems, Priority = 900)]
     
    226230        IntPtr y, // N_Vector
    227231        IntPtr ydot, // N_Vector
     232        IntPtr user_data
     233      );
     234
     235    [UnmanagedFunctionPointer(CallingConvention.Cdecl)]
     236    public delegate int CVRhsFuncB(
     237        double t, // realtype
     238        IntPtr y, // N_Vector
     239        IntPtr yB, // N_Vector
     240        IntPtr yBdot, // N_Vector
    228241        IntPtr user_data
    229242      );
     
    241254      );
    242255
     256    [UnmanagedFunctionPointer(CallingConvention.Cdecl)]
     257    public delegate int CVDlsJacFuncB(
     258      double t,
     259      IntPtr y, // N_Vector
     260      IntPtr yB, // N_Vector
     261      IntPtr fyB, // N_Vector
     262      IntPtr Jac, // SUNMatrix
     263      IntPtr user_data,
     264      IntPtr tmp1, // N_Vector
     265      IntPtr tmp2, // N_Vector
     266      IntPtr tmp3 // N_Vector
     267    );
     268
     269    [UnmanagedFunctionPointer(CallingConvention.Cdecl)]
     270    public delegate int CVQuadRhsFn(
     271      double t,
     272      IntPtr y, // N_Vector
     273      IntPtr yQdot, // N_Vector
     274      IntPtr user_data
     275      );
     276
     277    [UnmanagedFunctionPointer(CallingConvention.Cdecl)]
     278    public delegate int CVQuadRhsFnB(
     279      double t,
     280      IntPtr y, // N_Vector
     281      IntPtr yB, // N_Vector
     282      IntPtr qBdot, // N_Vector
     283      IntPtr user_data
     284      );
    243285
    244286    [DllImport("sundials_cvodes.dll", EntryPoint = "CVodeCreate", ExactSpelling = true, CallingConvention = CallingConvention.Cdecl)]
     
    274316      IntPtr cvode_mem,
    275317      CVDlsJacFunc jacFunc
     318      );
     319
     320    [DllImport("sundials_cvodes.dll", EntryPoint = "CVodeQuadInit", ExactSpelling = true, CallingConvention = CallingConvention.Cdecl)]
     321    public static extern int CVodeQuadInit(
     322      IntPtr cvode_mem,
     323      CVQuadRhsFn qF,
     324      IntPtr yQ0 // N_Vector, initial value of yQ
     325      );
     326
     327    [DllImport("sundials_cvodes.dll", EntryPoint = "CVodeQuadInitB", ExactSpelling = true, CallingConvention = CallingConvention.Cdecl)]
     328    public static extern int CVodeQuadInitB(
     329      IntPtr cvode_mem,
     330      int indexB,
     331      CVQuadRhsFnB rhsQB,
     332      IntPtr yQB0 // N_Vector, initial value of yQB
     333    );
     334
     335
     336    [DllImport("sundials_cvodes.dll", EntryPoint = "CVodeAdjInit", ExactSpelling = true, CallingConvention = CallingConvention.Cdecl)]
     337    public static extern int CVodeAdjInit(
     338      IntPtr cvode_mem,
     339      int nd,  // integration steps between checkpoints
     340      int interpType // either CV_POLYNOMIAL or CV_HERMITE
     341      );
     342
     343    [DllImport("sundials_cvodes.dll", EntryPoint = "CVodeF", ExactSpelling = true, CallingConvention = CallingConvention.Cdecl)]
     344    public static extern int CVodeF(
     345        IntPtr cvode_mem,
     346        double t_out, // the next time at which a computed solution is desired
     347        IntPtr y, // N_Vector, the solution vector y
     348        ref double t_ret, // the time reached by the solver (output)
     349        int itask, // CV_NORMAL or CV_ONE_STEP
     350        ref int ncheck  // the number of internal checkpoints stored so far.
     351      );
     352
     353    [DllImport("sundials_cvodes.dll", EntryPoint = "CVodeGetNumSteps", ExactSpelling = true, CallingConvention = CallingConvention.Cdecl)]
     354    public static extern int CVodeGetNumSteps(
     355      IntPtr cvode_mem,
     356      ref long ncheck  // the number of steps taken
     357    );
     358
     359    [DllImport("sundials_cvodes.dll", EntryPoint = "CVodeCreateB", ExactSpelling = true, CallingConvention = CallingConvention.Cdecl)]
     360    public static extern int CVodeCreateB(
     361        IntPtr cvode_mem,
     362        MultistepMethod lmmB,
     363        NonlinearSolverIteration iterB,
     364        ref int which
     365      );
     366
     367    [DllImport("sundials_cvodes.dll", EntryPoint = "CVodeInitB", ExactSpelling = true, CallingConvention = CallingConvention.Cdecl)]
     368    public static extern int CVodeInitB(
     369        IntPtr cvode_mem,
     370        int which,
     371        CVRhsFuncB rightHandSide,
     372        double tB0, // endpoint T where final conditions are provided (equal to endpoint of forward integration)
     373        IntPtr yB0 // N_Vector inital value at t=tb0 of backward solution
    276374      );
    277375
     
    285383      );
    286384
     385
     386    [DllImport("sundials_cvodes.dll", EntryPoint = "CVodeSStolerancesB", ExactSpelling = true, CallingConvention = CallingConvention.Cdecl)]
     387    public static extern int CVodeSStolerancesB(
     388      IntPtr cvode_mem,
     389      int which,
     390      double reltol,
     391      double abstol
     392    );
     393
     394    [DllImport("sundials_cvodes.dll", EntryPoint = "CVDlsSetLinearSolverB", ExactSpelling = true, CallingConvention = CallingConvention.Cdecl)]
     395    public static extern int CVDlsSetLinearSolverB(
     396      IntPtr cvode_mem,
     397      int which,
     398      IntPtr linearSolver, // SUNLinearSolver
     399      IntPtr j // SUNMatrix
     400    );
     401
     402    [DllImport("sundials_cvodes.dll", EntryPoint = "CVDlsSetJacFnB", ExactSpelling = true, CallingConvention = CallingConvention.Cdecl)]
     403    public static extern int CVDlsSetJacFnB(
     404      IntPtr cvode_mem,
     405      int which,
     406      CVDlsJacFuncB jacFunc
     407    );
     408
     409    [DllImport("sundials_cvodes.dll", EntryPoint = "CVodeB", ExactSpelling = true, CallingConvention = CallingConvention.Cdecl)]
     410    public static extern int CVodeB(
     411      IntPtr cvode_mem,
     412      double tout, // next time at which a solution is desired
     413      int itask // flag indicating the job of the solver for the next step.
     414    );
     415
     416    [DllImport("sundials_cvodes.dll", EntryPoint = "CVodeGetB", ExactSpelling = true, CallingConvention = CallingConvention.Cdecl)]
     417    public static extern int CVodeGetB(
     418          IntPtr cvode_mem,
     419          int which,
     420          ref double tret,
     421          IntPtr yB // the backward solution at time tret
     422        );
     423    [DllImport("sundials_cvodes.dll", EntryPoint = "CVodeGetAdjY", ExactSpelling = true, CallingConvention = CallingConvention.Cdecl)]
     424    public static extern int CVodeGetAdjY(
     425          IntPtr cvode_mem,
     426          double t,
     427          IntPtr y // the forward solution y(t)
     428        );
     429
     430    [DllImport("sundials_cvodes.dll", EntryPoint = "CVodeGetAdjCVodeBmem", ExactSpelling = true, CallingConvention = CallingConvention.Cdecl)]
     431    public static extern IntPtr CVodeGetAdjCVodeBmem(
     432              IntPtr cvode_mem,
     433              int which
     434            );
     435
     436    [DllImport("sundials_cvodes.dll", EntryPoint = "CVodeGetQuadB", ExactSpelling = true, CallingConvention = CallingConvention.Cdecl)]
     437    public static extern int CVodeGetQuadB(
     438        IntPtr cvode_mem,
     439        int which,
     440        ref double tret,
     441        IntPtr yQB // N_Vector
     442      );
     443
    287444    [DllImport("sundials_cvodes.dll", EntryPoint = "CVodeFree", ExactSpelling = true, CallingConvention = CallingConvention.Cdecl)]
    288445
     
    357514    public unsafe static extern double* N_VGetArrayPointer_Serial(IntPtr vec);
    358515
     516    [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);
     518
    359519    /*
    360520#define NV_CONTENT_S(v)  ( (N_VectorContent_Serial)(v->content) )
     
    368528#define NV_Ith_S(v,i)    ( NV_DATA_S(v)[i] )
    369529*/
    370 // methods for macros
     530    // methods for macros
    371531    public unsafe static int* NV_CONTENT_S(IntPtr v) {
    372532      int* content = (int*)*(int*)v.ToPointer();
Note: See TracChangeset for help on using the changeset viewer.