Re: r2390 - in /branches/c_modules/maths_fns: c_chi2.c exponential.c relax_fit.c relax_fit.h -- March 16, 2006 - 06:12

Thank you Gary for the help with the C code.  Your fixes have been
merged into the main relax development line, they are much
appreciated.  Apart from the memory leak, which I think is because I
need to call the Numeric-Python C function Py_DECREF() more often, the
code works beautifully.

Thanks,

Edward


On 3/16/06, c.a.macraild@xxxxxxxxxxx <c.a.macraild@xxxxxxxxxxx> wrote:
> Author: macraild
> Date: Wed Mar 15 16:17:06 2006
> New Revision: 2390
>
> URL: http://svn.gna.org/viewcvs/relax?rev=2390&view=rev
> Log:
> A fix for the -fPIC compilation bug (#5473). Some problems remain.
>
> Modified:
>     branches/c_modules/maths_fns/c_chi2.c
>     branches/c_modules/maths_fns/exponential.c
>     branches/c_modules/maths_fns/relax_fit.c
>     branches/c_modules/maths_fns/relax_fit.h
>
> Modified: branches/c_modules/maths_fns/c_chi2.c
> URL: 
> http://svn.gna.org/viewcvs/relax/branches/c_modules/maths_fns/c_chi2.c?rev=2390&r1=2389&r2=2390&view=diff
> ==============================================================================
> --- branches/c_modules/maths_fns/c_chi2.c (original)
> +++ branches/c_modules/maths_fns/c_chi2.c Wed Mar 15 16:17:06 2006
> @@ -20,15 +20,6 @@
>
>  #include <stdio.h>
>  #include <math.h>
> -
> -/* This include must come first */
> -#include <Python.h>
> -
> -/* The header for all functions which will be called */
> -#include "relax_fit.h"
> -
> -/* The Numeric array object header file, must come second */
> -#include <Numeric/arrayobject.h>
>
>  #define square(x) (x)*(x)
>
> @@ -53,15 +44,15 @@
>         */
>
>      /* Declarations */
> -    extern int *num_times;
> -    extern PyArrayObject *values, *sd;
> +    extern int num_times;
> +    extern double *values, *sd;
>      extern double back_calc[];
>         int i;
>         double chi2 = 0.0;
>
>      /* Loop over the time points and sum the chi-squared components */
> -       for (i = 0; i < *num_times; ++i) {
> -               chi2 = chi2 + square((*(double *)values->data[i] - 
> back_calc[i]) / *(double *)sd->data[i]);
> +       for (i = 0; i < num_times; ++i) {
> +               chi2 = chi2 + square((values[i] - back_calc[i]) / sd[i]);
>         }
>
>         return chi2;
>
> Modified: branches/c_modules/maths_fns/exponential.c
> URL: 
> http://svn.gna.org/viewcvs/relax/branches/c_modules/maths_fns/exponential.c?rev=2390&r1=2389&r2=2390&view=diff
> ==============================================================================
> --- branches/c_modules/maths_fns/exponential.c (original)
> +++ branches/c_modules/maths_fns/exponential.c Wed Mar 15 16:17:06 2006
> @@ -22,11 +22,11 @@
>  /* This include must come first */
>  #include <Python.h>
>
> +/* The Numeric array object header file, must come second */
> +#include <Numeric/arrayobject.h>
> +
>  /* The header for all functions which will be called */
>  #include "relax_fit.h"
> -
> -/* The Numeric array object header file, must come second */
> -#include <Numeric/arrayobject.h>
>
>
>
> @@ -35,22 +35,23 @@
>         */
>
>      /* Declarations */
> -    extern int *num_times;
> -    extern PyArrayObject *params, *relax_times;
> +    extern int num_times;
> +    extern double *params, *relax_times;
>      extern double back_calc[];
>      double Rx, I0;
>      int i;
>
>
>      /* Loop over the time points */
> -    for (i = 0; i < *num_times; i++) {
> +    /* for (i = 0; i < num_times; i++) { */
> +    for (i = 0; i < num_times; i++) {
>          /* Zero Rx value */
>          if (params[0] == 0.0)
>              back_calc[i] = 0.0;
>
>          /* Back calculate */
>          else
> -            back_calc[i] = (double *)params->data[1] * exp(-(double 
> *)relax_times->data[i] * (double *)params->data[0]);
> +            back_calc[i] = params[1] * exp(-relax_times[i] * params[0]);
>
>      }
>  }
>
> Modified: branches/c_modules/maths_fns/relax_fit.c
> URL: 
> http://svn.gna.org/viewcvs/relax/branches/c_modules/maths_fns/relax_fit.c?rev=2390&r1=2389&r2=2390&view=diff
> ==============================================================================
> --- branches/c_modules/maths_fns/relax_fit.c (original)
> +++ branches/c_modules/maths_fns/relax_fit.c Wed Mar 15 16:17:06 2006
> @@ -21,58 +21,46 @@
>  /* This include must come first */
>  #include <Python.h>
>
> +/* The Numeric array object header file, must come second */
> +#include <Numeric/arrayobject.h>
> +
>  /* The header for all functions which will be called */
>  #include "relax_fit.h"
> -
> -/* The Numeric array object header file, must come second */
> -#include <Numeric/arrayobject.h>
>
>
>
>  static PyObject *
>  setup(PyObject *self, PyObject *args, PyObject *keywords) {
> -    /* Debugging */
> -    printf("Setting up\n");
> -
>      /* Python declarations */
>      PyObject *values_arg, *sd_arg, *relax_times_arg, *scaling_matrix_arg;
> -    extern PyArrayObject *values, *sd, *relax_times, *scaling_matrix;
> +    PyArrayObject *numpy_values, *numpy_sd, *numpy_relax_times, 
> *numpy_scaling_matrix;
>
>      /* Normal declarations */
> -    extern double *relax_time_array;
> -    extern int *num_params, *num_times;
> +    extern double *values, *sd, *relax_times, *scaling_matrix;
> +    extern double relax_time_array;
> +    extern int num_params, num_times;
>
>      /* The keyword list */
>      static char *keyword_list[] = {"num_params", "num_times", "values", 
> "sd", "relax_times", "scaling_matrix", NULL};
>
> -
>      /* Parse the function arguments */
> -    if (!PyArg_ParseTupleAndKeywords(args, keywords, "iiOOOO", 
> keyword_list, &num_params, &num_times, &values_arg, &sd, &relax_times_arg, 
> &scaling_matrix_arg))
> +    if (!PyArg_ParseTupleAndKeywords(args, keywords, "iiOOOO", 
> keyword_list, &num_params, &num_times, &values_arg, &sd_arg, 
> &relax_times_arg, &scaling_matrix_arg))
>          return NULL;
>
> -    /* Debugging */
> -    printf("Function args have been parsed\n");
> +    /* Make the Numeric arrays contiguous */
> +    numpy_values = (PyArrayObject *) 
> PyArray_ContiguousFromObject(values_arg, PyArray_DOUBLE, 1, 1);
> +    numpy_sd = (PyArrayObject *) PyArray_ContiguousFromObject(sd_arg, 
> PyArray_DOUBLE, 1, 1);
> +    numpy_relax_times = (PyArrayObject *) 
> PyArray_ContiguousFromObject(relax_times_arg, PyArray_DOUBLE, 1, 1);
> +    numpy_scaling_matrix = (PyArrayObject *) 
> PyArray_ContiguousFromObject(scaling_matrix_arg, PyArray_DOUBLE, 2, 2);
>
> -    /* Make the Numeric arrays contiguous */
> -    values = (PyArrayObject *) PyArray_ContiguousFromObject(values_arg, 
> PyArray_DOUBLE, 1, 1);
> -    sd = (PyArrayObject *) PyArray_ContiguousFromObject(sd_arg, 
> PyArray_DOUBLE, 1, 1);
> -    relax_times = (PyArrayObject *) 
> PyArray_ContiguousFromObject(relax_times_arg, PyArray_DOUBLE, 1, 1);
> -    scaling_matrix = (PyArrayObject *) 
> PyArray_ContiguousFromObject(scaling_matrix_arg, PyArray_DOUBLE, 2, 2);
> -
> -    /* Debugging */
> -    printf("Py_INCREF\n");
> -
> -    /* ??? */
> -    Py_INCREF(Py_None);
> -    Py_INCREF(values);
> -    Py_INCREF(sd);
> -    Py_INCREF(relax_times);
> -    Py_INCREF(scaling_matrix);
> -
> -    /* Debugging */
> -    printf("Finished setup\n");
> +    /* Pointers to the Numeric arrays */
> +    values = (double *) numpy_values->data;
> +    sd = (double *) numpy_sd->data;
> +    relax_times = (double *) numpy_relax_times->data;
> +    scaling_matrix = (double *) numpy_scaling_matrix->data;
>
>      /* Return nothing */
> +    Py_INCREF(Py_None);
>      return Py_None;
>  }
>
> @@ -85,13 +73,9 @@
>       * calculated
>       */
>
> -    /* Debugging */
> -    printf("Inside func\n");
> -
>      /* Declarations */
>      PyObject *arg1;
> -    extern PyArrayObject *params;
> -    extern double *params;
> +    extern PyArrayObject *numpy_params;
>      double chi2(void);
>
>
> @@ -100,13 +84,13 @@
>          return NULL;
>
>      /* Convert the Numeric array to be contiguous */
> -    params = (PyArrayObject *) PyArray_ContiguousFromObject(arg1, 
> PyArray_DOUBLE, 1, 1);
> +    numpy_params = (PyArrayObject *) PyArray_ContiguousFromObject(arg1, 
> PyArray_DOUBLE, 1, 1);
> +
> +    /* Pointers to the Numeric arrays */
> +    params = (double *) numpy_params->data;
>
>      /* Back calculated the peak intensities */
>      exponential();
> -
> -    /* ??? */
> -    Py_INCREF(params);
>
>      /* Calculate and return the chi-squared value */
>      return Py_BuildValue("f", chi2());
> @@ -119,8 +103,7 @@
>
>      /* Declarations */
>      PyObject *arg1;
> -    extern PyArrayObject *params;
> -    extern double *params;
> +    extern PyArrayObject *numpy_params;
>
>      /* Temp Declarations */
>      PyArrayObject *aaa_numpy;
> @@ -133,7 +116,10 @@
>          return NULL;
>
>      /* Convert the Numeric array to be contiguous */
> -    params = (PyArrayObject *) PyArray_ContiguousFromObject(arg1, 
> PyArray_DOUBLE, 1, 1);
> +    numpy_params = (PyArrayObject *) PyArray_ContiguousFromObject(arg1, 
> PyArray_DOUBLE, 1, 1);
> +
> +    /* Pointers to the Numeric arrays */
> +    params = (double *) numpy_params->data;
>
>      /* Back calculated the peak intensities */
>      exponential();
> @@ -141,7 +127,8 @@
>
>      /* Test code (convert aaa to a Numeric array */
>      /* aaa_numpy = (PyArrayObject *) PyArray_FromDimsAndData(1, 
> num_params, PyArray_DOUBLE, aaa_pointer); */
> -    aaa_numpy = (PyArrayObject *) PyArray_FromDims(1, (int *) *num_params, 
> PyArray_DOUBLE);
> +    aaa_numpy = (PyArrayObject *) PyArray_FromDims(1, &num_params, 
> PyArray_DOUBLE);
> +    aaa_pointer = (double *) aaa_numpy->data;
>
>      /* Fill the Numeric array */
>      for (i = 0; i < 2; i++)
> @@ -164,14 +151,14 @@
>
>      /* Declarations */
>      extern double back_calc[];
> -    extern int *num_times;
> +    extern int num_times;
>      int i;
>
> -    PyObject *back_calc_py = PyList_New(*num_times);
> +    PyObject *back_calc_py = PyList_New(num_times);
>      assert(PyList_Check(back_calc_py));
>
>      /* Copy the values out of the C array into the Python array */
> -    for (i = 0; i < *num_times; i++)
> +    for (i = 0; i < num_times; i++)
>          PyList_SetItem(back_calc_py, i, Py_BuildValue("f", back_calc[i]));
>
>      /* Return the Numeric array */
>
> Modified: branches/c_modules/maths_fns/relax_fit.h
> URL: 
> http://svn.gna.org/viewcvs/relax/branches/c_modules/maths_fns/relax_fit.h?rev=2390&r1=2389&r2=2390&view=diff
> ==============================================================================
> --- branches/c_modules/maths_fns/relax_fit.h (original)
> +++ branches/c_modules/maths_fns/relax_fit.h Wed Mar 15 16:17:06 2006
> @@ -36,11 +36,16 @@
>  /****************************************/
>
>  /* Variables sent to the setup function to be stored for later use */
> -PyArrayObject *values, *sd, *relax_times, *scaling_matrix;
> -int *num_params, *num_times;
> +PyArrayObject *numpy_values, *numpy_sd, *numpy_relax_times, 
> *numpy_scaling_matrix;
> +int num_params, num_times;
> +double *sd;
>
>  /* Variables sent to 'func', 'dfunc', and 'd2func' during optimisation */
> -PyArrayObject *params;
> +PyArrayObject *numpy_params;
> +
> +/* Pointers to contiguous PyArrayObjects */
> +double *values, *sd, *relax_times, *scaling_matrix;
> +double *params;
>
>  /* Variables used for storage during the function calls of optimisation */
>  double back_calc[MAXTIMES];
>
>
> _______________________________________________
> Relax-commits mailing list
> Relax-commits@xxxxxxx
>
> To unsubscribe from this list, get a password reminder, or change your 
> subscription options, visit the list information page at 
> https://mail.gna.org/listinfo/relax-commits