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

using System;
using HeuristicLab.Problems.DynamicalSystemsModelling;
using Microsoft.VisualStudio.TestTools.UnitTesting;

namespace AutoDiffForDynamicalModelsTest {
  [TestClass]
  public class TestCvodes {
    [TestMethod]
    public unsafe void TestCvodesMethod() {
      // test vectors
      var arr = new double[] { 3.14, 2.71 };
      var vec = Problem.N_VMake_Serial(2, arr);
      Problem.N_VDestroy_Serial(vec);

      vec = Problem.N_VNew_Serial(10);

      unsafe {
        int* content = (int*)*(int*)vec.ToPointer();
        long length = *content;
        Console.WriteLine(*content);
        int own_data = *(content + 2);
        Console.WriteLine(own_data);
        double* data = (double*)*(content + 3);
        double data0 = *data;
      }
      Problem.N_VConst_Serial(2.0, vec);
      Problem.N_VPrint_Serial(vec);
      Assert.AreEqual(20, Problem.N_VL1Norm_Serial(vec));
      Problem.N_VDestroy_Serial(vec);


      // linear oscillator
      int numberOfEquations = 2;
      var y = Problem.N_VNew_Serial(numberOfEquations);
      // y must be initialized before calling CVodeInit
      // Problem.N_VConst_Serial(100.0, y);
      Problem.NV_Set_Ith_S(y, 0, 0.5); // x
      Problem.NV_Set_Ith_S(y, 1, 1);  // v

      var cvode_mem = Problem.CVodeCreate(Problem.MultistepMethod.CV_ADAMS, Problem.NonlinearSolverIteration.CV_FUNCTIONAL);

      var flag = Problem.CVodeInit(cvode_mem, F, 0.0, y);
      Assert.AreEqual(Problem.CV_SUCCESS, flag);

      flag = Problem.CVodeSStolerances(cvode_mem, 1E-4, 1.0);
      Assert.AreEqual(Problem.CV_SUCCESS, flag);

      var A = Problem.SUNDenseMatrix(numberOfEquations, numberOfEquations);
      Assert.AreNotSame(A, IntPtr.Zero);


      var linearSolver = Problem.SUNDenseLinearSolver(y, A);
      Assert.AreNotSame(linearSolver, IntPtr.Zero);

      flag = Problem.CVDlsSetLinearSolver(cvode_mem, linearSolver, A);
      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
      int nout = 100; // number of output times
      for (int iout = 0; iout < nout; iout++) {
        flag = Problem.CVode(cvode_mem, tout, y, ref t, Problem.CV_NORMAL);
        Assert.AreEqual(Problem.CV_SUCCESS, flag);
        Console.WriteLine("{0} {1} {2}", t, Problem.NV_Get_Ith_S(y, 0), Problem.NV_Get_Ith_S(y, 1));
        // Problem.N_VPrint_Serial(y);
        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);
      Problem.CVodeFree(cvode_mem);
      Problem.SUNLinSolFree(linearSolver);
      Problem.SUNMatDestroy(A);
    }


    private int JacF(double t, IntPtr y, IntPtr fy, IntPtr Jac, IntPtr user_data, IntPtr tmp1, IntPtr tmp2, IntPtr tmp3) {
      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;
    }

    public static int F(
      double t, // realtype
      IntPtr y, // N_Vector
      IntPtr ydot, // N_Vector
      IntPtr user_data) {

      Problem.NV_Set_Ith_S(ydot, 0, Problem.NV_Get_Ith_S(y, 1));
      Problem.NV_Set_Ith_S(ydot, 1, -0.3 * Problem.NV_Get_Ith_S(y, 0));
      return 0;
      ;
    }
    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();
    }
  }
}