Changeset 5da3cc5 in sasview for sansmodels

Timestamp:

Jul 10, 2012 7:59:33 AM (12 years ago)

Author:

Kieran Campbell <kieranrcampbell@…>

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:

83d9120

Parents:

c2e5898

Message:

Implemented erf(x) in libfunc.c and added pass down to C++ of multfactor

Location:

sansmodels

Files:

• include/dabmodel.h (modified) (1 diff)
• include/libmultifunc/libfunc.h (modified) (1 diff)
• src/c_models/libfunc.c (modified) (2 diffs)
• src/python_wrapper/WrapperGenerator.py (modified) (4 diffs)
• src/python_wrapper/classTemplate.txt (modified) (1 diff)
• src/python_wrapper/modelTemplate.txt (modified) (3 diffs)

sansmodels/include/dabmodel.h

@@ -4 +4 @@
 #include "parameters.hh"
 
+/**
+        This software was developed by Institut Laue-Langevin as part of
+        Distributed Data Analysis of Neutron Scattering Experiments (DANSE).
+
+        Copyright 2012 Institut Laue-Langevin
+
+**/
 
 

sansmodels/include/libmultifunc/libfunc.h

@@ -10 +10 @@
 double gamln(double x);
 
+void gser(float *gamser, float a, float x, float *gln);
+
+void gcf(float *gammcf, float a, float x, float *gln);
+
+float gammp(float a,float x);
+
+float erff(float x);
+
 #endif

sansmodels/src/c_models/libfunc.c

@@ -11 +11 @@
 //used in Si func
 
-int factorial(int i)
+int factorial(int i) {
+
+        int k, j;
+        if (i<2){
+                return 1;
+        }
+
+        k=1;
+
+        for(j=1;j<i;j++) {
+                k=k*(j+1);
+        }
+
+        return k;
+
+}
+
+
+
+// Used in pearl nec model
+
+// Sine integral function: approximated within 1%!!!
+
+// integral of sin(x)/x up to namx term nmax=6 looks the best.
+
+double Si(double x)
 
 {
-
-        int k, j;
-
-        if (i<2){
-
-                return 1;
-
-        }
-
-        k=1;
-
-        for(j=1;j<i;j++)
-
-        {
-
-                k=k*(j+1);
-
-        }
-
-        return k;
-
-}
-
-
-
-// Used in pearl nec model
-
-// Sine integral function: approximated within 1%!!!
-
-// integral of sin(x)/x up to namx term nmax=6 looks the best.
-
-double Si(double x)
-
-{
-
         int i;
-
         int nmax=6;
-
         double out;
-
         long double power;
-
         double pi = 4.0*atan(1.0);
 
         if (x >= pi*6.2/4.0){
-
-
-
                 double out_sin = 0.0;
-
                 double out_cos = 0.0;
-
                 out = pi/2.0;
 
-                for (i=0; i<nmax-2; i+=1){
-
+                for (i=0; i<nmax-2; i+=1) {
                         out_cos += pow(-1.0, i) * (double)factorial(2*i) / pow(x, 2*i+1);
-
                         out_sin += pow(-1.0, i) * (double)factorial(2*i+1) / pow(x, 2*i+2);
-
                 }
 
                 out -= cos(x) * out_cos;
-
                 out -= sin(x) * out_sin;
-
                 return out;
-
         }
 
         out = 0.0;
 
-        for (i=0; i<nmax; i+=1)
-
-        {
-
-                if (i==0){
-
+        for (i=0; i<nmax; i+=1) {
+                if (i==0) {
                         out += x;
-
                         continue;
-
                 }
 
                 power = pow(x,(2 * i + 1));
-
                 out += (double)pow(-1, i) * power / ((2.0 * (double)i + 1.0) * (double)factorial(2 * i + 1));
 
                 //printf ("Si=%g %g %d\n", x, out, i);
-
         }
 
         return out;
-
 }
 
@@ -138 +104 @@
 }
 
+/** Modifications below by kieranrcampbell@gmail.com
+    Institut Laue-Langevin, July 2012
+**/
+
+/**
+   Implements eq 6.2.5 (small gamma) of Numerical Recipes in C, essentially
+   the incomplete gamma function multiplied by the gamma function.
+   Required for implementation of fast error function (erf)
+**/
+
+
+#define ITMAX 100
+#define EPS 3.0e-7
+#define FPMIN 1.0e-30
+
+void gser(float *gamser, float a, float x, float *gln) {
+  int n;
+  float sum,del,ap;
+
+  *gln = gamln(a);
+  if(x <= 0.0) {
+    if (x < 0.0) printf("Error: x less than 0 in routine gser");
+    *gamser = 0.0;
+    return;
+  } else {
+    ap = a;
+    del = sum = 1.0/a;
+    
+    for(n=1;n<=ITMAX;n++) {
+      ++ap;
+      del *= x/ap;
+      sum += del;
+      if(fabs(del) < fabs(sum)*EPS) {
+        *gamser = sum * exp(-x + a * log(x) - (*gln));
+        return;
+      }
+    }
+    printf("a too large, ITMAX too small in routine gser");
+    return;
+
+  }
+
+
+}
+
+/**
+   Implements the incomplete gamma function Q(a,x) evaluated by its continued fraction 
+   representation 
+**/
+
+void gcf(float *gammcf, float a, float x, float *gln) {
+  int i;
+  float an,b,c,d,del,h;
+
+  *gln = gamln(a);
+  b = x+1.0-a;
+  c = 1.0/FPMIN;
+  d = 1.0/b;
+  h=d;
+  for (i=1;i <= ITMAX; i++) {
+    an = -i*(i-a);
+    b += 2.0;
+    d = an*d + b;
+    if (fabs(d) < FPMIN) d = FPMIN;
+    c = b+an/c;
+    if (fabs(c) < FPMIN) c = FPMIN;
+    d = 1.0/d;
+    del = d*c;
+    h += del;
+    if (fabs(del-1.0) < EPS) break;
+  }
+  if (i > ITMAX) printf("a too large, ITMAX too small in gcf");
+  *gammcf = exp(-x+a*log(x)-(*gln))*h;
+  return;
+}
+
+
+/**
+   Represents incomplete error function, P(a,x)
+**/
+float gammp(float a, float x) {
+  float gamser,gammcf,gln;
+  if(x < 0.0 || a <= 0.0) printf("Invalid arguments in routine gammp");
+  if (x < (a+1.0)) {
+    gser(&gamser,a,x,&gln);
+    return gamser;
+  } else {
+    gcf(&gammcf,a,x,&gln);
+    return 1.0 - gammcf;
+  }
+}
+
+/** 
+    Implementation of the error function, erf(x)
+**/
+
+float erff(float x) {
+  return x < 0.0 ? -gammp(0.5,x*x) : gammp(0.5,x*x);
+}
+

sansmodels/src/python_wrapper/WrapperGenerator.py

@@ -92 +92 @@
         # Model category
         self.category = None
+        # Whether model belongs to multifunc
+        self.is_multifunc = False
         ## output directory for wrappers
         self.output_dir = output_dir
@@ -228 +230 @@
                     raise ValueError, "Could not parse file %s" % self.file
 
-
+            # is_multifunc
+            key = "[MULTIPLICITY_INFO]"
+            if line.count(key) > 0:
+                self.is_multifunc = True
+                try:
+                    index = line.index(key)
+                    toks = line[index:].split("=")
+                    self.multiplicity_info = toks[1].lstrip().rstrip()
+                except:
+                    raise ValueError, "Could not parse file %s" % self.file
 
             # Catch struct name
@@ -331 +342 @@
             newline = self.replaceToken(newline, 
                                         "[MODELSTRUCT]", self.structName)
+
+            # Sort model initialization based on multifunc
+            if(self.is_multifunc):
+                line = "int level = 1;\nPyArg_ParseTuple(args,\"i\",&level);\n"
+                line += "self->model = new " + self.pythonClass + "(level);"
+            else:
+                line = "self->model = new " + self.pythonClass + "();"
+    
+            newline = self.replaceToken(newline,"[INITIALIZE_MODEL]",
+                                            line)
             
             # Dictionary initialization
@@ -465 +486 @@
             newline = self.replaceToken(newline, 
                                         "[PAR_DETAILS]", self.details)
+            
+            # Call base constructor
+            if self.is_multifunc:
+                newline = self.replaceToken(newline,"[CALL_CPYTHON_INIT]",
+                                            'C' + self.pythonClass + ".__init__(self,multfactor)\n\tself.is_multifunc = True")
+                newline = self.replaceToken(newline,"[MULTIPLICITY_INFO]",self.multiplicity_info)
+            else:
+                newline = self.replaceToken(newline,"[CALL_CPYTHON_INIT]",
+                                            'C' + self.pythonClass + ".__init__(self)\n\tself.is_multifunc = False")
+                newline = self.replaceToken(newline,"[MULTIPLICITY_INFO]","None")
+
            
             # fixed list  details

sansmodels/src/python_wrapper/classTemplate.txt

@@ -86 +86 @@
         self->params = PyDict_New();
         self->dispersion = PyDict_New();
-        self->model = new [CMODEL]();
-        
+
+        [INITIALIZE_MODEL]
+
         [INITDICTIONARY]
          

sansmodels/src/python_wrapper/modelTemplate.txt

@@ -43 +43 @@
     """
         
-    def __init__(self):
+    def __init__(self,multfactor=1):
         """ Initialization """
         self.__dict__ = {}
@@ -50 +50 @@
         BaseComponent.__init__(self)
         #apply([CPYTHONCLASS].__init__, (self,)) 
-        [CPYTHONCLASS].__init__(self)
-        
+
+        [CALL_CPYTHON_INIT]
+                        
         ## Name of the model
         self.name = "[PYTHONCLASS]"
@@ -70 +71 @@
 
         self.category = [CATEGORY]
+        self.multiplicity_info = [MULTIPLICITY_INFO]
+        
 
     def __setstate__(self, state):