Changeset 5da3cc5 in sasview for sansmodels


Ignore:
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:
6 edited

Legend:

Unmodified
Added
Removed
  • sansmodels/include/dabmodel.h

    rfa6db8b r5da3cc5  
    44#include "parameters.hh" 
    55 
     6/** 
     7        This software was developed by Institut Laue-Langevin as part of 
     8        Distributed Data Analysis of Neutron Scattering Experiments (DANSE). 
     9 
     10        Copyright 2012 Institut Laue-Langevin 
     11 
     12**/ 
    613 
    714 
  • sansmodels/include/libmultifunc/libfunc.h

    r08648c0 r5da3cc5  
    1010double gamln(double x); 
    1111 
     12void gser(float *gamser, float a, float x, float *gln); 
     13 
     14void gcf(float *gammcf, float a, float x, float *gln); 
     15 
     16float gammp(float a,float x); 
     17 
     18float erff(float x); 
     19 
    1220#endif 
  • sansmodels/src/c_models/libfunc.c

    rf8644dd r5da3cc5  
    1111//used in Si func 
    1212 
    13 int factorial(int i) 
     13int factorial(int i) { 
     14 
     15        int k, j; 
     16        if (i<2){ 
     17                return 1; 
     18        } 
     19 
     20        k=1; 
     21 
     22        for(j=1;j<i;j++) { 
     23                k=k*(j+1); 
     24        } 
     25 
     26        return k; 
     27 
     28} 
     29 
     30 
     31 
     32// Used in pearl nec model 
     33 
     34// Sine integral function: approximated within 1%!!! 
     35 
     36// integral of sin(x)/x up to namx term nmax=6 looks the best. 
     37 
     38double Si(double x) 
    1439 
    1540{ 
    16  
    17         int k, j; 
    18  
    19         if (i<2){ 
    20  
    21                 return 1; 
    22  
    23         } 
    24  
    25         k=1; 
    26  
    27         for(j=1;j<i;j++) 
    28  
    29         { 
    30  
    31                 k=k*(j+1); 
    32  
    33         } 
    34  
    35         return k; 
    36  
    37 } 
    38  
    39  
    40  
    41 // Used in pearl nec model 
    42  
    43 // Sine integral function: approximated within 1%!!! 
    44  
    45 // integral of sin(x)/x up to namx term nmax=6 looks the best. 
    46  
    47 double Si(double x) 
    48  
    49 { 
    50  
    5141        int i; 
    52  
    5342        int nmax=6; 
    54  
    5543        double out; 
    56  
    5744        long double power; 
    58  
    5945        double pi = 4.0*atan(1.0); 
    6046 
    6147        if (x >= pi*6.2/4.0){ 
    62  
    63  
    64  
    6548                double out_sin = 0.0; 
    66  
    6749                double out_cos = 0.0; 
    68  
    6950                out = pi/2.0; 
    7051 
    71                 for (i=0; i<nmax-2; i+=1){ 
    72  
     52                for (i=0; i<nmax-2; i+=1) { 
    7353                        out_cos += pow(-1.0, i) * (double)factorial(2*i) / pow(x, 2*i+1); 
    74  
    7554                        out_sin += pow(-1.0, i) * (double)factorial(2*i+1) / pow(x, 2*i+2); 
    76  
    7755                } 
    7856 
    7957                out -= cos(x) * out_cos; 
    80  
    8158                out -= sin(x) * out_sin; 
    82  
    8359                return out; 
    84  
    8560        } 
    8661 
    8762        out = 0.0; 
    8863 
    89         for (i=0; i<nmax; i+=1) 
    90  
    91         { 
    92  
    93                 if (i==0){ 
    94  
     64        for (i=0; i<nmax; i+=1) { 
     65                if (i==0) { 
    9566                        out += x; 
    96  
    9767                        continue; 
    98  
    9968                } 
    10069 
    10170                power = pow(x,(2 * i + 1)); 
    102  
    10371                out += (double)pow(-1, i) * power / ((2.0 * (double)i + 1.0) * (double)factorial(2 * i + 1)); 
    10472 
    10573                //printf ("Si=%g %g %d\n", x, out, i); 
    106  
    10774        } 
    10875 
    10976        return out; 
    110  
    11177} 
    11278 
     
    138104} 
    139105 
     106/** Modifications below by kieranrcampbell@gmail.com 
     107    Institut Laue-Langevin, July 2012 
     108**/ 
     109 
     110/** 
     111   Implements eq 6.2.5 (small gamma) of Numerical Recipes in C, essentially 
     112   the incomplete gamma function multiplied by the gamma function. 
     113   Required for implementation of fast error function (erf) 
     114**/ 
     115 
     116 
     117#define ITMAX 100 
     118#define EPS 3.0e-7 
     119#define FPMIN 1.0e-30 
     120 
     121void gser(float *gamser, float a, float x, float *gln) { 
     122  int n; 
     123  float sum,del,ap; 
     124 
     125  *gln = gamln(a); 
     126  if(x <= 0.0) { 
     127    if (x < 0.0) printf("Error: x less than 0 in routine gser"); 
     128    *gamser = 0.0; 
     129    return; 
     130  } else { 
     131    ap = a; 
     132    del = sum = 1.0/a; 
     133     
     134    for(n=1;n<=ITMAX;n++) { 
     135      ++ap; 
     136      del *= x/ap; 
     137      sum += del; 
     138      if(fabs(del) < fabs(sum)*EPS) { 
     139        *gamser = sum * exp(-x + a * log(x) - (*gln)); 
     140        return; 
     141      } 
     142    } 
     143    printf("a too large, ITMAX too small in routine gser"); 
     144    return; 
     145 
     146  } 
     147 
     148 
     149} 
     150 
     151/** 
     152   Implements the incomplete gamma function Q(a,x) evaluated by its continued fraction  
     153   representation  
     154**/ 
     155 
     156void gcf(float *gammcf, float a, float x, float *gln) { 
     157  int i; 
     158  float an,b,c,d,del,h; 
     159 
     160  *gln = gamln(a); 
     161  b = x+1.0-a; 
     162  c = 1.0/FPMIN; 
     163  d = 1.0/b; 
     164  h=d; 
     165  for (i=1;i <= ITMAX; i++) { 
     166    an = -i*(i-a); 
     167    b += 2.0; 
     168    d = an*d + b; 
     169    if (fabs(d) < FPMIN) d = FPMIN; 
     170    c = b+an/c; 
     171    if (fabs(c) < FPMIN) c = FPMIN; 
     172    d = 1.0/d; 
     173    del = d*c; 
     174    h += del; 
     175    if (fabs(del-1.0) < EPS) break; 
     176  } 
     177  if (i > ITMAX) printf("a too large, ITMAX too small in gcf"); 
     178  *gammcf = exp(-x+a*log(x)-(*gln))*h; 
     179  return; 
     180} 
     181 
     182 
     183/** 
     184   Represents incomplete error function, P(a,x) 
     185**/ 
     186float gammp(float a, float x) { 
     187  float gamser,gammcf,gln; 
     188  if(x < 0.0 || a <= 0.0) printf("Invalid arguments in routine gammp"); 
     189  if (x < (a+1.0)) { 
     190    gser(&gamser,a,x,&gln); 
     191    return gamser; 
     192  } else { 
     193    gcf(&gammcf,a,x,&gln); 
     194    return 1.0 - gammcf; 
     195  } 
     196} 
     197 
     198/**  
     199    Implementation of the error function, erf(x) 
     200**/ 
     201 
     202float erff(float x) { 
     203  return x < 0.0 ? -gammp(0.5,x*x) : gammp(0.5,x*x); 
     204} 
     205 
  • sansmodels/src/python_wrapper/WrapperGenerator.py

    rfa6db8b r5da3cc5  
    9292        # Model category 
    9393        self.category = None 
     94        # Whether model belongs to multifunc 
     95        self.is_multifunc = False 
    9496        ## output directory for wrappers 
    9597        self.output_dir = output_dir 
     
    228230                    raise ValueError, "Could not parse file %s" % self.file 
    229231 
    230  
     232            # is_multifunc 
     233            key = "[MULTIPLICITY_INFO]" 
     234            if line.count(key) > 0: 
     235                self.is_multifunc = True 
     236                try: 
     237                    index = line.index(key) 
     238                    toks = line[index:].split("=") 
     239                    self.multiplicity_info = toks[1].lstrip().rstrip() 
     240                except: 
     241                    raise ValueError, "Could not parse file %s" % self.file 
    231242 
    232243            # Catch struct name 
     
    331342            newline = self.replaceToken(newline,  
    332343                                        "[MODELSTRUCT]", self.structName) 
     344 
     345            # Sort model initialization based on multifunc 
     346            if(self.is_multifunc): 
     347                line = "int level = 1;\nPyArg_ParseTuple(args,\"i\",&level);\n" 
     348                line += "self->model = new " + self.pythonClass + "(level);" 
     349            else: 
     350                line = "self->model = new " + self.pythonClass + "();" 
     351     
     352            newline = self.replaceToken(newline,"[INITIALIZE_MODEL]", 
     353                                            line) 
    333354             
    334355            # Dictionary initialization 
     
    465486            newline = self.replaceToken(newline,  
    466487                                        "[PAR_DETAILS]", self.details) 
     488             
     489            # Call base constructor 
     490            if self.is_multifunc: 
     491                newline = self.replaceToken(newline,"[CALL_CPYTHON_INIT]", 
     492                                            'C' + self.pythonClass + ".__init__(self,multfactor)\n\tself.is_multifunc = True") 
     493                newline = self.replaceToken(newline,"[MULTIPLICITY_INFO]",self.multiplicity_info) 
     494            else: 
     495                newline = self.replaceToken(newline,"[CALL_CPYTHON_INIT]", 
     496                                            'C' + self.pythonClass + ".__init__(self)\n\tself.is_multifunc = False") 
     497                newline = self.replaceToken(newline,"[MULTIPLICITY_INFO]","None") 
     498 
    467499            
    468500            # fixed list  details 
  • sansmodels/src/python_wrapper/classTemplate.txt

    re08bd5b r5da3cc5  
    8686        self->params = PyDict_New(); 
    8787        self->dispersion = PyDict_New(); 
    88         self->model = new [CMODEL](); 
    89          
     88 
     89        [INITIALIZE_MODEL] 
     90 
    9091        [INITDICTIONARY] 
    9192          
  • sansmodels/src/python_wrapper/modelTemplate.txt

    rfa6db8b r5da3cc5  
    4343    """ 
    4444         
    45     def __init__(self): 
     45    def __init__(self,multfactor=1): 
    4646        """ Initialization """ 
    4747        self.__dict__ = {} 
     
    5050        BaseComponent.__init__(self) 
    5151        #apply([CPYTHONCLASS].__init__, (self,))  
    52         [CPYTHONCLASS].__init__(self) 
    53          
     52 
     53        [CALL_CPYTHON_INIT] 
     54                         
    5455        ## Name of the model 
    5556        self.name = "[PYTHONCLASS]" 
     
    7071 
    7172        self.category = [CATEGORY] 
     73        self.multiplicity_info = [MULTIPLICITY_INFO] 
     74         
    7275 
    7376    def __setstate__(self, state): 
Note: See TracChangeset for help on using the changeset viewer.