TestAdaoHelper.cxx

   1 #include <vector>
   2 #include "PyObjectRAII.hxx"
   3
   4 #include <cmath>
   5
   6 /* func pour test3DVar testBlue et testNonLinearLeastSquares*/
   7 std::vector<double> funcBase(const std::vector<double>& vec)
   8 {
   9   return {vec[0],2.*vec[1],3.*vec[2],vec[0]+2.*vec[1]+3.*vec[2]};
  10 }
  11
  12 double funcCrueInternal(double Q, double K_s)
  13 {
  14   constexpr double L(5.0e3);
  15   constexpr double B(300.);
  16   constexpr double Z_v(49.);
  17   constexpr double Z_m(51.);
  18   constexpr double alpha( (Z_m - Z_v)/L );
  19   double H(pow((Q/(K_s*B*sqrt(alpha))),(3.0/5.0)));
  20   return H;
  21 }
  22
  23 /* func pour testCasCrue*/
  24 std::vector<double> funcCrue(const std::vector<double>& vec)
  25 {
  26   double K_s(vec[0]);
  27   constexpr double Qs[]={10.,20.,30.,40.};
  28   constexpr size_t LEN(sizeof(Qs)/sizeof(double));
  29   std::vector<double> ret(LEN);
  30   for(std::size_t i=0;i<LEN;++i)
  31     {
  32       ret[i] = funcCrueInternal(Qs[i],K_s);
  33     }
  34   return ret;
  35 }
  36
  37 /* Visitor commun pour test3DVar testBlue et testNonLinearLeastSquares*/
  38 class Visitor2 : public AdaoModel::PythonLeafVisitor
  39 {
  40 public:
  41   Visitor2(PyObject *context):_context(context)
  42   {
  43     std::vector< std::vector<double> > bounds{ {0., 10.}, {3., 13.}, {1.5, 15.5} };
  44     std::vector< double > Xb{5.,7.,9.};
  45     py2cpp::PyPtr boundsPy(py2cpp::toPyPtr(bounds));
  46     _bounds = boundsPy.get();
  47     Py_XINCREF(_bounds);
  48     py2cpp::PyPtr XbPy(py2cpp::toPyPtr(Xb));
  49     _Xb = XbPy.get();
  50     Py_XINCREF(_Xb);
  51     std::vector<double> observation{2., 6., 12., 20.};
  52     py2cpp::PyPtr observationPy(py2cpp::toPyPtr(observation));
  53     _observation = observationPy.get();
  54     Py_XINCREF(_observation);
  55   }
  56
  57   void visit(AdaoModel::MainModel *godFather, AdaoModel::PyObjKeyVal *obj) override
  58   {
  59     if(obj->getKey()=="Bounds")
  60       {
  61         std::ostringstream oss; oss << "___" << _cnt++;
  62         std::string varname(oss.str());
  63         obj->setVal(_bounds);
  64         PyDict_SetItemString(_context,varname.c_str(),_bounds);
  65         obj->setVarName(varname);
  66         return ;
  67       }
  68     if(godFather->findPathOf(obj)=="Background/Vector")
  69       {
  70         std::ostringstream oss; oss << "___" << _cnt++;
  71         std::string varname(oss.str());
  72         obj->setVal(_Xb);
  73         PyDict_SetItemString(_context,varname.c_str(),_Xb);
  74         obj->setVarName(varname);
  75       }
  76     if(godFather->findPathOf(obj)=="Observation/Vector")
  77       {
  78         std::ostringstream oss; oss << "____" << _cnt++;
  79         std::string varname(oss.str());
  80         obj->setVal(_observation);
  81         PyDict_SetItemString(_context,varname.c_str(),_observation);
  82         obj->setVarName(varname);
  83       }
  84   }
  85 private:
  86   unsigned int _cnt = 0;
  87   PyObject *_bounds = nullptr;
  88   PyObject *_Xb = nullptr;
  89   PyObject *_observation = nullptr;
  90   PyObject *_context = nullptr;
  91 };
  92
  93 /* Visitor pour testCasCrue */
  94 class VisitorCruePython : public AdaoModel::PythonLeafVisitor
  95 {
  96 public:
  97   VisitorCruePython(PyObject *context):_context(context)
  98   {
  99     {//case.set( 'Background',          Vector=thetaB)
 100       std::vector< double > Xb{ 20. };//thetaB
 101       py2cpp::PyPtr XbPy(py2cpp::toPyPtr(Xb));
 102       _Xb = XbPy.get();
 103       Py_XINCREF(_Xb);
 104     }
 105     {//case.set( 'BackgroundError',     DiagonalSparseMatrix=sigmaTheta )
 106       std::vector< double > sigmaTheta{ 5.e10 };
 107       py2cpp::PyPtr sigmaThetaPy(py2cpp::toPyPtr(sigmaTheta));
 108       _sigmaTheta = sigmaThetaPy.get();
 109       Py_XINCREF(_sigmaTheta);
 110     }
 111     {//case.set( 'Observation',         Vector=Hobs)
 112       std::vector<double> observation{0.19694513, 0.298513, 0.38073079, 0.45246109};
 113       py2cpp::PyPtr observationPy(py2cpp::toPyPtr(observation));
 114       _observation = observationPy.get();
 115       Py_XINCREF(_observation);
 116     }
 117     {//case.set( 'ObservationError',    ScalarSparseMatrix=sigmaH )
 118       double sigmaH( 0.5);
 119       py2cpp::PyPtr sigmaHPy(py2cpp::toPyPtr(sigmaH));
 120       _sigmaH = sigmaHPy.get();
 121       Py_XINCREF(_sigmaH);
 122     }
 123   }
 124
 125
 126   void visit(AdaoModel::MainModel *godFather, AdaoModel::PyObjKeyVal *obj) override
 127   {
 128     if(godFather->findPathOf(obj)=="Background/Vector")
 129       {
 130         std::ostringstream oss; oss << "___" << _cnt++;
 131         std::string varname(oss.str());
 132         obj->setVal(_Xb);
 133         PyDict_SetItemString(_context,varname.c_str(),_Xb);
 134         obj->setVarName(varname);
 135       }
 136     if(godFather->findPathOf(obj)=="BackgroundError/Matrix")
 137       {
 138         std::ostringstream oss; oss << "___" << _cnt++;
 139         std::string varname(oss.str());
 140         obj->setVal(_sigmaTheta);
 141         PyDict_SetItemString(_context,varname.c_str(),_Xb);
 142         obj->setVarName(varname);
 143       }
 144     if(godFather->findPathOf(obj)=="Observation/Vector")
 145       {
 146         std::ostringstream oss; oss << "____" << _cnt++;
 147         std::string varname(oss.str());
 148         obj->setVal(_observation);
 149         PyDict_SetItemString(_context,varname.c_str(),_observation);
 150         obj->setVarName(varname);
 151       }
 152     if(godFather->findPathOf(obj)=="ObservationError/Matrix")
 153       {
 154         std::ostringstream oss; oss << "____" << _cnt++;
 155         std::string varname(oss.str());
 156         obj->setVal(_sigmaH);
 157         PyDict_SetItemString(_context,varname.c_str(),_sigmaH);
 158         obj->setVarName(varname);
 159       }
 160   }
 161 private:
 162   unsigned int _cnt = 0;
 163   PyObject *_Xb = nullptr;
 164   PyObject *_sigmaH = nullptr;
 165   PyObject *_sigmaTheta = nullptr;
 166   PyObject *_observation = nullptr;
 167   PyObject *_context = nullptr;
 168 };
 169