Changeset abad620 in sasview for pr_inversion

Timestamp:

May 8, 2008 10:55:47 AM (16 years ago)

Author:

Mathieu Doucet <doucetm@…>

Branches:

master, ESS_GUI, ESS_GUI_Docs, ESS_GUI_batch_fitting, ESS_GUI_bumps_abstraction, ESS_GUI_iss1116, ESS_GUI_iss879, ESS_GUI_iss959, ESS_GUI_opencl, ESS_GUI_ordering, ESS_GUI_sync_sascalc, costrafo411, magnetic_scatt, release-4.1.1, release-4.1.2, release-4.2.2, release_4.0.1, ticket-1009, ticket-1094-headless, ticket-1242-2d-resolution, ticket-1243, ticket-1249, ticket885, unittest-saveload

Children:

4241e48

Parents:

2e07e8f

Message:

Added alpha estimate

Location:

pr_inversion

Files:

• c_extensions/Cinvertor.c (modified) (6 diffs)
• c_extensions/invertor.c (modified) (2 diffs)
• c_extensions/invertor.h (modified) (1 diff)
• invertor.py (modified) (2 diffs)

pr_inversion/c_extensions/Cinvertor.c

@@ -319 +319 @@
         double regterm = 0.0;
         double tmp = 0.0;
+        // Number of slices in regularization term estimate
+        int nslice = 25;
         
         PyObject *data_obj;
@@ -331 +333 @@
     residuals = PyList_New(self->params.npoints);
 
-    regterm = reg_term(pars, self->params.d_max, npars);
+    regterm = reg_term(pars, self->params.d_max, npars, nslice);
     
     for(i=0; i<self->params.npoints; i++) {
@@ -367 +369 @@
         double regterm = 0.0;
         double tmp = 0.0;
+        // Number of slices in regularization term estimate
+        int nslice = 25;
         
         PyObject *data_obj;
@@ -378 +382 @@
     residuals = PyList_New(self->params.npoints);
 
-    regterm = reg_term(pars, self->params.d_max, npars);
+    regterm = reg_term(pars, self->params.d_max, npars, nslice);
 
     
@@ -465 +469 @@
         if (!PyArg_ParseTuple(args, "did", &d_max, &n, &q)) return NULL;
         return Py_BuildValue("f", ortho_transformed(d_max, n, q));      
+        
+}
+
+static PyObject * oscillations(Cinvertor *self, PyObject *args) {
+        double *pars;
+        PyObject *data_obj;
+        Py_ssize_t npars;
+        double oscill, norm;
+        
+        if (!PyArg_ParseTuple(args, "O", &data_obj)) return NULL;
+        OUTVECTOR(data_obj,pars,npars);
+        
+        oscill = reg_term(pars, self->params.d_max, npars, 100);
+        norm   = int_p2(pars, self->params.d_max, npars, 100);
+        return Py_BuildValue("f", sqrt(oscill/norm)/acos(-1.0)*self->params.d_max );    
         
 }
@@ -489 +508 @@
                    {"is_valid", (PyCFunction)is_valid, METH_VARARGS, ""},
                    {"basefunc_ft", (PyCFunction)basefunc_ft, METH_VARARGS, ""},
+                   {"oscillations", (PyCFunction)oscillations, METH_VARARGS, ""},
    
    {NULL}

pr_inversion/c_extensions/invertor.c

@@ -1 +1 @@
 #include <math.h>
 #include "invertor.h"
+#include <memory.h>
+#include <stdio.h>
+#include <stdlib.h>
 
 double pi = 3.1416;
@@ -120 +123 @@
  * regularization term calculated from the expansion.
  */
-double reg_term(double *pars, double d_max, int n_c) {
+double reg_term(double *pars, double d_max, int n_c, int nslice) {
     double sum = 0.0; 
     double r;
     double deriv;
         int i;
-    for (i=0; i<25; i++) {
-        r = d_max/25.0*i;
+    for (i=0; i<nslice; i++) {
+        r = d_max/(1.0*nslice)*i;
         deriv = dprdr(pars, d_max, n_c, r);
         sum += deriv*deriv;
     }
-    return sum/25.0*d_max;
+    return sum/(1.0*nslice)*d_max;
 }
 
+/**
+ * regularization term calculated from the expansion.
+ */
+double int_p2(double *pars, double d_max, int n_c, int nslice) {
+    double sum = 0.0; 
+    double r; 
+    double value;
+        int i;
+    for (i=0; i<nslice; i++) {
+        r = d_max/(1.0*nslice)*i;
+        value = pr(pars, d_max, n_c, r);
+        sum += value*value;
+    }
+    return sum/(1.0*nslice)*d_max;
+}
+

pr_inversion/c_extensions/invertor.h

@@ -26 +26 @@
 double pr(double *c, double d_max, int n_c, double r);
 double dprdr(double *pars, double d_max, int n_c, double r);
-double reg_term(double *pars, double d_max, int n_c);
+double reg_term(double *pars, double d_max, int n_c, int nslice);
+double int_p2(double *pars, double d_max, int n_c, int nslice);
 void pr_err(double *pars, double *err, double d_max, int n_c, 
                 double r, double *pr_value, double *pr_value_err);

pr_inversion/invertor.py

@@ -8 +8 @@
     ## Time elapsed for last computation
     elapsed = 0
+    ## Alpha to get the reg term the same size as the signal
+    suggested_alpha = 0
     
     def __init__(self):
@@ -191 +193 @@
                     
         new_alpha = sum_sig/(sum_reg/self.alpha)
-        print "Suggested alpha =", 0.1*new_alpha
+        self.suggested_alpha = new_alpha
         
         try: