From 3a322a64957a19c0ed44ebf86545dea5fff5230f Mon Sep 17 00:00:00 2001 From: Anthony Geay Date: Wed, 27 Feb 2019 09:34:35 +0100 Subject: [PATCH] First implementation --- AdaoExchangeLayer4Quintet.cxx | 423 ++++++++++++++++++++ AdaoExchangeLayer4Quintet.hxx | 48 +++ AdaoModelKeyVal.cxx | 439 ++++++++++++++++++++ AdaoModelKeyVal.hxx | 466 ++++++++++++++++++++++ CMakeLists.txt | 76 ++++ Cases/Definition_complete_de_cas_3DVAR.py | 169 ++++++++ Cases/Definition_complete_de_cas_Blue.py | 166 ++++++++ Cases/Definition_complete_de_cas_NLLS.py | 164 ++++++++ PyObjectRAII.hxx | 66 +++ TestAdaoExchange.cxx | 293 ++++++++++++++ TestAdaoExchange.hxx | 40 ++ ThreeDVarCase.py | 96 +++++ start.tar.gz | Bin 0 -> 13713 bytes 13 files changed, 2446 insertions(+) create mode 100644 AdaoExchangeLayer4Quintet.cxx create mode 100644 AdaoExchangeLayer4Quintet.hxx create mode 100644 AdaoModelKeyVal.cxx create mode 100644 AdaoModelKeyVal.hxx create mode 100644 CMakeLists.txt create mode 100644 Cases/Definition_complete_de_cas_3DVAR.py create mode 100644 Cases/Definition_complete_de_cas_Blue.py create mode 100644 Cases/Definition_complete_de_cas_NLLS.py create mode 100644 PyObjectRAII.hxx create mode 100644 TestAdaoExchange.cxx create mode 100644 TestAdaoExchange.hxx create mode 100644 ThreeDVarCase.py create mode 100644 start.tar.gz diff --git a/AdaoExchangeLayer4Quintet.cxx b/AdaoExchangeLayer4Quintet.cxx new file mode 100644 index 0000000..0825d1b --- /dev/null +++ b/AdaoExchangeLayer4Quintet.cxx @@ -0,0 +1,423 @@ +// Copyright (C) 2019 EDF R&D +// +// This library is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 2.1 of the License. +// +// This library is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// Lesser General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License along with this library; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +// +// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com +// +// Author: Anthony Geay, anthony.geay@edf.fr, EDF R&D + +#include "AdaoExchangeLayer4Quintet.hxx" +#include "AdaoExchangeLayerException.hxx" +#include "AdaoModelKeyVal.hxx" +#include "PyObjectRAII.hxx" +#include "Python.h" + +#include "py2cpp/py2cpp.hxx" + +#include + +#include +#include +#include +#include +#include +#include + +struct DataExchangedBetweenThreads // data written by subthread and read by calling thread +{ +public: + DataExchangedBetweenThreads(); + ~DataExchangedBetweenThreads(); +public: + sem_t _sem; + sem_t _sem_result_is_here; + volatile bool _finished = false; + volatile PyObject *_data = nullptr; +}; + +///////////////////////////////////////////// + +struct AdaoCallbackSt +{ + PyObject_HEAD + DataExchangedBetweenThreads *_data; +}; + +static PyObject *adaocallback_call(AdaoCallbackSt *self, PyObject *args, PyObject *kw) +{ + if(!PyTuple_Check(args)) + throw AdaoExchangeLayerException("Input args is not a tuple as expected !"); + if(PyTuple_Size(args)!=1) + throw AdaoExchangeLayerException("Input args is not a tuple of size 1 as expected !"); + PyObjectRAII zeobj(PyObjectRAII::FromBorrowed(PyTuple_GetItem(args,0))); + if(zeobj.isNull()) + throw AdaoExchangeLayerException("Retrieve of elt #0 of input tuple has failed !"); + volatile PyObject *ret(nullptr); + PyThreadState *tstate(PyEval_SaveThread());// GIL is acquired (see ExecuteAsync). Before entering into non python section. Release lock + { + self->_data->_finished = false; + self->_data->_data = zeobj; + sem_post(&self->_data->_sem); + sem_wait(&self->_data->_sem_result_is_here); + ret = self->_data->_data; + } + PyEval_RestoreThread(tstate);//End of parallel section. Reaquire the GIL and restore the thread state + return (PyObject *)ret; +} + +static int adaocallback___init__(PyObject *self, PyObject *args, PyObject *kwargs) { return 0; } + +static PyObject *adaocallback___new__(PyTypeObject *type, PyObject *args, PyObject *kwargs) +{ + return (PyObject *)( type->tp_alloc(type, 0) ); +} + +static void adaocallback_dealloc(PyObject *self) +{ + Py_TYPE(self)->tp_free(self); +} + +PyTypeObject AdaoCallbackType = { + PyVarObject_HEAD_INIT(&PyType_Type, 0) + "adaocallbacktype", + sizeof(AdaoCallbackSt), + 0, + adaocallback_dealloc, /*tp_dealloc*/ + 0, /*tp_print*/ + 0, /*tp_getattr*/ + 0, /*tp_setattr*/ + 0, /*tp_compare*/ + 0, /*tp_repr*/ + 0, /*tp_as_number*/ + 0, /*tp_as_sequence*/ + 0, /*tp_as_mapping*/ + 0, /*tp_hash*/ + (ternaryfunc)adaocallback_call, /*tp_call*/ + 0, /*tp_str*/ + 0, /*tp_getattro*/ + 0, /*tp_setattro*/ + 0, /*tp_as_buffer*/ + Py_TPFLAGS_DEFAULT | Py_TPFLAGS_HAVE_GC | Py_TPFLAGS_BASETYPE, /*tp_flags*/ + 0, /*tp_doc*/ + 0, /*tp_traverse*/ + 0, /*tp_clear*/ + 0, /*tp_richcompare*/ + 0, /*tp_weaklistoffset*/ + 0, /*tp_iter*/ + 0, /*tp_iternext*/ + 0, /*tp_methods*/ + 0, /*tp_members*/ + 0, /*tp_getset*/ + 0, /*tp_base*/ + 0, /*tp_dict*/ + 0, /*tp_descr_get*/ + 0, /*tp_descr_set*/ + 0, /*tp_dictoffset*/ + adaocallback___init__, /*tp_init*/ + PyType_GenericAlloc, /*tp_alloc*/ + adaocallback___new__, /*tp_new*/ + PyObject_GC_Del, /*tp_free*/ +}; + +///////////////////////////////////////////// + +DataExchangedBetweenThreads::DataExchangedBetweenThreads() +{ + if(sem_init(&_sem,0,0)!=0)// put value to 0 to lock by default + throw AdaoExchangeLayerException("Internal constructor : Error on initialization of semaphore !"); + if(sem_init(&_sem_result_is_here,0,0)!=0)// put value to 0 to lock by default + throw AdaoExchangeLayerException("Internal constructor : Error on initialization of semaphore !"); +} + +DataExchangedBetweenThreads::~DataExchangedBetweenThreads() +{ + sem_destroy(&_sem); + sem_destroy(&_sem_result_is_here); +} + +class AdaoCallbackKeeper +{ +public: + void assign(AdaoCallbackSt *pt, DataExchangedBetweenThreads *data) + { + release(); + _pt = pt; + _pt->_data = data; + } + PyObject *getPyObject() const { return reinterpret_cast(_pt); } + ~AdaoCallbackKeeper() { release(); } +private: + void release() { if(_pt) { Py_XDECREF(_pt); } } +private: + AdaoCallbackSt *_pt = nullptr; +}; + +class AdaoExchangeLayer4Quintet::Internal +{ +public: + Internal():_context(PyObjectRAII::FromNew(PyDict_New())) + { + PyObject *mainmod(PyImport_AddModule("__main__")); + PyObject *globals(PyModule_GetDict(mainmod)); + PyObject *bltins(PyEval_GetBuiltins()); + PyDict_SetItemString(_context,"__builtins__",bltins); + } +public: + PyObjectRAII _context; + PyObjectRAII _generate_case_func; + PyObjectRAII _decorator_func; + PyObjectRAII _adao_case; + PyObjectRAII _execute_func; + AdaoCallbackKeeper _py_call_back; + std::future< void > _fut; + PyThreadState *_tstate = nullptr; + DataExchangedBetweenThreads _data_btw_threads; +}; + +wchar_t **ConvertToWChar(int argc, const char *argv[]) +{ + wchar_t **ret(new wchar_t*[argc]); + for(int i=0;i(TABW[0])); + Py_Initialize(); // Initialize the interpreter + PySys_SetArgv(1,TABW); + FreeWChar(1,TABW); + PyEval_InitThreads(); + } + delete _internal; + _internal = new Internal; +} + +class Visitor1 : public AdaoModel::PythonLeafVisitor +{ +public: + Visitor1(PyObjectRAII func, PyObject *context):_func(func),_context(context) +{ + std::vector< std::vector > bounds{ {0., 10.}, {3., 13.}, {1.5, 15.5} }; + std::vector< double > Xb{5.,7.,9.}; + py2cpp::PyPtr boundsPy(py2cpp::toPyPtr(bounds)); + _bounds = boundsPy.get(); + Py_XINCREF(_bounds); + py2cpp::PyPtr XbPy(py2cpp::toPyPtr(Xb)); + _Xb = XbPy.get(); + Py_XINCREF(_Xb); + std::vector observation{2., 6., 12., 20.}; + py2cpp::PyPtr observationPy(py2cpp::toPyPtr(observation)); + _observation = observationPy.get(); + Py_XINCREF(_observation); + +} + void visit(AdaoModel::MainModel *godFather, AdaoModel::PyObjKeyVal *obj) override + { + if(obj->getKey()=="Matrix" || obj->getKey()=="DiagonalSparseMatrix") + { + std::ostringstream oss; oss << "__" << _cnt++; + std::string varname(oss.str()); + obj->setVal(Py_None); + PyDict_SetItemString(_context,varname.c_str(),Py_None); + obj->setVarName(varname); + return ; + } + if(obj->getKey()=="Bounds") + { + std::ostringstream oss; oss << "__" << _cnt++; + std::string varname(oss.str()); + obj->setVal(_bounds); + PyDict_SetItemString(_context,varname.c_str(),_bounds); + obj->setVarName(varname); + return ; + } + if(godFather->findPathOf(obj)=="Background/Vector") + { + std::ostringstream oss; oss << "__" << _cnt++; + std::string varname(oss.str()); + obj->setVal(_Xb); + PyDict_SetItemString(_context,varname.c_str(),_Xb); + obj->setVarName(varname); + } + if(obj->getKey()=="OneFunction") + { + std::ostringstream oss; oss << "__" << _cnt++; + std::string varname(oss.str()); + obj->setVal(_func); + PyDict_SetItemString(_context,varname.c_str(),_func); + obj->setVarName(varname); + return ; + } + if(godFather->findPathOf(obj)=="Observation/Vector") + { + std::ostringstream oss; oss << "__" << _cnt++; + std::string varname(oss.str()); + obj->setVal(_observation); + PyDict_SetItemString(_context,varname.c_str(),_observation); + obj->setVarName(varname); + } + } +private: + unsigned int _cnt = 0; + PyObject *_bounds = nullptr; + PyObject *_Xb = nullptr; + PyObject *_observation = nullptr; + PyObjectRAII _func; + PyObject *_context; +}; + +void AdaoExchangeLayer4Quintet::loadTemplate(AdaoModel::MainModel *model) +{ + const char DECORATOR_FUNC[]="def DecoratorAdao(cppFunc):\n" + " def evaluator( xserie ):\n" + " import numpy as np\n" + " yserie = [np.array(elt) for elt in cppFunc(xserie)]\n" + " return yserie\n" + " return evaluator\n"; + this->_internal->_py_call_back.assign(PyObject_GC_New(AdaoCallbackSt,&AdaoCallbackType), + &this->_internal->_data_btw_threads); + PyObject *callbackPyObj(this->_internal->_py_call_back.getPyObject()); + // + { + PyObjectRAII res(PyObjectRAII::FromNew(PyRun_String(DECORATOR_FUNC,Py_file_input,this->_internal->_context,this->_internal->_context))); + PyObjectRAII decoratorGenerator( PyObjectRAII::FromBorrowed(PyDict_GetItemString(this->_internal->_context,"DecoratorAdao")) ); + if(decoratorGenerator.isNull()) + throw AdaoExchangeLayerException("Fail to locate DecoratorAdao function !"); + PyObjectRAII args(PyObjectRAII::FromNew(PyTuple_New(1))); + { PyTuple_SetItem(args,0,callbackPyObj); Py_XINCREF(callbackPyObj); } + this->_internal->_decorator_func = PyObjectRAII::FromNew(PyObject_CallObject(decoratorGenerator,args)); + if(this->_internal->_decorator_func.isNull()) + throw AdaoExchangeLayerException("Fail to generate result of DecoratorAdao function !"); + } + // + Visitor1 visitor(this->_internal->_decorator_func,this->_internal->_context); + model->visitPythonLeaves(&visitor); + // + { + std::string sciptPyOfModelMaker(model->pyStr()); + PyObjectRAII res(PyObjectRAII::FromNew(PyRun_String(sciptPyOfModelMaker.c_str(),Py_file_input,this->_internal->_context,this->_internal->_context))); + _internal->_adao_case = PyObjectRAII::FromNew(PyDict_GetItemString(this->_internal->_context,"case")); + } + if(_internal->_adao_case.isNull()) + throw AdaoExchangeLayerException("Fail to generate ADAO case object !"); + // + _internal->_execute_func=PyObjectRAII::FromNew(PyObject_GetAttrString(_internal->_adao_case,"execute")); + if(_internal->_execute_func.isNull()) + throw AdaoExchangeLayerException("Fail to locate execute function of ADAO case object !"); +} + +void ExecuteAsync(PyObject *pyExecuteFunction, DataExchangedBetweenThreads *data) +{ + AutoGIL gil; // launched in a separed thread -> protect python calls + PyObjectRAII args(PyObjectRAII::FromNew(PyTuple_New(0))); + PyObjectRAII nullRes(PyObjectRAII::FromNew(PyObject_CallObject(pyExecuteFunction,args)));// go to adaocallback_call + PyErr_Print(); + data->_finished = true; + data->_data = nullptr; + sem_post(&data->_sem); +} + +void AdaoExchangeLayer4Quintet::execute() +{ + _internal->_tstate=PyEval_SaveThread(); // release the lock acquired in AdaoExchangeLayer4Quintet::initPythonIfNeeded by PyEval_InitThreads() + _internal->_fut = std::async(std::launch::async,ExecuteAsync,_internal->_execute_func,&_internal->_data_btw_threads); +} + +bool AdaoExchangeLayer4Quintet::next(PyObject *& inputRequested) +{ + sem_wait(&_internal->_data_btw_threads._sem); + if(_internal->_data_btw_threads._finished) + { + inputRequested = nullptr; + return false; + } + else + { + inputRequested = (PyObject *)_internal->_data_btw_threads._data; + return true; + } +} + +void AdaoExchangeLayer4Quintet::setResult(PyObject *outputAssociated) +{ + _internal->_data_btw_threads._data = outputAssociated; + _internal->_data_btw_threads._finished = false; + sem_post(&_internal->_data_btw_threads._sem_result_is_here); +} + +PyObject *AdaoExchangeLayer4Quintet::getResult() +{ + _internal->_fut.wait(); + PyEval_RestoreThread(_internal->_tstate); + AutoGIL gil; + // now retrieve case.get("Analysis")[-1] + PyObjectRAII get_func_of_adao_case(PyObjectRAII::FromNew(PyObject_GetAttrString(_internal->_adao_case,"get"))); + if(get_func_of_adao_case.isNull()) + throw AdaoExchangeLayerException("Fail to locate \"get\" method from ADAO case !"); + PyObjectRAII all_intermediate_results; + {// retrieve return data from case.get("Analysis") + PyObjectRAII args(PyObjectRAII::FromNew(PyTuple_New(1))); + PyTuple_SetItem(args,0,PyUnicode_FromString("Analysis")); + all_intermediate_results=PyObjectRAII::FromNew(PyObject_CallObject(get_func_of_adao_case,args)); + if(all_intermediate_results.isNull()) + throw AdaoExchangeLayerException("Fail to retrieve result of case.get(\"Analysis\") !"); + } + PyObjectRAII optimum; + { + PyObjectRAII param(PyObjectRAII::FromNew(PyLong_FromLong(-1))); + optimum=PyObjectRAII::FromNew(PyObject_GetItem(all_intermediate_results,param)); + if(optimum.isNull()) + throw AdaoExchangeLayerException("Fail to retrieve result of last element of case.get(\"Analysis\") !"); + } + /*PyObjectRAII code(PyObjectRAII::FromNew(Py_CompileString("case.get(\"Analysis\")[-1]","retrieve result",Py_file_input))); + if(code.isNull()) + throw AdaoExchangeLayerException("Fail to compile code to retrieve result after ADAO computation !"); + PyObjectRAII res(PyObjectRAII::FromNew(PyEval_EvalCode(code,_internal->_context,_internal->_context)));*/ + return optimum.retn(); +} diff --git a/AdaoExchangeLayer4Quintet.hxx b/AdaoExchangeLayer4Quintet.hxx new file mode 100644 index 0000000..bb11f33 --- /dev/null +++ b/AdaoExchangeLayer4Quintet.hxx @@ -0,0 +1,48 @@ +// Copyright (C) 2019 EDF R&D +// +// This library is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 2.1 of the License. +// +// This library is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// Lesser General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License along with this library; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +// +// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com +// +// Author: Anthony Geay, anthony.geay@edf.fr, EDF R&D + +#pragma once + +#include "Python.h" + +class AdaoCallbackSt; + +namespace AdaoModel +{ + class MainModel; +} + +class AdaoExchangeLayer4Quintet +{ + class Internal; +public: + AdaoExchangeLayer4Quintet(); + ~AdaoExchangeLayer4Quintet(); + void init(AdaoModel::MainModel *model); + void execute(); + bool next(PyObject *& inputRequested); + void setResult(PyObject *outputAssociated); + PyObject *getResult(); +private: + void initPythonIfNeeded(); + void loadTemplate(AdaoModel::MainModel *model); +private: + Internal *_internal = nullptr; +}; diff --git a/AdaoModelKeyVal.cxx b/AdaoModelKeyVal.cxx new file mode 100644 index 0000000..b38c56f --- /dev/null +++ b/AdaoModelKeyVal.cxx @@ -0,0 +1,439 @@ +// Copyright (C) 2019 EDF R&D +// +// This library is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 2.1 of the License. +// +// This library is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// Lesser General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License along with this library; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +// +// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com +// +// Author: Anthony Geay, anthony.geay@edf.fr, EDF R&D + +#include "AdaoModelKeyVal.hxx" +#include "AdaoExchangeLayerException.hxx" + +#include + +using namespace AdaoModel; + +const char CostDecrementTolerance::KEY[]="CostDecrementTolerance"; + +const char StringObserver::KEY[]="String"; + +const char InfoObserver::KEY[]="Info"; + +const char StoreSupplKeyVal::KEY[]="StoreSupplementaryCalculations"; + +const char *StoreSupplKeyVal::DFTL[]={"CostFunctionJAtCurrentOptimum","CostFunctionJoAtCurrentOptimum","CurrentOptimum","SimulatedObservationAtCurrentOptimum","SimulatedObservationAtOptimum"}; + +const char EnumAlgoKeyVal::KEY[]="Algorithm"; + +const char ParametersOfAlgorithmParameters::KEY[]="Parameters"; + +const char Bounds::KEY[]="Bounds"; + +const char MaximumNumberOfSteps::KEY[]="MaximumNumberOfSteps"; + +const char VectorBackground::KEY[]="Vector"; + +const char StoreBackground::KEY[]="Stored"; + +const char MatrixBackgroundError::KEY[]="Matrix"; + +const char ScalarSparseMatrixError::KEY[]="ScalarSparseMatrix"; + +const char DiagonalSparseMatrixError::KEY[]="DiagonalSparseMatrix"; + +const char OneFunction::KEY[]="OneFunction"; + +const char DifferentialIncrement::KEY[]="DifferentialIncrement"; + +const char ObservationOperatorParameters::KEY[]="Parameters"; + +const char CenteredFiniteDifference::KEY[]="CenteredFiniteDifference"; + +const char InputFunctionAsMulti::KEY[]="InputFunctionAsMulti"; + +const char VariableKV::KEY[]="Variable"; + +const char TemplateKV::KEY[]="Template"; + +const char AlgorithmParameters::KEY[]="AlgorithmParameters"; + +const char Background::KEY[]="Background"; + +const char BackgroundError::KEY[]="BackgroundError"; + +const double BackgroundError::BACKGROUND_SCALAR_SPARSE_DFT = 1.e10; + +const char ObservationError::KEY[]="ObservationError"; + +const char ObservationOperator::KEY[]="ObservationOperator"; + +const char Observation::KEY[]="Observation"; + +const double ObservationError::BACKGROUND_SCALAR_SPARSE_DFT = 1.; + +const char ObserverEntry::KEY[]="Observer"; + +std::string TopEntry::getParamForSet(const GenericKeyVal& entry) const +{ + std::ostringstream oss; + oss << "case.set(\'" << entry.getKey() << "\' , **" << entry.pyStr() << ")"; + return oss.str(); +} + +std::string GenericKeyVal::pyStrKeyVal() const +{ + std::ostringstream oss; + oss << "\"" << this->getKey() << "\" : " << this->pyStr(); + return oss.str(); +} + +void GenericKeyVal::visitAll(MainModel *godFather, RecursiveVisitor *visitor) +{ + visitor->visit(this); +} + +std::string DoubleKeyVal::pyStr() const +{ + std::ostringstream oss; + oss << std::scientific << _val; + return oss.str(); +} + +std::string BoolKeyVal::pyStr() const +{ + return _val?"True":"False"; +} + +std::string StringKeyVal::pyStr() const +{ + std::ostringstream oss; + oss << "\"" << _val << "\""; + return oss.str(); +} + +std::string NoneKeyVal::pyStr() const +{ + return "None"; +} + +std::string ListStringsKeyVal::pyStr() const +{ + std::ostringstream oss; + oss << "[ "; + std::size_t len(_val.size()); + for(std::size_t i=0;i EnumAlgoKeyVal::generateDftParameters() const +{ + switch(_enum) + { + case EnumAlgo::ThreeDVar: + case EnumAlgo::NonLinearLeastSquares: + { + return templateForOthers(); + } + case EnumAlgo::Blue: + { + return templateForBlue(); + } + default: + throw AdaoExchangeLayerException("EnumAlgoKeyVal::generateDftParameters : Unrecognized Algo !"); + } +} + +std::string EnumAlgoKeyVal::getRepresentation() const +{ + switch(_enum) + { + case EnumAlgo::ThreeDVar: + return "3DVAR"; + case EnumAlgo::NonLinearLeastSquares: + return "NonLinearLeastSquares"; + case EnumAlgo::Blue: + return "Blue"; + default: + throw AdaoExchangeLayerException("EnumAlgoKeyVal::getRepresentation : Unrecognized Algo !"); + } +} + +std::string EnumAlgoKeyVal::pyStr() const +{ + std::ostringstream oss; + oss << "\"" << this->getRepresentation() << "\""; + return oss.str(); +} + +std::shared_ptr EnumAlgoKeyVal::templateForBlue() const +{ + std::shared_ptr ret(std::make_shared()); + std::shared_ptr v(std::make_shared()); + ret->pushBack(std::static_pointer_cast(v)); + return ret; +} + +std::shared_ptr EnumAlgoKeyVal::templateForOthers() const +{ + std::shared_ptr ret(std::make_shared()); + std::shared_ptr v0(std::make_shared()); + std::shared_ptr v1(std::make_shared()); + std::shared_ptr v2(std::make_shared()); + std::shared_ptr v3(std::make_shared()); + ret->pushBack(std::static_pointer_cast(v0)); + ret->pushBack(std::static_pointer_cast(v1)); + ret->pushBack(std::static_pointer_cast(v2)); + ret->pushBack(std::static_pointer_cast(v3)); + return ret; +} + +std::string DictKeyVal::pyStr() const +{ + std::ostringstream oss; + oss << "{ "; + std::size_t len(_pairs.size()); + for(std::size_t i=0;ipyStrKeyVal(); + if(i!=len-1) + oss << ", "; + } + oss << " }"; + return oss.str(); +} + +void DictKeyVal::visitPython(MainModel *godFather, PythonLeafVisitor *visitor) +{ + for(auto elt : _pairs) + { + elt->visitPython(godFather, visitor); + } +} + +void DictKeyVal::visitAll(MainModel *godFather, RecursiveVisitor *visitor) +{ + visitor->enterSubDir(this); + for(auto elt : _pairs) + { + elt->visitAll(godFather, visitor); + } + visitor->exitSubDir(this); +} + +std::string PyObjKeyVal::pyStr() const +{ + std::ostringstream oss; + oss << _var_name; + return oss.str(); +} + +void PyObjKeyVal::visitPython(MainModel *godFather, PythonLeafVisitor *visitor) +{ + visitor->visit(godFather,this); +} + +std::string UnsignedIntKeyVal::pyStr() const +{ + std::ostringstream oss; + oss << _val; + return oss.str(); +} + +void InputFunctionAsMulti::setVal(bool val) +{ + if(!val) + throw AdaoExchangeLayerException("InputFunctionAsMulti : value has to remain to true !"); +} + +ObservationOperatorParameters::ObservationOperatorParameters():DictKeyVal(KEY) +{ + std::shared_ptr v0(std::make_shared()); + std::shared_ptr v1(std::make_shared()); + _pairs.push_back(std::static_pointer_cast(v0)); + _pairs.push_back(std::static_pointer_cast(v1)); +} + +AlgorithmParameters::AlgorithmParameters():DictKeyVal(KEY) +{ + std::shared_ptr v0(std::make_shared()); + std::shared_ptr v1(v0->generateDftParameters()); + _pairs.push_back(std::static_pointer_cast(v0)); + _pairs.push_back(std::static_pointer_cast(v1)); +} + +Background::Background():DictKeyVal(KEY) +{ + std::shared_ptr v0(std::make_shared()); + std::shared_ptr v1(std::make_shared()); + v1->setVal(true); + _pairs.push_back(std::static_pointer_cast(v0)); + _pairs.push_back(std::static_pointer_cast(v1)); +} + +GenericError::GenericError(const std::string& key, double dftValForScalarSparseMatrix):DictKeyVal(key) +{ + std::shared_ptr v0(std::make_shared()); + std::shared_ptr v1(std::make_shared(dftValForScalarSparseMatrix)); + std::shared_ptr v2(std::make_shared()); + _pairs.push_back(std::static_pointer_cast(v0)); + _pairs.push_back(std::static_pointer_cast(v1)); + _pairs.push_back(std::static_pointer_cast(v2)); +} + +Observation::Observation():DictKeyVal(KEY) +{ + std::shared_ptr v0(std::make_shared()); + std::shared_ptr v1(std::make_shared()); + v1->setVal(false); + _pairs.push_back(std::static_pointer_cast(v0)); + _pairs.push_back(std::static_pointer_cast(v1)); +} + +ObservationOperator::ObservationOperator():DictKeyVal(KEY) +{ + std::shared_ptr v0(std::make_shared()); + std::shared_ptr v1(std::make_shared()); + std::shared_ptr v2(std::make_shared()); + _pairs.push_back(std::static_pointer_cast(v0)); + _pairs.push_back(std::static_pointer_cast(v1)); + _pairs.push_back(std::static_pointer_cast(v2)); +} + +ObserverEntry::ObserverEntry():DictKeyVal(KEY) +{ + std::shared_ptr v0(std::make_shared()); + std::shared_ptr v1(std::make_shared()); + std::shared_ptr v2(std::make_shared()); + std::shared_ptr v3(std::make_shared()); + _pairs.push_back(std::static_pointer_cast(v0)); + _pairs.push_back(std::static_pointer_cast(v1)); + _pairs.push_back(std::static_pointer_cast(v2)); + _pairs.push_back(std::static_pointer_cast(v3)); +} + +MainModel::MainModel():_algo(std::make_shared()), + _bg(std::make_shared()), + _bg_err(std::make_shared()), + _obs(std::make_shared()), + _obs_err(std::make_shared()), + _observ_op(std::make_shared()), + _observ_entry(std::make_shared()) +{ +} + +std::string MainModel::pyStr() const +{ + std::ostringstream oss; + oss << "from adao import adaoBuilder" << std::endl; + oss << "case = adaoBuilder.New()" << std::endl; + oss << _algo->getParamForSet(*_algo) << std::endl; + oss << _bg->getParamForSet(*_bg) << std::endl; + oss << _bg_err->getParamForSet(*_bg_err) << std::endl; + oss << _obs->getParamForSet(*_obs) << std::endl; + oss << _obs_err->getParamForSet(*_obs_err) << std::endl; + oss << _observ_op->getParamForSet(*_observ_op) << std::endl; + oss << _observ_entry->getParamForSet(*_observ_entry) << std::endl; + return oss.str(); +} + +std::vector< std::shared_ptr > MainModel::toVect() const +{ + return { + std::static_pointer_cast(_algo), + std::static_pointer_cast(_bg), + std::static_pointer_cast(_bg_err), + std::static_pointer_cast(_obs), + std::static_pointer_cast(_obs_err), + std::static_pointer_cast(_observ_op), + std::static_pointer_cast(_observ_entry) + }; +} + +class MyAllVisitor : public RecursiveVisitor +{ +public: + MyAllVisitor(GenericKeyVal *elt):_elt_to_find(elt) { } + void visit(GenericKeyVal *elt) override + { + if(_found) return ; + if(elt != _elt_to_find) return ; + _path.push_back(elt->getKey()); + _found=true; + } + void enterSubDir(DictKeyVal *subdir) override + { + if(_found) return ; + _path.push_back(subdir->getKey()); + } + void exitSubDir(DictKeyVal *subdir) override + { + if(_found) return ; + _path.pop_back(); + } + std::string getPath() const + { + std::ostringstream oss; + std::size_t len(_path.size()),ii(0); + for(auto elt : _path) + { + oss << elt; + if(ii != len-1) + oss << "/"; + ii++; + } + return oss.str(); + } +private: + GenericKeyVal *_elt_to_find = nullptr; + std::vector _path; + bool _found = false; +}; + +std::string MainModel::findPathOf(GenericKeyVal *elt) +{ + MyAllVisitor vis(elt); + this->visitAll(&vis); + return vis.getPath(); +} + +void MainModel::visitPythonLeaves(PythonLeafVisitor *visitor) +{ + std::vector< std::shared_ptr > sons(toVect()); + for(auto elt : sons) + { + elt->visitPython(this, visitor); + } +} + +void MainModel::visitAll(RecursiveVisitor *visitor) +{ + std::vector< std::shared_ptr > sons(toVect()); + for(auto elt : sons) + { + elt->visitAll(this, visitor); + } +} diff --git a/AdaoModelKeyVal.hxx b/AdaoModelKeyVal.hxx new file mode 100644 index 0000000..f9c8132 --- /dev/null +++ b/AdaoModelKeyVal.hxx @@ -0,0 +1,466 @@ +// Copyright (C) 2019 EDF R&D +// +// This library is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 2.1 of the License. +// +// This library is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// Lesser General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License along with this library; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +// +// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com +// +// Author: Anthony Geay, anthony.geay@edf.fr, EDF R&D + +#pragma once + +#include "PyObjectRAII.hxx" + +#include +#include +#include + +namespace AdaoModel +{ + enum class Type + { + Double, + UnsignedInt, + Bool, + String, + PyObj, + EnumAlgo, + ListStrings, + None, + Child + }; + + enum class EnumAlgo + { + ThreeDVar, + Blue, + NonLinearLeastSquares + }; + + class GenericKeyVal; + class MainModel; + + class TopEntry + { + public: + std::string getParamForSet(const GenericKeyVal& entry) const; + }; + + class PythonLeafVisitor; + class RecursiveVisitor; + + class GenericKeyVal + { + public: + virtual bool isNeeded() const = 0; + virtual Type getType() const = 0; + virtual std::string pyStr() const = 0; + virtual void visitPython(MainModel *godFather, PythonLeafVisitor *visitor) { } + virtual void visitAll(MainModel *godFather, RecursiveVisitor *visitor); + virtual ~GenericKeyVal() { } + std::string pyStrKeyVal() const; + std::string getKey() const { return _key; } + protected: + GenericKeyVal(const std::string& key):_key(key) { } + private: + std::string _key; + }; + + class NeededGenericKeyVal : public GenericKeyVal + { + protected: + NeededGenericKeyVal(const std::string& key):GenericKeyVal(key) { } + public: + bool isNeeded() const override { return true; } + }; + + class NotNeededGenericKeyVal : public GenericKeyVal + { + public: + bool isNeeded() const override { return false; } + }; + + class DoubleKeyVal : public NeededGenericKeyVal + { + public: + DoubleKeyVal(const std::string& key):NeededGenericKeyVal(key) { } + Type getType() const override { return Type::Double; } + void setVal(double val) { _val=val; } + double getVal() const { return _val; } + std::string pyStr() const override; + private: + double _val = 0.; + }; + + class BoolKeyVal : public NeededGenericKeyVal + { + public: + BoolKeyVal(const std::string& key):NeededGenericKeyVal(key) { } + Type getType() const override { return Type::Bool; } + virtual void setVal(bool val) { _val=val; } + bool getVal() const { return _val; } + std::string pyStr() const override; + private: + bool _val = false; + }; + + class CostDecrementTolerance : public DoubleKeyVal + { + public: + CostDecrementTolerance():DoubleKeyVal(KEY) { setVal(1e-7); } + public: + static const char KEY[]; + }; + + class StringKeyVal : public NeededGenericKeyVal + { + public: + StringKeyVal(const std::string& key):NeededGenericKeyVal(key) { } + Type getType() const override { return Type::String; } + std::string pyStr() const override; + void setVal(const std::string& val) { _val = val; } + std::string getVal() const { return _val; } + private: + std::string _val; + }; + + class NoneKeyVal : public NeededGenericKeyVal + { + public: + NoneKeyVal(const std::string& key):NeededGenericKeyVal(key) { } + Type getType() const override { return Type::None; } + std::string pyStr() const override; + }; + + class StringObserver : public NoneKeyVal + { + public: + StringObserver():NoneKeyVal(KEY) { } + public: + static const char KEY[]; + }; + + class InfoObserver : public NoneKeyVal + { + public: + InfoObserver():NoneKeyVal(KEY) { } + public: + static const char KEY[]; + }; + + class ListStringsKeyVal : public NeededGenericKeyVal + { + protected: + ListStringsKeyVal(const std::string& key):NeededGenericKeyVal(key) { } + public: + Type getType() const override { return Type::ListStrings; } + std::string pyStr() const override; + protected: + std::vector< std::string > _val; + }; + + class StoreSupplKeyVal : public ListStringsKeyVal + { + public: + StoreSupplKeyVal(); + public: + static const char *DFTL[]; + static const char KEY[]; + }; + + class DictKeyVal; + + class EnumAlgoKeyVal : public NeededGenericKeyVal + { + public: + EnumAlgoKeyVal():NeededGenericKeyVal(KEY),_enum(EnumAlgo::ThreeDVar) { } + Type getType() const override { return Type::EnumAlgo; } + std::shared_ptr generateDftParameters() const; + std::string getRepresentation() const; + std::string pyStr() const override; + void setVal(EnumAlgo newAlgo) { _enum = newAlgo; } + EnumAlgo getVal() const { return _enum; } + private: + std::shared_ptr templateForBlue() const; + std::shared_ptr templateForOthers() const; + private: + static const char KEY[]; + EnumAlgo _enum; + }; + + class DictKeyVal : public NeededGenericKeyVal + { + public: + DictKeyVal(const std::string& key):NeededGenericKeyVal(key) { } + void pushBack(std::shared_ptr elt) { _pairs.push_back(elt); } + Type getType() const override { return Type::Child; } + std::string pyStr() const override; + void visitPython(MainModel *godFather, PythonLeafVisitor *visitor) override; + void visitAll(MainModel *godFather, RecursiveVisitor *visitor) override; + protected: + std::vector< std::shared_ptr > _pairs; + }; + + class ParametersOfAlgorithmParameters : public DictKeyVal + { + public: + ParametersOfAlgorithmParameters():DictKeyVal(KEY) { } + public: + static const char KEY[]; + }; + + class PyObjKeyVal : public NeededGenericKeyVal + { + public: + PyObjKeyVal(const std::string& key):NeededGenericKeyVal(key) { } + Type getType() const override { return Type::PyObj; } + void setVal(PyObject *obj) { _val = PyObjectRAII::FromBorrowed(obj); } + // borrowed ref + PyObject *getVal() const { return _val; } + std::string pyStr() const; + void setVarName(const std::string& vn) { _var_name = vn; } + void visitPython(MainModel *godFather, PythonLeafVisitor *visitor) override; + private: + PyObjectRAII _val; + std::string _var_name; + }; + + class Bounds : public PyObjKeyVal + { + public: + Bounds():PyObjKeyVal(KEY) { } + public: + static const char KEY[]; + }; + + class UnsignedIntKeyVal : public NeededGenericKeyVal + { + public: + UnsignedIntKeyVal(const std::string& key):NeededGenericKeyVal(key) { } + void setVal(unsigned int val) { _val=val; } + unsigned int getVal() const { return _val; } + Type getType() const override { return Type::UnsignedInt; } + std::string pyStr() const; + private: + unsigned int _val = -1; + }; + + class MaximumNumberOfSteps : public UnsignedIntKeyVal + { + public: + MaximumNumberOfSteps():UnsignedIntKeyVal(KEY) { setVal(100); } + public: + static const char KEY[]; + }; + + class VectorBackground : public PyObjKeyVal + { + public: + VectorBackground():PyObjKeyVal(KEY) { } + public: + static const char KEY[]; + }; + + class StoreBackground : public BoolKeyVal + { + public: + StoreBackground():BoolKeyVal(KEY) { } + public: + static const char KEY[]; + }; + + class MatrixBackgroundError : public PyObjKeyVal + { + public: + MatrixBackgroundError():PyObjKeyVal(KEY) { } + public: + static const char KEY[]; + }; + + class ScalarSparseMatrixError : public DoubleKeyVal + { + public: + ScalarSparseMatrixError(double val):DoubleKeyVal(KEY) { setVal(val); } + public: + static const char KEY[]; + }; + + class DiagonalSparseMatrixError : public PyObjKeyVal + { + public: + DiagonalSparseMatrixError():PyObjKeyVal(KEY) { } + public: + static const char KEY[]; + }; + + class OneFunction : public PyObjKeyVal + { + public: + OneFunction():PyObjKeyVal(KEY) { } + public: + static const char KEY[]; + }; + + class DifferentialIncrement : public DoubleKeyVal + { + public: + DifferentialIncrement():DoubleKeyVal(KEY) { setVal(0.0001); } + public: + static const char KEY[]; + }; + + class ObservationOperatorParameters : public DictKeyVal + { + public: + ObservationOperatorParameters(); + public: + static const char KEY[]; + }; + + class CenteredFiniteDifference : public BoolKeyVal + { + public: + CenteredFiniteDifference():BoolKeyVal(KEY) { setVal(false); } + public: + static const char KEY[]; + }; + + class InputFunctionAsMulti : public BoolKeyVal + { + public: + InputFunctionAsMulti():BoolKeyVal(KEY) { BoolKeyVal::setVal(true); } + void setVal(bool val) override; + public: + static const char KEY[]; + }; + + class VariableKV : public StringKeyVal + { + public: + VariableKV():StringKeyVal(KEY) { setVal("CurrentState"); } + public: + static const char KEY[]; + }; + + class TemplateKV : public StringKeyVal + { + public: + TemplateKV():StringKeyVal(KEY) { setVal("CurrentState"); } + public: + static const char KEY[]; + }; + + ////////////// + + class AlgorithmParameters : public DictKeyVal, public TopEntry + { + public: + AlgorithmParameters(); + public: + static const char KEY[]; + }; + + class Background : public DictKeyVal, public TopEntry + { + public: + Background(); + public: + static const char KEY[]; + }; + + class GenericError : public DictKeyVal, public TopEntry + { + protected: + GenericError(const std::string& key, double dftValForScalarSparseMatrix); + }; + + class BackgroundError : public GenericError + { + public: + BackgroundError():GenericError(KEY,BACKGROUND_SCALAR_SPARSE_DFT) { } + public: + static const char KEY[]; + static const double BACKGROUND_SCALAR_SPARSE_DFT; + }; + + class ObservationError : public GenericError + { + public: + ObservationError():GenericError(KEY,BACKGROUND_SCALAR_SPARSE_DFT) { } + public: + static const char KEY[]; + static const double BACKGROUND_SCALAR_SPARSE_DFT; + }; + + class ObservationOperator : public DictKeyVal, public TopEntry + { + public: + ObservationOperator(); + public: + static const char KEY[]; + }; + + class Observation : public DictKeyVal, public TopEntry + { + public: + Observation(); + public: + static const char KEY[]; + }; + + class ObserverEntry : public DictKeyVal, public TopEntry + { + public: + ObserverEntry(); + public: + static const char KEY[]; + }; + + /////////////// + + class PythonLeafVisitor + { + public: + virtual void visit(MainModel *godFather, PyObjKeyVal *obj) = 0; + virtual ~PythonLeafVisitor() { } + }; + + class RecursiveVisitor + { + public: + virtual ~RecursiveVisitor() { } + virtual void visit(GenericKeyVal *obj) = 0; + virtual void enterSubDir(DictKeyVal *subdir) = 0; + virtual void exitSubDir(DictKeyVal *subdir) = 0; + }; + + class MainModel + { + public: + MainModel(); + std::string findPathOf(GenericKeyVal *elt); + std::string pyStr() const; + std::vector< std::shared_ptr > toVect() const; + void visitPythonLeaves(PythonLeafVisitor *visitor); + void visitAll(RecursiveVisitor *visitor); + private: + std::shared_ptr _algo; + std::shared_ptr _bg; + std::shared_ptr _bg_err; + std::shared_ptr _obs; + std::shared_ptr _obs_err; + std::shared_ptr _observ_op; + std::shared_ptr _observ_entry; + }; +} diff --git a/CMakeLists.txt b/CMakeLists.txt new file mode 100644 index 0000000..cedb667 --- /dev/null +++ b/CMakeLists.txt @@ -0,0 +1,76 @@ +# Copyright (C) 2019 EDF R&D +# +# This library is free software; you can redistribute it and/or +# modify it under the terms of the GNU Lesser General Public +# License as published by the Free Software Foundation; either +# version 2.1 of the License. +# +# This library is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +# Lesser General Public License for more details. +# +# You should have received a copy of the GNU Lesser General Public +# License along with this library; if not, write to the Free Software +# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +# +# See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com +# +# Author: Anthony Geay, anthony.geay@edf.fr, EDF R&D + +cmake_minimum_required(VERSION 2.6) +project(AdaoExchangeLayer4Quintet) + +## + +set(BUILD_SHARED_LIBS TRUE) +set(CONFIGURATION_ROOT_DIR $ENV{CONFIGURATION_ROOT_DIR} CACHE PATH "Path to the Salome CMake configuration files") +if(EXISTS ${CONFIGURATION_ROOT_DIR}) + list(APPEND CMAKE_MODULE_PATH "${CONFIGURATION_ROOT_DIR}/cmake") + include(SalomeMacros) + include(SalomeSetupPlatform) +else(EXISTS ${CONFIGURATION_ROOT_DIR}) + message(FATAL_ERROR "We absolutely need the Salome CMake configuration files, please define CONFIGURATION_ROOT_DIR !") +endif(EXISTS ${CONFIGURATION_ROOT_DIR}) + +## + +option(AEL_ENABLE_TESTS "Build tests (default ON)." ON) +if(AEL_ENABLE_TESTS) + if(EXISTS ${PY2CPP_ROOT_DIR}) + set(PY2CPP_ROOT_DIR $ENV{PY2CPP_ROOT_DIR} CACHE PATH "Path to Py2cpp") + list(APPEND CMAKE_MODULE_PATH "${PY2CPP_ROOT_DIR}/lib/cmake/py2cpp") + find_package(Py2cpp REQUIRED) + get_target_property(PY2CPP_INCLUDE_DIR py2cpp INTERFACE_INCLUDE_DIRECTORIES) + get_target_property(PY2CPP_LIB py2cpp IMPORTED_LOCATION_RELEASE) + find_package(SalomeCppUnit) + else(EXISTS ${PY2CPP_ROOT_DIR}) + message(FATAL_ERROR "PY2CPP declared as enabled whereas not found !") + endif(EXISTS ${PY2CPP_ROOT_DIR}) +endif(AEL_ENABLE_TESTS) + +## + +find_package(SalomePythonInterp REQUIRED) +find_package(SalomePythonLibs REQUIRED) + +## + +include_directories( + ${PYTHON_INCLUDE_DIRS} + ${PY2CPP_INCLUDE_DIR} + ${CPPUNIT_INCLUDE_DIRS} + ) +set(adaoexchange_SOURCES AdaoExchangeLayer4Quintet.cxx AdaoModelKeyVal.cxx) +add_library(adaoexchange ${adaoexchange_SOURCES}) +target_link_libraries(adaoexchange ${PYTHON_LIBRARIES}) +install(FILES AdaoExchangeLayer4Quintet.hxx PyObjectRAII.hxx AdaoExchangeLayerException.hxx AdaoModelKeyVal.hxx DESTINATION include) +install(TARGETS adaoexchange DESTINATION lib) + +## + +if(AEL_ENABLE_TESTS) + add_executable(TestAdaoExchange TestAdaoExchange.cxx) + target_link_libraries(TestAdaoExchange adaoexchange ${PY2CPP_LIB} ${CPPUNIT_LIBRARIES}) + install(TARGETS TestAdaoExchange DESTINATION bin) +endif(AEL_ENABLE_TESTS) diff --git a/Cases/Definition_complete_de_cas_3DVAR.py b/Cases/Definition_complete_de_cas_3DVAR.py new file mode 100644 index 0000000..35d3135 --- /dev/null +++ b/Cases/Definition_complete_de_cas_3DVAR.py @@ -0,0 +1,169 @@ +# -*- coding: utf-8 -*- +# +# Copyright (C) 2008-2019 EDF R&D +# +# This library is free software; you can redistribute it and/or +# modify it under the terms of the GNU Lesser General Public +# License as published by the Free Software Foundation; either +# version 2.1 of the License. +# +# This library is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +# Lesser General Public License for more details. +# +# You should have received a copy of the GNU Lesser General Public +# License along with this library; if not, write to the Free Software +# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +# +# See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com +# +# Author: Jean-Philippe Argaud, jean-philippe.argaud@edf.fr, EDF R&D +"Full definition of a use case for the standard user" + +import sys +import unittest +import numpy + +# ============================================================================== +# +# Artificial user data example +# ---------------------------- +alpha = 5. +beta = 7 +gamma = 9.0 +# +alphamin, alphamax = 0., 10. +betamin, betamax = 3, 13 +gammamin, gammamax = 1.5, 15.5 +# +def simulation(x): + "Observation operator H for Y=H(X)" + import numpy + __x = numpy.ravel(x) + __H = numpy.array([[1, 0, 0], [0, 2, 0], [0, 0, 3], [1, 2, 3]]) + return numpy.dot(__H,__x) +# +def multisimulation( xserie ): + yserie = [] + for x in xserie: + yserie.append( simulation( x ) ) + return yserie +# +# Twin experiment observations +# ---------------------------- +observations = simulation((2, 3, 4)) + +# ============================================================================== +class InTest(unittest.TestCase): + def test1(self): + print(""" + Full definition of a use case for the standard user + +++++++++++++++++++++++++++++++++++++++++++++++++++ + """) + # + import numpy + from adao import adaoBuilder + # + # Mise en forme des entrees + # ------------------------- + Xb = (alpha, beta, gamma) + Bounds = ( + (alphamin, alphamax), + (betamin, betamax ), + (gammamin, gammamax)) + # + # TUI ADAO + # -------- + case = adaoBuilder.New() + case.set( 'AlgorithmParameters', + Algorithm = '3DVAR', # Mots-clé réservé + Parameters = { # Dictionnaire + "Bounds":Bounds, # Liste de paires de Real ou de None + "MaximumNumberOfSteps":100, # Int >= 0 + "CostDecrementTolerance":1.e-7, # Real > 0 + "StoreSupplementaryCalculations":[# Liste de mots-clés réservés + "CostFunctionJAtCurrentOptimum", + "CostFunctionJoAtCurrentOptimum", + "CurrentOptimum", + "SimulatedObservationAtCurrentOptimum", + "SimulatedObservationAtOptimum", + ], + } + ) + case.set( 'Background', + Vector = Xb, # array, list, tuple, matrix + Stored = True, # Bool + ) + case.set( 'Observation', + Vector = observations, # array, list, tuple, matrix + Stored = False, # Bool + ) + case.set( 'BackgroundError', + Matrix = None, # None ou matrice carrée + ScalarSparseMatrix = 1.0e10, # None ou Real > 0 + DiagonalSparseMatrix = None, # None ou vecteur + ) + case.set( 'ObservationError', + Matrix = None, # None ou matrice carrée + ScalarSparseMatrix = 1.0, # None ou Real > 0 + DiagonalSparseMatrix = None, # None ou vecteur + ) + case.set( 'ObservationOperator', + OneFunction = multisimulation, # MultiFonction [Y] = F([X]) + Parameters = { # Dictionnaire + "DifferentialIncrement":0.0001, # Real > 0 + "CenteredFiniteDifference":False, # Bool + }, + InputFunctionAsMulti = True, # Bool + ) + case.set( 'Observer', + Variable = "CurrentState", # Mot-clé + Template = "ValuePrinter", # Mot-clé + String = None, # None ou code Python + Info = None, # None ou string + + ) + case.execute() + # + # Exploitation independante + # ------------------------- + Xbackground = case.get("Background") + Xoptimum = case.get("Analysis")[-1] + FX_at_optimum = case.get("SimulatedObservationAtOptimum")[-1] + J_values = case.get("CostFunctionJAtCurrentOptimum")[:] + print("") + print("Number of internal iterations...: %i"%len(J_values)) + print("Initial state...................: %s"%(numpy.ravel(Xbackground),)) + print("Optimal state...................: %s"%(numpy.ravel(Xoptimum),)) + print("Simulation at optimal state.....: %s"%(numpy.ravel(FX_at_optimum),)) + print("") + # + # Fin du cas + # ---------- + ecart = assertAlmostEqualArrays(Xoptimum, [ 2., 3., 4.]) + # + print(" The maximal absolute error in the test is of %.2e."%ecart) + print(" The results are correct.") + print("") + # + return Xoptimum + +# ============================================================================== +def assertAlmostEqualArrays(first, second, places=7, msg=None, delta=None): + "Compare two vectors, like unittest.assertAlmostEqual" + import numpy + if msg is not None: + print(msg) + if delta is not None: + if ( (numpy.asarray(first) - numpy.asarray(second)) > float(delta) ).any(): + raise AssertionError("%s != %s within %s places"%(first,second,delta)) + else: + if ( (numpy.asarray(first) - numpy.asarray(second)) > 10**(-int(places)) ).any(): + raise AssertionError("%s != %s within %i places"%(first,second,places)) + return max(abs(numpy.asarray(first) - numpy.asarray(second))) + +# ============================================================================== +if __name__ == '__main__': + print("\nAUTODIAGNOSTIC\n==============") + unittest.main() diff --git a/Cases/Definition_complete_de_cas_Blue.py b/Cases/Definition_complete_de_cas_Blue.py new file mode 100644 index 0000000..43c96da --- /dev/null +++ b/Cases/Definition_complete_de_cas_Blue.py @@ -0,0 +1,166 @@ +# -*- coding: utf-8 -*- +# +# Copyright (C) 2008-2019 EDF R&D +# +# This library is free software; you can redistribute it and/or +# modify it under the terms of the GNU Lesser General Public +# License as published by the Free Software Foundation; either +# version 2.1 of the License. +# +# This library is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +# Lesser General Public License for more details. +# +# You should have received a copy of the GNU Lesser General Public +# License along with this library; if not, write to the Free Software +# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +# +# See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com +# +# Author: Jean-Philippe Argaud, jean-philippe.argaud@edf.fr, EDF R&D +"Full definition of a use case for the standard user" + +import sys +import unittest +import numpy + +# ============================================================================== +# +# Artificial user data example +# ---------------------------- +alpha = 5. +beta = 7 +gamma = 9.0 +# +alphamin, alphamax = 0., 10. +betamin, betamax = 3, 13 +gammamin, gammamax = 1.5, 15.5 +# +def simulation(x): + "Observation operator H for Y=H(X)" + import numpy + __x = numpy.ravel(x) + __H = numpy.array([[1, 0, 0], [0, 2, 0], [0, 0, 3], [1, 2, 3]]) + return numpy.dot(__H,__x) +# +def multisimulation( xserie ): + yserie = [] + for x in xserie: + yserie.append( simulation( x ) ) + return yserie +# +# Twin experiment observations +# ---------------------------- +observations = simulation((2, 3, 4)) + +# ============================================================================== +class InTest(unittest.TestCase): + def test1(self): + print(""" + Full definition of a use case for the standard user + +++++++++++++++++++++++++++++++++++++++++++++++++++ + """) + # + import numpy + from adao import adaoBuilder + # + # Mise en forme des entrees + # ------------------------- + Xb = (alpha, beta, gamma) + Bounds = ( + (alphamin, alphamax), + (betamin, betamax ), + (gammamin, gammamax)) + # + # TUI ADAO + # -------- + case = adaoBuilder.New() + case.set( 'AlgorithmParameters', + Algorithm = 'Blue', # Mots-clé réservé + Parameters = { # Dictionnaire + "StoreSupplementaryCalculations":[# Liste de mots-clés réservés + "CostFunctionJAtCurrentOptimum", + "CostFunctionJoAtCurrentOptimum", + "CurrentOptimum", + "SimulatedObservationAtCurrentOptimum", + "SimulatedObservationAtOptimum", + ], + } + ) + case.set( 'Background', + Vector = Xb, # array, list, tuple, matrix + Stored = True, # Bool + ) + case.set( 'Observation', + Vector = observations, # array, list, tuple, matrix + Stored = False, # Bool + ) + case.set( 'BackgroundError', + Matrix = None, # None ou matrice carrée + ScalarSparseMatrix = 1.0e10, # None ou Real > 0 + DiagonalSparseMatrix = None, # None ou vecteur + ) + case.set( 'ObservationError', + Matrix = None, # None ou matrice carrée + ScalarSparseMatrix = 1.0, # None ou Real > 0 + DiagonalSparseMatrix = None, # None ou vecteur + ) + case.set( 'ObservationOperator', + OneFunction = multisimulation, # MultiFonction [Y] = F([X]) + Parameters = { # Dictionnaire + "DifferentialIncrement":0.0001, # Real > 0 + "CenteredFiniteDifference":False, # Bool + }, + InputFunctionAsMulti = True, # Bool + ) + case.set( 'Observer', + Variable = "CurrentState", # Mot-clé + Template = "ValuePrinter", # Mot-clé + String = None, # None ou code Python + Info = None, # None ou string + + ) + case.execute() + # + # Exploitation independante + # ------------------------- + Xbackground = case.get("Background") + Xoptimum = case.get("Analysis")[-1] + FX_at_optimum = case.get("SimulatedObservationAtOptimum")[-1] + J_values = case.get("CostFunctionJAtCurrentOptimum")[:] + print("") + print("Number of internal iterations...: %i"%len(J_values)) + print("Initial state...................: %s"%(numpy.ravel(Xbackground),)) + print("Optimal state...................: %s"%(numpy.ravel(Xoptimum),)) + print("Simulation at optimal state.....: %s"%(numpy.ravel(FX_at_optimum),)) + print("") + # + # Fin du cas + # ---------- + ecart = assertAlmostEqualArrays(Xoptimum, [ 2., 3., 4.]) + # + print("The maximal absolute error in the test is of %.2e."%ecart) + print("The results are correct.") + print("") + # + return Xoptimum + +# ============================================================================== +def assertAlmostEqualArrays(first, second, places=7, msg=None, delta=None): + "Compare two vectors, like unittest.assertAlmostEqual" + import numpy + if msg is not None: + print(msg) + if delta is not None: + if ( (numpy.asarray(first) - numpy.asarray(second)) > float(delta) ).any(): + raise AssertionError("%s != %s within %s places"%(first,second,delta)) + else: + if ( (numpy.asarray(first) - numpy.asarray(second)) > 10**(-int(places)) ).any(): + raise AssertionError("%s != %s within %i places"%(first,second,places)) + return max(abs(numpy.asarray(first) - numpy.asarray(second))) + +# ============================================================================== +if __name__ == '__main__': + print("\nAUTODIAGNOSTIC\n==============") + unittest.main() diff --git a/Cases/Definition_complete_de_cas_NLLS.py b/Cases/Definition_complete_de_cas_NLLS.py new file mode 100644 index 0000000..033e355 --- /dev/null +++ b/Cases/Definition_complete_de_cas_NLLS.py @@ -0,0 +1,164 @@ +# -*- coding: utf-8 -*- +# +# Copyright (C) 2008-2019 EDF R&D +# +# This library is free software; you can redistribute it and/or +# modify it under the terms of the GNU Lesser General Public +# License as published by the Free Software Foundation; either +# version 2.1 of the License. +# +# This library is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +# Lesser General Public License for more details. +# +# You should have received a copy of the GNU Lesser General Public +# License along with this library; if not, write to the Free Software +# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +# +# See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com +# +# Author: Jean-Philippe Argaud, jean-philippe.argaud@edf.fr, EDF R&D +"Full definition of a use case for the standard user" + +import sys +import unittest +import numpy + +# ============================================================================== +# +# Artificial user data example +# ---------------------------- +alpha = 5. +beta = 7 +gamma = 9.0 +# +alphamin, alphamax = 0., 10. +betamin, betamax = 3, 13 +gammamin, gammamax = 1.5, 15.5 +# +def simulation(x): + "Observation operator H for Y=H(X)" + import numpy + __x = numpy.ravel(x) + __H = numpy.array([[1, 0, 0], [0, 2, 0], [0, 0, 3], [1, 2, 3]]) + return numpy.dot(__H,__x) +# +def multisimulation( xserie ): + yserie = [] + for x in xserie: + yserie.append( simulation( x ) ) + return yserie +# +# Twin experiment observations +# ---------------------------- +observations = simulation((2, 3, 4)) + +# ============================================================================== +class InTest(unittest.TestCase): + def test1(self): + print(""" + Full definition of a use case for the standard user + +++++++++++++++++++++++++++++++++++++++++++++++++++ + """) + # + import numpy + from adao import adaoBuilder + # + # Mise en forme des entrees + # ------------------------- + Xb = (alpha, beta, gamma) + Bounds = ( + (alphamin, alphamax), + (betamin, betamax ), + (gammamin, gammamax)) + # + # TUI ADAO + # -------- + case = adaoBuilder.New() + case.set( 'AlgorithmParameters', + Algorithm = 'NonLinearLeastSquares', # Mots-clé réservé + Parameters = { # Dictionnaire + "Bounds":Bounds, # Liste de paires de Real ou de None + "MaximumNumberOfSteps":100, # Int >= 0 + "CostDecrementTolerance":1.e-7, # Real > 0 + "StoreSupplementaryCalculations":[# Liste de mots-clés réservés + "CostFunctionJAtCurrentOptimum", + "CostFunctionJoAtCurrentOptimum", + "CurrentOptimum", + "SimulatedObservationAtCurrentOptimum", + "SimulatedObservationAtOptimum", + ], + } + ) + case.set( 'Background', + Vector = Xb, # array, list, tuple, matrix + Stored = True, # Bool + ) + case.set( 'Observation', + Vector = observations, # array, list, tuple, matrix + Stored = False, # Bool + ) + case.set( 'ObservationError', + Matrix = None, # None ou matrice carrée + ScalarSparseMatrix = 1.0, # None ou Real > 0 + DiagonalSparseMatrix = None, # None ou vecteur + ) + case.set( 'ObservationOperator', + OneFunction = multisimulation, # MultiFonction [Y] = F([X]) + Parameters = { # Dictionnaire + "DifferentialIncrement":0.0001, # Real > 0 + "CenteredFiniteDifference":False, # Bool + }, + InputFunctionAsMulti = True, # Bool + ) + case.set( 'Observer', + Variable = "CurrentState", # Mot-clé + Template = "ValuePrinter", # Mot-clé + String = None, # None ou code Python + Info = None, # None ou string + + ) + case.execute() + # + # Exploitation independante + # ------------------------- + Xbackground = case.get("Background") + Xoptimum = case.get("Analysis")[-1] + FX_at_optimum = case.get("SimulatedObservationAtOptimum")[-1] + J_values = case.get("CostFunctionJAtCurrentOptimum")[:] + print("") + print("Number of internal iterations...: %i"%len(J_values)) + print("Initial state...................: %s"%(numpy.ravel(Xbackground),)) + print("Optimal state...................: %s"%(numpy.ravel(Xoptimum),)) + print("Simulation at optimal state.....: %s"%(numpy.ravel(FX_at_optimum),)) + print("") + # + # Fin du cas + # ---------- + ecart = assertAlmostEqualArrays(Xoptimum, [ 2., 3., 4.]) + # + print("The maximal absolute error in the test is of %.2e."%ecart) + print("The results are correct.") + print("") + # + return Xoptimum + +# ============================================================================== +def assertAlmostEqualArrays(first, second, places=7, msg=None, delta=None): + "Compare two vectors, like unittest.assertAlmostEqual" + import numpy + if msg is not None: + print(msg) + if delta is not None: + if ( (numpy.asarray(first) - numpy.asarray(second)) > float(delta) ).any(): + raise AssertionError("%s != %s within %s places"%(first,second,delta)) + else: + if ( (numpy.asarray(first) - numpy.asarray(second)) > 10**(-int(places)) ).any(): + raise AssertionError("%s != %s within %i places"%(first,second,places)) + return max(abs(numpy.asarray(first) - numpy.asarray(second))) + +# ============================================================================== +if __name__ == '__main__': + print("\nAUTODIAGNOSTIC\n==============") + unittest.main() diff --git a/PyObjectRAII.hxx b/PyObjectRAII.hxx new file mode 100644 index 0000000..6136719 --- /dev/null +++ b/PyObjectRAII.hxx @@ -0,0 +1,66 @@ +// Copyright (C) 2019 EDF R&D +// +// This library is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 2.1 of the License. +// +// This library is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// Lesser General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License along with this library; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +// +// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com +// +// Author: Anthony Geay, anthony.geay@edf.fr, EDF R&D + +#pragma once + +#include "Python.h" + +class PyObjectRAII +{ +public: + static PyObjectRAII FromBorrowed(PyObject *obj) { IncRef(obj); return PyObjectRAII(obj); } + static PyObjectRAII FromNew(PyObject *obj) { return PyObjectRAII(obj); } + PyObjectRAII():_obj(nullptr) { } + PyObjectRAII(PyObjectRAII&& other):_obj(other._obj) { other._obj=nullptr; } + PyObjectRAII(const PyObjectRAII& other):_obj(other._obj) { incRef(); } + PyObjectRAII& operator=(PyObjectRAII&& other) { unRef(); _obj=other._obj; other._obj=nullptr; } + PyObjectRAII& operator=(const PyObjectRAII& other) { if(_obj==other._obj) return *this; unRef(); _obj=other._obj; incRef(); return *this; } + ~PyObjectRAII() { unRef(); } + PyObject *retn() { incRef(); return _obj; } + operator PyObject *() const { return _obj; } + bool isNull() const { return _obj==nullptr; } +private: + PyObjectRAII(PyObject *obj):_obj(obj) { } + void unRef() { Py_XDECREF(_obj); } + void incRef() { IncRef(_obj); } + static void IncRef(PyObject *obj) { Py_XINCREF(obj); } +private: + PyObject *_obj; +}; + +class AutoGIL +{ +public: + AutoGIL():_gstate(PyGILState_Ensure()) { } + ~AutoGIL() { PyGILState_Release(_gstate); } +private: + PyGILState_STATE _gstate; +}; + +class AutoSaveThread +{ +public: + AutoSaveThread() { unlock(); } + ~AutoSaveThread() { lock(); } + void lock() { PyEval_RestoreThread(_tstate); } + void unlock() { _tstate = PyEval_SaveThread(); } +private: + PyThreadState *_tstate; +}; diff --git a/TestAdaoExchange.cxx b/TestAdaoExchange.cxx new file mode 100644 index 0000000..f038e1d --- /dev/null +++ b/TestAdaoExchange.cxx @@ -0,0 +1,293 @@ +// Copyright (C) 2019 EDF R&D +// +// This library is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 2.1 of the License. +// +// This library is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// Lesser General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License along with this library; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +// +// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com +// +// Author: Anthony Geay, anthony.geay@edf.fr, EDF R&D + +#include "TestAdaoExchange.hxx" + +#include "AdaoExchangeLayer4Quintet.hxx" +#include "AdaoExchangeLayerException.hxx" +#include "AdaoModelKeyVal.hxx" +#include "PyObjectRAII.hxx" + +#include "py2cpp/py2cpp.hxx" + +#include +#include + +std::vector func(const std::vector& vec) +{ + return {vec[0],2.*vec[1],3.*vec[2],vec[0]+2.*vec[1]+3.*vec[2]}; +} + +PyObject *multiFunc(PyObject *inp) +{ + PyGILState_STATE gstate(PyGILState_Ensure()); + PyObjectRAII iterator(PyObjectRAII::FromNew(PyObject_GetIter(inp))); + if(iterator.isNull()) + throw AdaoExchangeLayerException("Input object is not iterable !"); + PyObject *item(nullptr); + PyObjectRAII numpyModule(PyObjectRAII::FromNew(PyImport_ImportModule("numpy"))); + if(numpyModule.isNull()) + throw AdaoExchangeLayerException("Failed to load numpy"); + PyObjectRAII ravelFunc(PyObjectRAII::FromNew(PyObject_GetAttrString(numpyModule,"ravel"))); + std::vector< PyObjectRAII > pyrets; + while( item = PyIter_Next(iterator) ) + { + PyObjectRAII item2(PyObjectRAII::FromNew(item)); + { + PyObjectRAII args(PyObjectRAII::FromNew(PyTuple_New(1))); + PyTuple_SetItem(args,0,item2.retn()); + PyObjectRAII npyArray(PyObjectRAII::FromNew(PyObject_CallObject(ravelFunc,args))); + // Waiting management of npy arrays into py2cpp + PyObjectRAII lolistFunc(PyObjectRAII::FromNew(PyObject_GetAttrString(npyArray,"tolist"))); + PyObjectRAII listPy; + { + PyObjectRAII args2(PyObjectRAII::FromNew(PyTuple_New(0))); + listPy=PyObjectRAII::FromNew(PyObject_CallObject(lolistFunc,args2)); + } + std::vector vect; + { + py2cpp::PyPtr listPy2(listPy.retn()); + py2cpp::fromPyPtr(listPy2,vect); + } + // + PyGILState_Release(gstate); + std::vector res(func(vect)); + gstate=PyGILState_Ensure(); + // + py2cpp::PyPtr resPy(py2cpp::toPyPtr(res)); + PyObjectRAII resPy2(PyObjectRAII::FromBorrowed(resPy.get())); + pyrets.push_back(resPy2); + } + } + std::size_t len(pyrets.size()); + PyObjectRAII ret(PyObjectRAII::FromNew(PyList_New(len))); + for(std::size_t i=0;i vect; + { + py2cpp::PyPtr obj(optimum_4_py2cpp); + py2cpp::fromPyPtr(obj,vect); + } + CPPUNIT_ASSERT_EQUAL(3,(int)vect.size()); + CPPUNIT_ASSERT_DOUBLES_EQUAL(2.,vect[0],1e-7); + CPPUNIT_ASSERT_DOUBLES_EQUAL(3.,vect[1],1e-7); + CPPUNIT_ASSERT_DOUBLES_EQUAL(4.,vect[2],1e-7); +} + +void AdaoExchangeTest::testBlue() +{ + class TestBlueVisitor : public RecursiveVisitor + { + public: + void visit(GenericKeyVal *obj) + { + EnumAlgoKeyVal *objc(dynamic_cast(obj)); + if(objc) + objc->setVal(EnumAlgo::Blue); + } + void enterSubDir(DictKeyVal *subdir) { } + void exitSubDir(DictKeyVal *subdir) { } + }; + + MainModel mm; + // + TestBlueVisitor vis; + mm.visitAll(&vis); + // + AdaoExchangeLayer4Quintet adao; + adao.init(&mm); + { + std::string sciptPyOfModelMaker(mm.pyStr()); + std::cerr << sciptPyOfModelMaker << std::endl; + } + adao.execute(); + PyObject *listOfElts( nullptr ); + while( adao.next(listOfElts) ) + { + PyObject *resultOfChunk(multiFunc(listOfElts)); + adao.setResult(resultOfChunk); + } + PyObject *res(adao.getResult()); + PyObjectRAII optimum(PyObjectRAII::FromNew(res)); + PyObjectRAII optimum_4_py2cpp(NumpyToListWaitingForPy2CppManagement(optimum)); + std::vector vect; + { + py2cpp::PyPtr obj(optimum_4_py2cpp); + py2cpp::fromPyPtr(obj,vect); + } + CPPUNIT_ASSERT_EQUAL(3,(int)vect.size()); + CPPUNIT_ASSERT_DOUBLES_EQUAL(2.,vect[0],1e-7); + CPPUNIT_ASSERT_DOUBLES_EQUAL(3.,vect[1],1e-7); + CPPUNIT_ASSERT_DOUBLES_EQUAL(4.,vect[2],1e-7); +} + +void AdaoExchangeTest::testNonLinearLeastSquares() +{ + class TestNonLinearLeastSquaresVisitor : public RecursiveVisitor + { + public: + void visit(GenericKeyVal *obj) + { + EnumAlgoKeyVal *objc(dynamic_cast(obj)); + if(objc) + objc->setVal(EnumAlgo::NonLinearLeastSquares); + } + void enterSubDir(DictKeyVal *subdir) { } + void exitSubDir(DictKeyVal *subdir) { } + }; + + MainModel mm; + // + TestNonLinearLeastSquaresVisitor vis; + mm.visitAll(&vis); + // + AdaoExchangeLayer4Quintet adao; + adao.init(&mm); + { + std::string sciptPyOfModelMaker(mm.pyStr()); + std::cerr << sciptPyOfModelMaker << std::endl; + } + adao.execute(); + PyObject *listOfElts( nullptr ); + while( adao.next(listOfElts) ) + { + PyObject *resultOfChunk(multiFunc(listOfElts)); + adao.setResult(resultOfChunk); + } + PyObject *res(adao.getResult()); + PyObjectRAII optimum(PyObjectRAII::FromNew(res)); + PyObjectRAII optimum_4_py2cpp(NumpyToListWaitingForPy2CppManagement(optimum)); + std::vector vect; + { + py2cpp::PyPtr obj(optimum_4_py2cpp); + py2cpp::fromPyPtr(obj,vect); + } + CPPUNIT_ASSERT_EQUAL(3,(int)vect.size()); + CPPUNIT_ASSERT_DOUBLES_EQUAL(2.,vect[0],1e-7); + CPPUNIT_ASSERT_DOUBLES_EQUAL(3.,vect[1],1e-7); + CPPUNIT_ASSERT_DOUBLES_EQUAL(4.,vect[2],1e-7); +} + +CPPUNIT_TEST_SUITE_REGISTRATION( AdaoExchangeTest ); + +#include +#include +#include +#include +#include +#include +#include +#include + +int main(int argc, char* argv[]) +{ + // --- Create the event manager and test controller + CPPUNIT_NS::TestResult controller; + + // --- Add a listener that collects test result + CPPUNIT_NS::TestResultCollector result; + controller.addListener( &result ); + + // --- Add a listener that print dots as test run. +#ifdef WIN32 + CPPUNIT_NS::TextTestProgressListener progress; +#else + CPPUNIT_NS::BriefTestProgressListener progress; +#endif + controller.addListener( &progress ); + + // --- Get the top level suite from the registry + + CPPUNIT_NS::Test *suite = + CPPUNIT_NS::TestFactoryRegistry::getRegistry().makeTest(); + + // --- Adds the test to the list of test to run + + CPPUNIT_NS::TestRunner runner; + runner.addTest( suite ); + runner.run( controller); + + // --- Print test in a compiler compatible format. + std::ofstream testFile; + testFile.open("test.log", std::ios::out | std::ios::app); + testFile << "------ ADAO exchange test log:" << std::endl; + CPPUNIT_NS::CompilerOutputter outputter( &result, testFile ); + outputter.write(); + + // --- Run the tests. + + bool wasSucessful = result.wasSuccessful(); + testFile.close(); + + // --- Return error code 1 if the one of test failed. + + return wasSucessful ? 0 : 1; +} diff --git a/TestAdaoExchange.hxx b/TestAdaoExchange.hxx new file mode 100644 index 0000000..fd00d4b --- /dev/null +++ b/TestAdaoExchange.hxx @@ -0,0 +1,40 @@ +// Copyright (C) 2019 EDF R&D +// +// This library is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 2.1 of the License. +// +// This library is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// Lesser General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License along with this library; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +// +// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com +// +// Author: Anthony Geay, anthony.geay@edf.fr, EDF R&D + +#pragma once + +#include + +class AdaoExchangeTest : public CppUnit::TestFixture +{ + CPPUNIT_TEST_SUITE(AdaoExchangeTest); + CPPUNIT_TEST(test3DVar); + CPPUNIT_TEST(testBlue); + CPPUNIT_TEST(testNonLinearLeastSquares); + CPPUNIT_TEST_SUITE_END(); +public: + void setUp(); + void tearDown(); + void cleanUp(); +public: + void test3DVar(); + void testBlue(); + void testNonLinearLeastSquares(); +}; diff --git a/ThreeDVarCase.py b/ThreeDVarCase.py new file mode 100644 index 0000000..1b6b9f3 --- /dev/null +++ b/ThreeDVarCase.py @@ -0,0 +1,96 @@ +# coding: utf-8 +# +# Copyright (C) 2019 EDF R&D +# +# This library is free software; you can redistribute it and/or +# modify it under the terms of the GNU Lesser General Public +# License as published by the Free Software Foundation; either +# version 2.1 of the License. +# +# This library is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +# Lesser General Public License for more details. +# +# You should have received a copy of the GNU Lesser General Public +# License along with this library; if not, write to the Free Software +# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +# +# See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com +# +# Author: Anthony Geay, anthony.geay@edf.fr, EDF R&D + +# ***** +# ***** +# ***** +# ***** + +def BuildCase(cppFunc): + def evaluator( xserie ): + import numpy as np + yserie = [np.array(elt) for elt in cppFunc(xserie)] + return yserie + + from adao import adaoBuilder + # + Xb = (5.0, 7, 9.0) + observations = [2, 6, 12, 20] + alphamin, alphamax = 0., 10. + betamin, betamax = 3, 13 + gammamin, gammamax = 1.5, 15.5 + Bounds = ( + (alphamin, alphamax), + (betamin, betamax ), + (gammamin, gammamax)) + # + # TUI ADAO + # -------- + case = adaoBuilder.New() + case.set( 'AlgorithmParameters', + Algorithm = '3DVAR', # Mots-clé réservé + Parameters = { # Dictionnaire + "Bounds":Bounds, # Liste de paires de Real ou de None + "MaximumNumberOfSteps":100, # Int >= 0 + "CostDecrementTolerance":1.e-7, # Real > 0 + "StoreSupplementaryCalculations":[# Liste de mots-clés réservés + "CostFunctionJAtCurrentOptimum", + "CostFunctionJoAtCurrentOptimum", + "CurrentOptimum", + "SimulatedObservationAtCurrentOptimum", + "SimulatedObservationAtOptimum", + ], + } + ) + case.set( 'Background', + Vector = Xb, # array, list, tuple, matrix + Stored = True, # Bool + ) + case.set( 'Observation', + Vector = observations, # array, list, tuple, matrix + Stored = False, # Bool + ) + case.set( 'BackgroundError', + Matrix = None, # None ou matrice carrée + ScalarSparseMatrix = 1.0e10, # None ou Real > 0 + DiagonalSparseMatrix = None, # None ou vecteur + ) + case.set( 'ObservationError', + Matrix = None, # None ou matrice carrée + ScalarSparseMatrix = 1.0, # None ou Real > 0 + DiagonalSparseMatrix = None, # None ou vecteur + ) + case.set( 'ObservationOperator', + OneFunction = evaluator, # MultiFonction [Y] = F([X]) multisimulation + Parameters = { # Dictionnaire + "DifferentialIncrement":0.0001, # Real > 0 + "CenteredFiniteDifference":False, # Bool + }, + InputFunctionAsMulti = True, # Bool + ) + case.set( 'Observer', + Variable = "CurrentState", # Mot-clé + Template = "ValuePrinter", # Mot-clé + String = None, # None ou code Python + Info = None, # None ou string + ) + return case diff --git a/start.tar.gz b/start.tar.gz new file mode 100644 index 0000000000000000000000000000000000000000..6a81d0a8c00a49938d88e52c5b501c281c7023a3 GIT binary patch literal 13713 zcmZX2Q;aSQu;kdbjWf1wn`dm><{8_zZQHhO+n(?K+1!WCZl}_zbUjqRba&E(QBXjv z(r!9Hpx3wqFA_(?baY?=^CkBVD0P^6e-23jjW@a#q#SRM}&U&f8hq0X8U1Dy!GYRD3KmNHVrJbO5qlZlDXdk7JE zV7e`50&fasu9P^UEP0qUTWmP5K5ebcegZm1u>Y9t#loWVm>0V%Hl#)YFUkv<+)NaN zHoA0~;)u)I<6m_&7|roGgbQOZq49 zJujQ_^+?N6r->}I#W0zpt0-{)aobUp{_L(^{J;g|f%mr3hV5Azlg# ztPUh2=v2YQh_e9DSG}kDdD(iMdc!I)0*3I=ljmGQjyXINt^*D;9U+*FOL`IzVcTLgguPP*Q8T0|Z0P zD9~x!K3Qk4%Ooe1JwvlCf`vGj@amY`MNq(~kJm4p0#`T(&s$xBa?YOu6HQcnqQhEC zm9hv7*~Zq@wed`XNw^mqX>sTJ#p|l4=2qs>jTMPcM_=EaWZzDF4RkG9 z2%{;MP|q$u{{n7U=<_%;*4zxNl8OibUv&G|%p`;TM)OzJ0H2Jp*f0-{_x)-{2 z(%0NRN5f8Ml9f77&KW9{TO$z^tR6Bch)6p z|{|a9lBPI658c-kj&n z>>sQ2_PPTJXf*oBBn@Hm60j}dSFcQT#F1lgn+QcI?upQGdm>x#i(8mS#1({Ndyi9F z4@jn0oIbVsa*1#J&v{5{i)CE=_9s&Y}jClPZc0dUm#8Xm(C5WBOKKs4UlnSXe zvMXihCP{bn@*BgN?Pwy)LT>gP14=#Y=x?%?KI=c5wt~oggU_{MzfX%AV+Nq4Xu0{_ zU|kRuOj2eY({nkXX2DYhEt~AXR5CWF-{S1Wq@3NOh(N;_Akqu91T}NIkEG9|J+_we zJpq`Xja)(yPpdM#-qHt0S@USAo-Mm60!A==N^WrY1)M2MNkig)DX?56cCx<NfzxJuKojiETv4bVG((&OGSk{Cdit%Ih zH`R33RlFg=2qthKyPH;I8XbMOcnE2o%)-G8k^{1-@FD)9!`C2SWdAKTaM{va8Ua*( z7NY)cX@bUzZEQkV3t9QYaRx+Sae(%K`Lo6jne{LD1`Jq&{*=uGh~Th_T8oFKO^-Hn zd58ie`_pX3p;>_lX@F!_;RvellZC7d0KE5Yl+fuTi%wNC75&a4h8*b4Si%53<8>sQ zVkUVCRASp~$ukMc;5E{t{djA{{hQ(W%VXT3-txfuqeKV!cpAc=OX5J8op*F?n<_Be+*N}2&Tp^dpp z*{bk?LMJ$Hb~fzzm|G1vhs~BIA>$*y{kYE2K7WohIyx_5cvHha1~5P|3VBuq*$$1k zOW3=usAjNz{vdOX6h3+Fapu?R9b&982lxoCHxQ+tzhR1NG9w8BpuOvA3TaB*2X1^D zmu!R2!tJO%B%!zL&kz02Uk$g-mm;7MG>j2*Gt|;Fd&h>z1ju}ypnBaU`@2LhV=uV< z8*^Y!SEsQ4LX$lDVCeJ$q+9P&;pPsZ(qj-`(O3y>ezM=XdUVya+yVkyPcoH~okRRe zSo00(WzY09(YVs|s50S!i}cU^ekL4^r$Ee3zsL0tNRIen+Tf4ed1bNCp+AH}Sl z3(R+IBps6|;r7v2_x|{i3$3}a_3!MOjsMcWsUOidn85ee{gJE)@uDO-b z=hh2$Z?E({ctc|ZIEK2eUs`y1{`kfjU}l_+gdqK-F=qM#Ns)qL=JJ3 zY7JBRpapfMj5@l|YUor?S!Y&GVg~ob<+xZ6RFJ~*oFpZ{FC$HCrF7vwyn723hQ?cum1!$`GEA$?NK*_Y* zux3g~Isj46b^Oqcq+hrd2xD0e=WfNHEfrTYk$DRhfjK{Eu2A$&_BfM3R$C(Ck=fNn z6bkBCQ7&gg)o%w;W}T{Y1+A*>xS*^OV2D*OX#T4$gDtMsMdnGgK(GDdWmD(G(28hs zL9e;TM35rZcHFwPAXT^3u{3dm!Ee_vrgz50D{#jfX6oM+opc1udNP-po$xxUNsK;j*8J!8e_%s2`4$8fc5_Uq$!BH)1x4-?7B^pUqLgc$EfNcM--f z8!dz1rwh+sbX!n3q#4;duEk;?2^|aFKWFd)Mi+RB4981v)e7C~4CrqBDFe+Hr!P3~ zaoSvq0>SVS4UFv-Xs5V~*S8b}zn~p`{*zb4G=<9Yt1i^Brn!TUXsH-zSe>aUqzKz? zW0n`5DLOU+#Q49pH_V^slH5JnH{O%NLm9ycs38 zO;^o&9?AH6jYRQtZVTWd(9edBE^ICN-j`{!=fs0zvy)Q8W4L5H-2773k&~?t0+139 zG4O<~jv>4*6_r7ZwI*N5^8rl;v19LgiomRAP;TQA4Ce;BYj~tR>{H`|&aPvmLFh3{ zR{pj9BxcYO>gAZOR7;96%42xdFFVUWArfr$cj3M)uHKj&#*f+}pH>g7iTK!u_ASwU zSW|6J!3=F#eCwz%r*qP?X-Ws{nMuem zZDzLO08?ETZ?W*I#a1i~KEs=l>77N3^6}JUSx-DqSyYDYt@?|oue zUaXWDm(VeOl$~ZH!8+3!3e>GGN1wO4X0w#`o<&MClrU2=$r9iD*oq|5MUCdAbE2@0 z{$W7ShPKE`b<)fMg*2^#BwR>AY@dH7H$S{{$S#yddb)a%JnNDpu7@atoYSPAV zSZ;L{D07yS(yB$E_KhDihA5z(%J@CWMpw%W1((h?pd-apM`O|Ry6_b`!p90ZS51tf zzpSxb%{11Sab!J53HL=j+Uj3G+w>{)2Q<0)=IFgTbn|u2&emMr!^7N}J>M9}+D!kc z@H~v^LT4?{ln<(@Toqn=wg_U8hD6!%CevGt?J$k$KNpSQu0?`Phf_D_l@Omvu1Xo}J5(bc9()u*y zNen2JS?EKbFB($fQWZ&Z=6pCgGz=)o7``4~3?0tiSK?Ki{*C>C>2RH>@x#yE?demR z_Qua$Zad)Td+g`V{^83J@bw=deCz8Q`MD!vWiY3B;ESqN-Rq+-%{W#f>R$)uVuZI?xF>yo-)>oj|$rzF4W_x;5H%fH%qUB#eAsPR{dgjcWT9_=1QAr++w z@sT@UeY{nH-Uo9OX*UO5q2M0xQEQL?a=E z1eNhe(-iy&8Vy>9l*Tb5R1}=n+LRSi0+W(9iJ*%juKw2@>5kG6>V5WIT})$csH-PeU(%CK>;*;idv zuzyKSg)^ZiGtAR?m#C`jy7Ozp1Z|1PT1oA%#T?61qt0A&jpRnv{6822b-OgpA)8~k zOkEpQ595wW@bA=%%!176YGWArrAOW+Gpk$k1hQRU$QmbEOdFEC zT$ql?O_Cnjw^@+fekVe<4f8y;f>-3v2~`wTcnNsH3A}(XdLFLu-2jf;KV}6we^tp@ zT&5dj{C(4?68=yNO8e^uCX0cOi#?SV)^Y|;W@e(V44%G}DCV}HXs}A(yh42S^%CUv z@bNK`@X^)S~+Xftu7usTM>1eFcCjtpj+djk}6 z8eP9_qWq@POre4jYTZ;2tf)vNmh_Qdyd@ag*H1A%&U&?N3519XRUn@ze-YNNU!DB}Uh z_r!@7k-atbN7q?i;`>cIPJTPQXM)#O8n!D3T=AwPe}=N4g9CO4mE38Vv#Z+OrGlw_yBB^2*}eD?~F zeDgW$!>$|k^6u$RNBX_Wl$GbJG>n&z4^Wcg5aZV6r~hCRau&cP-cIm8)(5z!tj3_8 zNue*ReUJK{qy1{dPFIDwtYs0Xo36F9G zIaUUw!opa&2QK2cA{N)lEbN;c9y!gcehjmf5x)X|c->JT86`ovB1|Vl=QQdxPQBa> zMdqOUvRH?=|LmWfJ^9VA{zzV`@&L=(kW<|op{q9%Uzf<>Qad$c;T_a7f`70Wp=fuX zFZZ&qMj)>iw-0E0tDdpln%tfHMKAZ+RIHem^c&`T<#boomfybG&fW~5Xyp8RHsf_CEkS{jG0Dwuk793@z= z>3E%;XW{%LcE2#WOYuRH;A@ zxT7kCvXc;3+j^kqrmDjswmbCUzGNp-8}fMgXJv$niuSbL1uZ6G-%`-*)vA{%Opa2o zrWmMs{A9Q3g)k_#Sp>+iWEvi;yHsd2ajaG;R78wPgXyI5qNH&rIIOH_@`SJAy2OmY z`Bs=$9i?lYRQpkL=mwNyb+Cu{4B9;&1f&f!}; z4Y;b$DRE!3Y{3{i%gnNp`M5HdS!N^| z_JxrFQM_DY8*s|$hLLRWxEJijCAh&Z=|)W6gA_#U^xpyvfV2^8(-YU@FMw_Dwr9mR zz*_AV(5AZ$_$5;~0BF1e(0~2&mH`6399_5jV`=D{X?{TX@A>`*D(f`v*+9I+)Y-pG zZcmw+1@DfY9V_#*Oyde3#k#%De5v1J7E^?aiPZkqkXa{`ye#y4h zXh~01(9%6p`RJK^!I@~&Q>KWsAVV55n3i{zq2`KPE7oIqbmMmzT z5{i^4y!vpSpHxjVC^C0hzStr5dzmVDezgQ35P88`!jExWJm%8bw&dpO!$SQh@@W~D zS}DwaR23o@Jngfx!c6S#Pw-I1V7pP3Skht}Mq}LQle5&tST}u9n-{VTdHtg>%5&^z zm{+cHMJttI#ZxO?H7TKOk5{24u8UZr1Ssxq9;i*9ov}&5r_jc*?^&^c4yc+yc`YgM z2|X_P=CV%A>8i95B#U2{XlbM=ImnPO{4z~b#+JLu)Z-OI4XPdXD6$sM?r23zns@sD4FzJakda`XuPhbt%#fP@n$l(lZu2;Gs!HJVA5m2^YL4xkAH^Dk{ea7|9+mj}aEJ=(7@zJq4vw z>ee=kaFpg!L|DLDMc{zBu?K5vw6whBacq@BFYusKhkq#w^t&A%vc1=hKtIT}Fe^SM3N z32ZdS_2+-_&@hw18p)c9Zo(p*e{^c9e7Zo**srP+^^Bs@;#Vh55twt``XL{aY$Ov` zs7{b)sZMXu96S(i=RhsIW=++}xNOSVB1!WyKG|5-x})q|yR8nW!o8?ZT2Gpxp>WMl zu@L{7VjP$#N6|J%YuZiSP&*{Rzhy77j_nU#WMb&4@Y1)K+72RM=w&ENzt-;_qE`A* zUS;Fe1oAIbMp`x7lC=EE@~E+Cq@U@Sw<+Hl4f%9@z^A^hT%8!tj=GCpByp>*t(K=N zhXk>|RSb=~3|v@VfhQ5c%SqoulPtUsGJbVMmR_;WgZynh0AIr{lrh%x+ghG_kNQ=&J$S zkGE+ncDrA>+{S$eyL#^U`TM~E-=-`tkN=~=NhOzPEoWy?OqH7|GJZM6OjA+aUsaN} zLuh8!lK`rxA#qv}&**Ci;nKN1Ev;t`Ohevq`pHORgflL!)WBUZyS=7i3V@5x=g-sp{yG`vSt5G$wqM&|v_qVX5BxCq;?ICVm| zVlhT+1q4;;5V2w54Sn3W=zY)wZIQ*U zoBB~DR&A&AORK)iw&}}f{hZ||Nb;P1)VaX01(V^X;<{7`Mxufn$I@wXtnR+~f@YLD@=YDQknJbEJ@>kt|4SXz6eHCW@FADqtdIKIw1m)RwdA9@f zVzlspvga#Nvu!_VvP?XA`qTBw~a>IJ)lAzqTaYF#Nff z1J%5`gr_+m5qFbM8t=DI?spaE=+2A&AR(}4AN}eZs$KK&UhKAN^ly)p;NZlnns1JZ z8r?uf(*Wd&Uu``VT`h@6X@;nT~zbT=kU~Ju?`A7?`;qaG&)3U5-2Ci(60!W1C$yPzpzW5 zC)5-qWKIF2yoTuMGwKD3#9QSLzuaX-q!SUu7l0LFndM31$Ic@_^F{uR`Sk~0=r`0(@{^2Ilki@cR)4^7y0Ryl>R9XxB#b{uIxJ6Gs>t&5hC8OtOXn2j2c8aU+3;6NM*%2s+1dJTAMXEiA(hz`_N7rV_= ztb)9=k}qO;Y0sF`U^m4CmAk$hSdOogI3BHOmGEJu-$y5~jd6aBh~Sh-+-3;^O051% zq|T*o;QV5O$pvEecJ#(X8t~k_3c^wxj`jKR2=uUxcByQS@7-XacYx33%GaYpW*|XP zEY(-vjjb!d+SWblyRV9eA?jyOg*OD?egM#Y_Xi~0`md2UA_qyWV@HFAj3!J?Q^^tS zo6e(4FpWI0;p^}eUr^%>O_)y3FigJd>{<&&QW~8}9t9+b0Uky$ZNiV#0~WgWP*y7S zA(=50T3PKvnI@;Hn2*sK4G65q^))#?Uf*hn^B%8Nnq-4`I| z9OG}!3<|X@go)VuGIdMZBqdV$Am;#VvI@3EhtPA;Z=&nr-&iK+Cu`BS&TuZLh+$e z>451Zo6bamBf-cKRAdhiu98%Yi=*s0J#HK9w%-QrxhEE^;0h0CMXcK+X zfG@~+l|(U$$jqNA{Z0z(k7|sa!x-X{K0HDHXj75RF!r6?%ZahU0ACS}wUuhy5#Ud#WiO z@(eyhpKmT#?HYdlBg-@am{$O!9{}1Fz+@y3xuc?7p+qM|pNgTV;9BNP_TopH1)_l3 z(X~x+(p*?|;$Gzv_8CEd-f)!V)O<^?;nu{VjL?!bsW0HbOM@K*YbY%1($;jzKi0lw zYZAkMtGnjDm?7yukP+Of({+7%uwpAvThp{TFdnj*t_pKi)&4A_Bu|{;bb;N%mD%TN z)djzUHF{~%FhCq-z!>P){yXbzjDq)+5w3kI8@&I{Z@*Nvw#0i*;?pbz50>(O2HIx2 z0G}1lR#bh&0FXYWLe&>*6``}b2ge$M6?WJ{MGd-zMHI5M#|rVkMEkiUbrFmi`w41| zRVx=VlWmR?(-u`~z@3f|l|G3I^g&+p7CYo&WUnXHH0_IIKO1CZEC?*|F7CM)L|an7 z)`2OY&lBuP%wp)v339jI_o%CrFI1)}XAXP%BEnXv()^2_YjY8|%OMv>bWP78C{_BX z-?@6wV7SPDSv4jCA$j_8!#ZMW_w%0bJ-a1mtP*5{Gam?p(mekOe8*90FJn zY8(8B3F5NXqn+TNdp&{WO3_u-Q5|oR;#43z5jGj0`%O`pybyBG)t~D#z#6L!)$^!CHyjcD0cNN`%d$qJz0$s*T$K zC5Xy;seetei)ouhrFb z4p43ZS2!t4Esa7g4q#K4PqznO6Z;YhbPa6lqU%Z%$l}XFSaU!tV3A6k_^J1Is6?>t zhsYJb#drJjqSadlT3dAd@E%PQ(7>s>?4^`C@_(E+iAP*2^%kzLuVFG)rXy7o&af_L z%#I!_SKilLf0nl@`}IiZjm3RSB~SRc`+h&07tk29_|c=ks0d2zar#q%%4cBe>44Z# z1Gy<)MgZ^CYdNAen$T$4hlyYpltFBq_>ay_ zM;}3$ko)Vkq8pe4<&g-wFrkWWzE4SDFb;4&ue_Q+=VPyk(ey08V053gI)7|0-FHrp< z8&UR;gYcJp8gE(>TpC^JC}`sWGiwa1Bb&$!#AvYuS-FRcU4jYo-yxc>^V{GJhf^jh zf(R85#lU*>q+D3Jg+tu=tx!aBGg*AF~!uTR))iS+b7AmjO?@M>^cQWV*s0gG$*~8a><~>{D zS8`fir*Ip`)kb7O+)$GylObuAU9>3M%6S zsAU7FXZ?M-0f<|lzjE`q;KsN2@=;EXI6OU9QfT~xb(+=fgNTfv5Y38nx5eEqRHFp& zPo5D!Q@P3KO)OM);IXC>Yn(iT9ePwHCQy>iJ2z8%ljmtZ(r`Z51~U$Ab^b>b5JyQO zJFjm%$hkUwj_VMfoL(N$+@;fiuBBccZq! z-wTZnI-_!<53%r7vyh(ZdYtc#>c`5ho`o+r0WkhU`x`_TLvoYK&*TGqICSYV`*;NG zRwCc1#JjB@euXo8F#RWYqM_J6k-x*H1OfVsnhe^u9fitGu9p=X)FXk+@-XqHYm~)g zC;}YAA_8}eS_^^lVuW-=CuGPm3yw z;#NdD=~0UihHw1 zW>d0Fm^0{hX5O;UF`P`*c_9wvQ

l=w5FVW4zDmhv5rMMIjCCXmAqxD7c$4+KUNf zWG9XSh`Hev8p>Ww$*K)*W~b=XG819un1^=&?3|@Vsc?%``6JM=BLB|!EyWP_g(D?i zPZtmjA*f2#Lpr>bd}0gEawWPpSH^*9*EDHVDmBne!hi{{${Me05lyU9Zu0^&-hcnq zAL=B~7CT8LAPNxXPeifVEJC{kKzCMKo2mF2FbV%9@dCTRu`H73Q}Byd;vzv`5Sao+ z*X`k-n+T>P9Q=15Rndmsssv=%Cg>o+vexcm5EvE)PsVe(C3-nj==^SO5a7OGewguA z^yo5>N0poC6X#YV%EsmUaS9x+P7$w+1~|V*axtt1UeV9*o@ds$m;FTnuZMQm_(XH_ zu@b;Q*t0-O&dtqDPFlR)R<#`fbIVhw`TdDw4RHS}FxXx~JXwLE)#6QLubl&=@YZ7S zw=dg>Q6HeOwYyVW^6Bs1RuABt0sR_0!!1QX>eGA)SgB)&F*PWFqfHKEH{ItD@~(G$ z;J*$bHHi3p@8!EtT?ZVlCA?p@U092AX>#WR9l4-jM=#rwVmh+Hzak|D0+xum_Dp;3i;^HZgU+pkr+! zBlHh#3_=vFiRXE?mX5HbjaQzqQuC4wgxR7`BMfvti3xK-QAEHb^gj@mT=G47%ftll zsV#jf70J`9P?hebWpqMaXE^gm(d-xcP^s2Vd7K!l=yYy^OLr0O$1f#Ei{xZ5`1gQy z>{NW>_*OkVvc!a+Y`P!@Re7~gL2c(mY2UuSsrzF%If6due#Jt zu}Ey_Zv_-nQ|W)%aAEVv*y>gVb~3VZv!=s+Mk>nk>fx`Zu|j{8#m&Qq1&SH*%};VH z)Z((fW>VbT;HHm{dcne(iI!wloY$er^f;O|G+JC4H3m zc7+k#zN#{_s+DoUJFX*BBg8xT>5 z>0=XOGI9$4+w)I@_@@)OHc5-&EJ{gTO-Z{M$+RZVX>YgqLF!?2Z&Qlpj#mu}5@&Vr zrTiB^)dyBTFUpY{$cvd{&stco&A0jDtx|gDA4~a=oq9P!_Mv+~tG{)XDiK&CkXweY zj+4mIF}&HJHMQ#P-}2#qx4+waaNmn(|0GHsyJ>H0(V3QWcHzf7#FFw(h-ira3#AV<>B`$67kVFz;!z=4y^B!p`)FM zJwfC*T{c3kRnFyk<4R)c>3ecDi$HERI_M^WAB_YzFcVJ36(oa9wrF7$R<8sVZZxUG z^3`kCm}i3-)&F=a_VscLa98y4^VsX<;Qe|VV30x>4>lZ+@V!J7c3`0%)%NLw*qR{!eY4`Bk8blz9^r_DlbN z6?>s(*e(sjH|lG>ke#Y}Z*;f1zUm!oSM|EVQ?>POn5?zE;Zrr<%M2^P(7t7y4PdL* zHtV%d&IOC^|8+}p|Nq@|^tQGExzoQI=~FPy8C~st%^lSaX*AQJAsh9KxU)by>F}^k khKB64)1j$}F#pT`(>x-}{x6)x0lozadQ4w|EP#Rj2mbTa5C8xG literal 0 HcmV?d00001 -- 2.39.2