From: Anthony Geay Date: Wed, 27 Feb 2019 08:34:35 +0000 (+0100) Subject: First implementation X-Git-Tag: V9_4_0~14 X-Git-Url: http://git.salome-platform.org/gitweb/?a=commitdiff_plain;h=3a322a64957a19c0ed44ebf86545dea5fff5230f;p=tools%2Fadao_interface.git First implementation --- 3a322a64957a19c0ed44ebf86545dea5fff5230f 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 0000000..6a81d0a Binary files /dev/null and b/start.tar.gz differ