Homec4science
Diffusion LibMultiScale (features/components) Browse

libmultiscale/src/filter/compute_python.cc
b05d5f2aba78features/components

Last Change
Raw File

compute_python.cc
View Options

/**
* @file compute_python.cc
*
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
*
* @date Tue Jul 22 14:47:56 2014
*
* @brief This compute allows to use Python to produce computed values as input to LM
*
* @section LICENSE
*
* Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* LibMultiScale 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 3 of the License, or (at your option) any
* later version.
*
* LibMultiScale 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 LibMultiScale. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "lm_common.hh"
#include "lib_md.hh"
#include "compute_python.hh"
#include "filter_manager.hh"
#include "communicator.hh"
#include "units.hh"
#include "lm_parser.hh"
/* -------------------------------------------------------------------------- */
#define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION
#include <numpy/arrayobject.h>
/* -------------------------------------------------------------------------- */
__BEGIN_LIBMULTISCALE__
/* -------------------------------------------------------------------------- */
// PYTHON units accessors for the compute scripts
/* -------------------------------------------------------------------------- */
#include <Python.h>
static PyObject * python_lm_exit(PyObject *self, PyObject *args) {
UInt code = PyInt_AsLong(args);
DUMP("Terminating LibMultiScale simulation with error code " << code,DBG_MESSAGE);
switch (code){
case LM_EXIT_FAILURE:
MPI_Abort(MPI_COMM_WORLD, EXIT_FAILURE);
lm_exit(code);break;
case LM_EXIT_SUCCESS:
MPI_Finalize();
lm_exit(code);break;
default:
MPI_Abort(MPI_COMM_WORLD, code);
lm_exit(code);
}
}
#define DECLARE_UNITS_ACCESSOR(__NAME__) \
static PyObject * python_##__NAME__(PyObject *self, PyObject *args) { \
return PyFloat_FromDouble(code_unit_system.__NAME__); \
};
DECLARE_UNITS_ACCESSOR(mvv2e )
DECLARE_UNITS_ACCESSOR(f_m2v_t )
DECLARE_UNITS_ACCESSOR(ft_m2v )
DECLARE_UNITS_ACCESSOR(mv_t2f )
DECLARE_UNITS_ACCESSOR(kT2e )
DECLARE_UNITS_ACCESSOR(kT2fd )
DECLARE_UNITS_ACCESSOR(m_tt2f_d )
DECLARE_UNITS_ACCESSOR(e_dd2m_tt)
DECLARE_UNITS_ACCESSOR(e2fd )
DECLARE_UNITS_ACCESSOR(e_m2dd_tt)
DECLARE_UNITS_ACCESSOR(fd2e )
#undef DECLARE_UNITS_ACCESSOR
#define DECLARE_PYTHON_UNITS_FUNCTION(__NAME__) \
{#__NAME__, python_##__NAME__, METH_NOARGS, #__NAME__}
/* -------------------------------------------------------------------------- */
static PyMethodDef pyLMUnitsMethods[] = {
DECLARE_PYTHON_UNITS_FUNCTION(mvv2e ),
DECLARE_PYTHON_UNITS_FUNCTION(f_m2v_t ),
DECLARE_PYTHON_UNITS_FUNCTION(ft_m2v ),
DECLARE_PYTHON_UNITS_FUNCTION(mv_t2f ),
DECLARE_PYTHON_UNITS_FUNCTION(kT2e ),
DECLARE_PYTHON_UNITS_FUNCTION(kT2fd ),
DECLARE_PYTHON_UNITS_FUNCTION(m_tt2f_d ),
DECLARE_PYTHON_UNITS_FUNCTION(e_dd2m_tt),
DECLARE_PYTHON_UNITS_FUNCTION(e2fd ),
DECLARE_PYTHON_UNITS_FUNCTION(e_m2dd_tt),
DECLARE_PYTHON_UNITS_FUNCTION(fd2e ),
{"lm_exit", python_lm_exit, METH_O, "Clean termination of a libmultiscale simulation"},
{NULL, NULL, 0, NULL} /* Sentinel */
};
/* -------------------------------------------------------------------------- */
void appendLMConfEnvironment(PyObject * kwargs){
std::map<std::string, double> & vars = Parser::getAlgebraicVariables();
std::map<std::string, double>::iterator it = vars.begin();
std::map<std::string, double>::iterator end = vars.end();
while (it != end) {
std::string varname = it->first;
Real varval = it->second;
PyObject * pyvalue = PyFloat_FromDouble(varval);
PyDict_SetItemString(kwargs,varname.c_str(),pyvalue);
++it;
}
}
/* -------------------------------------------------------------------------- */
template <typename T>
void appendGeneralVariables(PyObject * kwargs,
const std::string & name, T & value){
LM_TOIMPLEMENT;
}
/* -------------------------------------------------------------------------- */
template <>
void appendGeneralVariables<std::string>(PyObject * kwargs,
const std::string & name,
std::string & value){
PyObject * pyvalue = PyString_FromString(value.c_str());
PyDict_SetItemString(kwargs,name.c_str(),pyvalue);
}
/* -------------------------------------------------------------------------- */
template <>
void appendGeneralVariables<const std::string>(PyObject * kwargs,
const std::string & name,
const std::string & value){
PyObject * pyvalue = PyString_FromString(value.c_str());
PyDict_SetItemString(kwargs,name.c_str(),pyvalue);
}
/* -------------------------------------------------------------------------- */
template <>
void appendGeneralVariables<IntegrationSchemeStage>(PyObject * kwargs,
const std::string & name,
IntegrationSchemeStage & value){
std::stringstream sstr;
sstr << value;
std::string str = sstr.str();
appendGeneralVariables<std::string>(kwargs,name,str);
}
/* -------------------------------------------------------------------------- */
template <>
void appendGeneralVariables<Real>(PyObject * kwargs,
const std::string & name, Real & value){
PyObject * pyvalue = PyFloat_FromDouble(value);
PyDict_SetItemString(kwargs,name.c_str(),pyvalue);
}
/* -------------------------------------------------------------------------- */
template <>
void appendGeneralVariables<UInt>(PyObject * kwargs,
const std::string & name, UInt & value){
PyObject * pyvalue = PyInt_FromLong(value);
PyDict_SetItemString(kwargs,name.c_str(),pyvalue);
}
/* -------------------------------------------------------------------------- */
template <typename Cont>
void appendGeneralArray(PyObject * kwargs,
Cont & cont,
bool gather_flag,
CommGroup group,
bool gather_root_proc){
const UInt Dim = cont.getDim();
std::vector<Real> * vdata_ptr = NULL;
if (gather_flag){
vdata_ptr = &cont.gatherAllData(gather_root_proc);
}
else {
std::vector<Real> & tmp = dynamic_cast<ComputeInterface & >(cont);
vdata_ptr = &tmp;
}
if (gather_flag){
Communicator & comm = Communicator::getCommunicator();
UInt my_rank = comm.groupRank(lm_my_proc_id, group);
if (my_rank != gather_root_proc) {
// this->setDim(0);
return;
}
}
LM_ASSERT(vdata_ptr,"null pointer detected: abort");
std::vector<Real> & vdata = *vdata_ptr;
UInt sz = vdata.size()/Dim;
npy_intp dims[2] = {sz,Dim};
PyObject * pyvalue = PyArray_SimpleNewFromData(2, dims,NPY_DOUBLE, &vdata[0]);
DUMP("adding the compute " << cont.getID(), DBG_INFO);
PyDict_SetItemString(kwargs,cont.getID().c_str(),pyvalue);
PyObject* pyNameComputed = PyDict_GetItemString(kwargs, "name_computed");
if (!pyNameComputed) pyNameComputed = PyDict_New();
PyObject* pyThisNameComputed = PyDict_GetItemString(pyNameComputed,
cont.getID().c_str());
if (pyThisNameComputed) LM_FATAL("the compute " << cont.getID()
<< " was already inputed: useless");
pyThisNameComputed = PyDict_New();
for (UInt i = 0; i < cont.name_computed.size(); ++i) {
PyObject * index = PyInt_FromLong(i);
PyDict_SetItemString(pyThisNameComputed,
cont.name_computed[i].c_str(),index);
}
PyDict_SetItemString(pyNameComputed,
cont.getID().c_str(),pyThisNameComputed);
PyDict_SetItemString(kwargs,
"name_computed",pyNameComputed);
}
/* -------------------------------------------------------------------------- */
void freeArgumentsDict(PyObject * kwargs){
PyObject *key, *value;
Py_ssize_t pos = 0;
std::vector<PyObject *> objects_to_destroy;
while (PyDict_Next(kwargs, &pos, &key, &value)) {
// objects_to_destroy.push_back(key);
objects_to_destroy.push_back(value);
}
PyDict_Clear(kwargs);
for (UInt i = 0; i < objects_to_destroy.size(); ++i) {
if (PyDict_Check(objects_to_destroy[i])){
if (PyString_Check(objects_to_destroy[i])){
char * tmp = PyString_AsString(objects_to_destroy[i]);
DUMP("wants to free " << tmp,DBG_MESSAGE);
}
}
if (PyDict_Check(objects_to_destroy[i]))
freeArgumentsDict(objects_to_destroy[i]);
Py_XDECREF(objects_to_destroy[i]);
}
}
/* -------------------------------------------------------------------------- */
template <typename _Input>
void ComputePython<_Input>::init(){
if (!Py_IsInitialized()){
DUMP("Initialize python",DBG_INFO_STARTUP);
Py_Initialize();
DUMP("import file into a module",DBG_INFO_STARTUP);
PyRun_SimpleString("import sys");
PyRun_SimpleString("sys.path.append(\"./\")");
// code_unit_system.print();
pUnitsModule = Py_InitModule("LMUnitsModule",pyLMUnitsMethods);
}
pModule = PyImport_ImportModule(filename.c_str());
if (pModule == NULL) {
DUMP(std::endl << "Python error in " << filename << std::endl << std::endl,DBG_MESSAGE);
PyErr_Print();
LM_FATAL("Failed to load " << filename);
}
pCompute = PyObject_GetAttrString(pModule, "compute");
if (!pCompute){
DUMP(std::endl << "Python error in" << filename << std::endl << std::endl,DBG_MESSAGE);
PyErr_Print();
LM_FATAL("Cannot find entry point: compute()");
}
if (!PyCallable_Check(pCompute)) LM_FATAL("Cannot find entry point: compute()");
pSetNameComputed = NULL;
if (PyObject_HasAttrString(pModule, "setNameComputed")){
pSetNameComputed = PyObject_GetAttrString(pModule, "setNameComputed");
if (pSetNameComputed){
if (!PyCallable_Check(pSetNameComputed))
LM_FATAL("Cannot find name function used to set column names");
}
}
}
/* -------------------------------------------------------------------------- */
template <typename Cont>
void ComputePython<Cont>::build(Cont & cont){
this->empty();
CommGroup group = this->getCommGroup();
PyObject * pArgs = PyTuple_New(0);
import_array();
PyObject * kwargs = NULL;
Communicator & comm = Communicator::getCommunicator();
UInt my_rank = comm.groupRank(lm_my_proc_id, group);
if (!this->gather_flag || (my_rank == gather_root_proc))
kwargs = PyDict_New();
appendGeneralArray(kwargs, cont, this->gather_flag,group,this->gather_root_proc);
if (cont.size()){
DUMP("adding the compute " << cont.getID() << " for " << this->getID(), DBG_INFO);
} else {
DUMP("not adding the compute " << cont.getID() << " for " << this->getID(), DBG_INFO);
DUMP("THIS REQUIRES FIXING BY GUILLAUME", DBG_MESSAGE);
}
for (UInt comp = 0; comp < compute_list.size(); ++comp) {
FilterID id = compute_list[comp];
ComputeInterface * my_compute
= dynamic_cast<ComputeInterface*>(FilterManager::getManager().getObject(id));
if (!my_compute)
LM_FATAL("invalid compute passed to dumper text " << compute_list[comp]);
my_compute->build();
appendGeneralArray(kwargs, *my_compute,
this->gather_flag,group,this->gather_root_proc);
if (my_compute->size()){
DUMP("adding the compute " << my_compute->getID() << " for " << this->getID(), DBG_INFO);
} else {
DUMP("not adding the compute " << my_compute->getID() << " for " << this->getID(), DBG_INFO);
DUMP("THIS REQUIRES FIXING BY GUILLAUME", DBG_MESSAGE);
}
}
if (this->gather_flag && (my_rank != gather_root_proc)) return;
appendLMConfEnvironment(kwargs);
appendGeneralVariables(kwargs,"step",current_step);
appendGeneralVariables(kwargs,"current_stage",current_stage);
appendGeneralVariables(kwargs,"compute_name",this->getID());
appendGeneralVariables(kwargs,"rank",lm_my_proc_id);
PyObject * pValue = PyObject_Call(pCompute, pArgs,(PyObject*)kwargs);
if (!pValue) {
DUMP(std::endl << "Python error in " << filename
<< std::endl << std::endl,DBG_MESSAGE);
PyErr_Print();
LM_FATAL("the value returned by the script is invalid");
}
if (!PyArray_Check(pValue)) {
DUMP("the object returned by the script is not a numpy array but is:"
<< std::endl,DBG_MESSAGE);
PyObject_Print(pValue,stderr,Py_PRINT_RAW);
std::cerr << std::endl << std::endl;
LM_FATAL("abort");
}
PyArray_Descr * numpy_type = (PyArray_Descr*)PyArray_DESCR((PyArrayObject*)pValue);
if (numpy_type->kind != 'f') LM_FATAL("You must return a float D-type kind for your numpy in your script");
UInt _n = PyArray_NDIM((PyArrayObject*)pValue);
if (_n == 1) this->setDim(1);
else if (_n == 2) {
npy_intp * ndims = PyArray_DIMS((PyArrayObject*)pValue);
this->setDim(ndims[1]);
}
else LM_FATAL("output vector should not have more than 2 dimensions: " << _n);
// DUMP("type name:" << pValue->ob_type->tp_name,DBG_MESSAGE);
PyArrayIterObject *iter = (PyArrayIterObject *)PyArray_IterNew(pValue);
if (iter == NULL) {
DUMP(std::endl << "Python error in " << filename << std::endl << std::endl,DBG_MESSAGE);
PyErr_Print();
LM_FATAL("the return value was not an iterable array");
}
setNameComputed(pArgs,kwargs);
while (iter->index < iter->size) {
// DUMP(*(Real *)iter->dataptr,DBG_MESSAGE);
this->add(*(Real *)iter->dataptr);
PyArray_ITER_NEXT(iter);
}
DUMP("compute python " << this->getID() << " computed "
<< this->size() << " values",DBG_INFO);
Py_XDECREF(iter);
Py_XDECREF(pArgs);
Py_XDECREF(pValue);
// PyDict_Clear(kwargs);
freeArgumentsDict(kwargs);
Py_XDECREF(kwargs);
}
/* -------------------------------------------------------------------------- */
template <typename Cont>
void ComputePython<Cont>::setNameComputed(PyObject* pArgs, PyObject* kwargs){
this->name_computed.clear();
if (pSetNameComputed){
PyObject * pValue = PyObject_Call(pSetNameComputed, pArgs,kwargs);
if (!pValue) {
DUMP(std::endl << "Python error in " << filename
<< std::endl << std::endl,DBG_MESSAGE);
PyErr_Print();
LM_FATAL("the value returned by the script is invalid");
}
if (!PyTuple_Check(pValue)){
DUMP("the object returned by the script is not a tuple array but is:"
<< std::endl,DBG_MESSAGE);
PyObject_Print(pValue,stderr,Py_PRINT_RAW);
std::cerr << std::endl << std::endl;
LM_FATAL("abort");
}
Py_ssize_t n_names = PyTuple_Size(pValue);
if (n_names != this->getDim())
LM_FATAL("The number of names provided is not correct "
<< n_names << " != " << this->getDim());
for (UInt i = 0; i < n_names; ++i) {
PyObject * entry = PyTuple_GetItem(pValue, i);
if (!PyString_Check(entry)) LM_FATAL("entry " << i
<< " of returned tuple is not a string");
char * name = PyString_AsString(entry);
DUMP("adding name_computed " << name,DBG_INFO);
this->name_computed.push_back(std::string(name));
}
}
else{
for (UInt i = 0; i < this->getDim(); ++i){
std::stringstream title;
title << "python" << i;
this->name_computed.push_back(title.str());
}
}
}
/* -------------------------------------------------------------------------- */
/* LMDESC PYTHON
This computes an array of reals based on a python script
*/
/* LMEXAMPLE
COMPUTE disp EXTRACT INPUT md FIELD displacement \\\\
COMPUTE disp2 PYTHON INPUT disp FILENAME script \\\\
with a python script stored in a file script.py with the content such as:\\\\
import numpy as np\\ \\
def compute(**kwargs):\\
~~vec = np.copy(kwargs["disp"])\\
~~vec *= 2. \\
~~return vec
*/
template <typename _Input>
void ComputePython<_Input>::declareParams(){
/* LMKEYWORD FILENAME
Specify the name of the python script to use
*/
this->parseKeyword("FILENAME",filename);
/* LMKEYWORD GATHER_ROOT
Specify the processor where to gather data before inputing to python script
*/
this->parseKeyword("GATHER_ROOT",gather_root_proc,0u);
/* LMKEYWORD GATHER
Specify the processor where to gather data before inputing to python script
*/
this->parseTag("GATHER",gather_flag,false);
/* LMKEYWORD ADD_COMPUTE
Add computes to input to the kwargs of the python script
*/
this->parseKeyword("ADD_COMPUTE",compute_list);
}
/* -------------------------------------------------------------------------- */
DECLARE_COMPUTE_REAL(ComputePython);
__END_LIBMULTISCALE__
c4science · Help