source: sasview/sansmodels/src/sans/models/c_extensions/disperser.c @ 87f51eb

/** 
 * Disperser
 * 
 * Python class that takes a model and averages its output
 * for a distribution of its parameters.
 * 
 * The parameters to be varied are specified at instantiation time.
 * The distributions are Gaussian, with std deviations specified for
 * each parameter at instantiation time.
 * 
 * @author   M.Doucet / UTK 
 */
 
#include <Python.h>
#include "structmember.h"
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>

/// Error object for raised exceptions
static PyObject * DisperserError = NULL;

// Class definition
typedef struct {
    PyObject_HEAD
    /// Dictionary of parameters for this class
    PyObject * params;
    /// Model object to disperse
    PyObject * model;
    /// Array of original values
    double * centers;
    /// Number of parameters to disperse
    int n_pars;
    /// Coordinate system flag
    int cylindrical;
} Disperser;

// Function definitions
double disperseParam(Disperser *self, int iPar, PyObject * pars);
double Disperser_readDouble(PyObject *p);

/**
 * Function to deallocate memory for the object
 */
static void
Disperser_dealloc(Disperser * self) {
    free(self->centers);
    self->ob_type->tp_free((PyObject*)self);
}

/**
 * Function to create a new Disperser object
 */
static PyObject *
Disperser_new(PyTypeObject *type, PyObject *args, PyObject *kwds)
{
    Disperser *self;
    
    self = (Disperser *)type->tp_alloc(type, 0);
   
    return (PyObject *)self;
}

/**
 * Initialization of the object (__init__)
 */
static int
Disperser_init(Disperser *self, PyObject *args, PyObject *kwds)
{
        PyObject * model;
        PyObject * paramList;
        PyObject * sigmaList;
        
    if (self != NULL) {
        
        // Create parameters
        self->params = PyDict_New();
        
                // Get input
                if ( !PyArg_ParseTuple(args,"OOO",&model, &paramList, &sigmaList) ) {
                    PyErr_SetString(DisperserError, 
                        "Disperser: could not parse to initialize. Need Disperser(model, paramList, sigmaList)");
                        return -1;
                }

                // Check that the lists are of the right type
                if( PyList_Check(paramList)!=1 || PyList_Check(sigmaList)!=1 ) {
                        PyErr_SetString(DisperserError,
                                "Disperser: bad input parameters; did not find lists.");
                        return -1;
                }
                
                // Check that the list lengths are compatible
                if( PyList_GET_SIZE(paramList) != PyList_GET_SIZE(sigmaList) ) {
                        PyErr_SetString(DisperserError,
                                "Disperser: supplied lists are not of the same length.");
                        return -1;
                } 
                
                // Set the data members
                self->model = model;
                
                // Number of parameters to disperse
                // Let run() find out since values can change
                self->n_pars = 0;

                // Set default coordinate system as cartesian
                self->cylindrical = 0;

        // Initialize parameter dictionary
        // Store the list of parameter names (strings)
        PyDict_SetItemString(self->params,"paramList", paramList);
        
        // Store the list of std deviations (floats)
        PyDict_SetItemString(self->params,"sigmaList", sigmaList);
        
        // Set the default number of points to average over
        PyDict_SetItemString(self->params,"n_pts",Py_BuildValue("i",25));
    }
    return 0;
}

/**
 * List of member to make public
 * Only the parameters and the model are available to python users
 */
static PyMemberDef Disperser_members[] = {
    {"params", T_OBJECT, offsetof(Disperser, params), 0,
     "Parameters"},
    {"model", T_OBJECT, offsetof(Disperser, model), 0,
     "Model to disperse"},
    {NULL}  /* Sentinel */
};

/** Read double from PyObject
 *      TODO: THIS SHOULD REALLY BE ELSEWHERE!
 *   @param p PyObject
 *   @return double
 */
double Disperser_readDouble(PyObject *p) {
        double value = 0.0;
        
    if (PyFloat_Check(p)==1) {
        value = (double)(((PyFloatObject *)(p))->ob_fval);
        //Py_DECREF(p);
        return value;
    } else if (PyInt_Check(p)==1) {
        value = (double)(((PyIntObject *)(p))->ob_ival);
        //Py_DECREF(p);
        return value;
    } else if (PyLong_Check(p)==1) {
        value = (double)PyLong_AsLong(p);
        //Py_DECREF(p);
        return value;
    } else {
        //Py_DECREF(p);         
        return 0.0;
    }
}


/**
 * Function to call to evaluate model
 * @param args: input q or [qx,qy]
 * @return: function value
 */
static PyObject * runXY(Disperser *self, PyObject *args) {
        PyObject* pars;
        int n_pts;
        int i;
        double value;
        PyObject * result;
        PyObject * temp;
        
        // Get parameters
        
        // Check that the lists are of the right type
        if( PyList_Check(PyDict_GetItemString(self->params, "paramList"))!=1 
                || PyList_Check(PyDict_GetItemString(self->params, "sigmaList"))!=1 ) {
                PyErr_SetString(DisperserError,
                        "Disperser: bad input parameters; did not find lists.");
                return NULL;
        }

        // Reader parameter dictionary
    n_pts = PyInt_AsLong( PyDict_GetItemString(self->params, "n_pts") );
    
        // Check that the list lengths are compatible
        if( PyList_GET_SIZE(PyDict_GetItemString(self->params, "paramList")) 
                != PyList_GET_SIZE(PyDict_GetItemString(self->params, "sigmaList")) ) {
                PyErr_SetString(DisperserError,
                        "Disperser: supplied lists are not of the same length.");
                return NULL;
        } 

    // Number of parameters to disperse
    self->n_pars = PyList_GET_SIZE(PyDict_GetItemString(self->params, "paramList"));
    
    // Allocate memory for centers
    free(self->centers);
    self->centers = (double *)malloc(self->n_pars * sizeof(double));
        if(self->centers==NULL) {
            PyErr_SetString(DisperserError, 
                "Disperser.run could not allocate memory.");
                return NULL;
        }
    
    // Store initial values of those parameters    
    for( i=0; i<self->n_pars; i++ ) {
        result = PyObject_CallMethod(self->model, "getParam", "(O)", 
                PyList_GetItem( PyDict_GetItemString(self->params, "paramList"), i ));
                
        if( result == NULL ) {
                    PyErr_SetString(DisperserError, 
                        "Disperser.run could not get current parameter values.");
                return NULL;
        } else {
                self->centers[i] = Disperser_readDouble(result);
                Py_DECREF(result);
        }; 
    }
        
        // Get input and determine whether we have to supply a 1D or 2D return value.
        if ( !PyArg_ParseTuple(args,"O",&pars) ) {
            PyErr_SetString(DisperserError, 
                "Disperser.run expects a q value.");
                return NULL;
        }    
        
    // Evaluate or calculate average
    if( self->n_pars > 0 ) {
                value = disperseParam(self, 0, pars);
    } else {
         if (self->cylindrical==1) {
                result = PyObject_CallMethod(self->model, "run", "(O)",pars);
         } else {
                result = PyObject_CallMethod(self->model, "runXY", "(O)",pars);
         }
         
         if( result == NULL ) {
                    PyErr_SetString(DisperserError, 
                        "Disperser.run failed.");
                return NULL;            
         } else {
                Py_DECREF(result);
                return result;
         }
    }
    
    // Put back the original model parameter values
    for( i=0; i<self->n_pars; i++ ) {
        // Maybe we need to decref the double here?
        temp = Py_BuildValue("d", self->centers[i]);
        result = PyObject_CallMethod(self->model, "setParam", "(OO)", 
                PyList_GetItem( PyDict_GetItemString(self->params, "paramList"), i ),
                temp );
        Py_DECREF(temp);
                
        if( result == NULL ) {
                    PyErr_SetString(DisperserError, 
                        "Disperser.run could not set back original parameter values.");
                return NULL;
        } else {
                Py_DECREF(result);
        }; 
    }
    
        return Py_BuildValue("d",value);        

}

/**
 * Function to call to evaluate model
 * @param args: input q or [q,phi]
 * @return: function value
 */
static PyObject * run(Disperser *self, PyObject *args) {
        PyObject * output;
        self->cylindrical = 1;
        output = runXY(self, args);
        self->cylindrical = 0;
        return output;
}




/**
 * Gaussian weight
 * @param mean: mean value of the Gaussian
 * @param sigma: standard deviation of the Gaussian
 * @param x: value at which the Gaussian is evaluated
 * @return: value of the Gaussian
 */
double gaussian_weight(double mean, double sigma, double x) {
        double vary, expo_value;
    vary = x-mean;
    expo_value = -vary*vary/(2*sigma*sigma);
    //return 1.0;
    return exp(expo_value);
}

/**
 * @param iPar: ID of parameter to disperse
 * @param pars: input parameter to the evaluation function
 */
double disperseParam(Disperser *self, int iPar, PyObject * pars) {
        double sigma;
        double min_value, max_value;
        double step;
        double prev_value;
        double value_sum;
        double gauss_sum;
        double gauss_value;
        double func_value;
        double error_sys;
        double value;
        int i, n_pts;
        PyObject *result;
        PyObject *temp;
        
        
    n_pts = PyInt_AsLong( PyDict_GetItemString(self->params, "n_pts") );
    
        // If we exhausted the parameter array, simply evaluate
    // the model
    if( iPar < self->n_pars ) {
            
                // Average over Gaussian distribution (2 sigmas)
        value_sum = 0.0;
        gauss_sum = 0.0;
            
        // Get the standard deviation for this parameter
        sigma = Disperser_readDouble(PyList_GetItem( PyDict_GetItemString(self->params, "sigmaList"), iPar ));
            
        // Average over 4 sigmas wide
        min_value = self->centers[iPar] - 2*sigma;
        max_value = self->centers[iPar] + 2*sigma;
            
        // Calculate step size
        step = (max_value - min_value)/(n_pts-1);
            
        // If we are not changing the parameter, just return the 
        // value of the function
        if (step == 0.0) {
            return disperseParam(self, iPar+1, pars);
        }
            
        // Compute average
        prev_value = 0.0;
        error_sys  = 0.0;
        for( i=0; i<n_pts; i++ ) {
            // Set the parameter value           
            value = min_value + (double)i*step;
                
            temp = Py_BuildValue("d", value);
                result = PyObject_CallMethod(self->model, "setParam", "(OO)",
                        PyList_GetItem( PyDict_GetItemString(self->params, "paramList"), iPar ), 
                                temp);
                Py_DECREF(temp);
                        
                    if( result == NULL ) {
                                printf("Could not set param %i\n", iPar);
                        // return a value that will create an NAN
                        return 0.0/(step-step);
                    } else {
                        Py_DECREF(result);
                    }
                
                        gauss_value = gaussian_weight(self->centers[iPar], sigma, value);
            func_value = disperseParam(self, iPar+1, pars);

            value_sum += gauss_value * func_value;
                gauss_sum += gauss_value;       
        }   
        return value_sum/gauss_sum;
        
    } else {
         if (self->cylindrical==1) {
                result = PyObject_CallMethod(self->model, "run", "(O)",pars);
         } else {
                result = PyObject_CallMethod(self->model, "runXY", "(O)",pars);
         }
         if( result == NULL ) {
                printf("Model.run() could not return\n");
                // return a value that will create an NAN
                return 0.0/(step-step);
         } else {
                value = Disperser_readDouble(result);                   
                Py_DECREF(result);
                return value;
         }
    }
        
}

static PyMethodDef Disperser_methods[] = {
    {"run",      (PyCFunction)run     , METH_VARARGS,
      "Evaluate the model at a given input value"},
    {"runXY",      (PyCFunction)runXY     , METH_VARARGS,
      "Evaluate the model at a given input value"},
   {NULL}
};

static PyTypeObject DisperserType = {
    PyObject_HEAD_INIT(NULL)
    0,                         /*ob_size*/
    "Disperser",             /*tp_name*/
    sizeof(Disperser),             /*tp_basicsize*/
    0,                         /*tp_itemsize*/
    (destructor)Disperser_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 */
    0,                         /*tp_call*/
    0,                         /*tp_str*/
    0,                         /*tp_getattro*/
    0,                         /*tp_setattro*/
    0,                         /*tp_as_buffer*/
    Py_TPFLAGS_DEFAULT | Py_TPFLAGS_BASETYPE, /*tp_flags*/
    "Disperser objects",           /* tp_doc */
    0,                         /* tp_traverse */
    0,                         /* tp_clear */
    0,                         /* tp_richcompare */
    0,                         /* tp_weaklistoffset */
    0,                         /* tp_iter */
    0,                         /* tp_iternext */
    Disperser_methods,             /* tp_methods */
    Disperser_members,             /* tp_members */
    0,                         /* tp_getset */
    0,                         /* tp_base */
    0,                         /* tp_dict */
    0,                         /* tp_descr_get */
    0,                         /* tp_descr_set */
    0,                         /* tp_dictoffset */
    (initproc)Disperser_init,      /* tp_init */
    0,                         /* tp_alloc */
    Disperser_new,                 /* tp_new */
};


static PyMethodDef module_methods[] = {
    {NULL} 
};

/**
 * Function used to add the model class to a module
 * @param module: module to add the class to
 */ 
void addDisperser(PyObject *module) {
        PyObject *d;
        
    if (PyType_Ready(&DisperserType) < 0)
        return;

    Py_INCREF(&DisperserType);
    PyModule_AddObject(module, "Disperser", (PyObject *)&DisperserType);
    
    d = PyModule_GetDict(module);
    DisperserError = PyErr_NewException("Disperser.error", NULL, NULL);
    PyDict_SetItemString(d, "DisperserError", DisperserError);
}