add statistical functions to numparam and nutmeg parsers

ChangeLog

@@ -1,3 +1,12 @@
+2010-12-28 Holger Vogt
+ * xpressn.c, cmath2.c, cmath2.h, randnumb.c, parse.c,
+   examples/Monte_Carlo/MonteCarlo.sp:
+    add poisson and exponential distribution to nutmeg parser,
+    add gauss, aunif, unif, limit to numparam parser
+ * spiceif.c: remove bug in experimental_code	
+ * control.c: no i/o redirection in define command 
+   (so > or < may be used safely)
+
 2010-12-27  Robert Larice
  * tests/bsim3soi/inv_dc.cir ,
  * tests/bsim3soi/inv_tr.cir ,

examples/Monte_Carlo/MonteCarlo.sp

@@ -17,30 +17,34 @@ R2 0 OUT 141
   set scratch=$curplot     $ store its name to 'scratch'
   setplot $scratch         $ make 'scratch' the active plot 
   let bwh=unitvec(mc_runs) $ create a vector in plot 'scratch' to store bandwidth data
-*
-  dowhile run < mc_runs    $ loop starts here
-*
-*    alter c1 = unif(1e-09, 0.1)
-*    alter c1 = aunif(1e-09, 100e-12)
-*    alter c1 = gauss(1e-09, 0.1, 3)
-*    alter c1 = agauss(1e-09, 100e-12, 3)
-*
+  
 * define distributions for random numbers:
 * unif: uniform distribution, deviation relativ to nominal value
 * aunif: uniform distribution, deviation absolut
 * gauss: Gaussian distribution, deviation relativ to nominal value
 * agauss: Gaussian distribution, deviation absolut
-  define unif(nom, var) (nom + (nom*var) * sunif(0))
+* limit: if unif. distributed value >=0 then add +avar to nom, else -avar
+  define unif(nom, rvar) (nom + (nom*rvar) * sunif(0))
   define aunif(nom, avar) (nom + avar * sunif(0))
-  define gauss(nom, var, sig) (nom + (nom*var)/sig * sgauss(0))
+  define gauss(nom, rvar, sig) (nom + (nom*rvar)/sig * sgauss(0))
   define agauss(nom, avar, sig) (nom + avar/sig * sgauss(0))
+*  define limit(nom, avar) (nom + ((sgauss(0) ge 0) ? avar : -avar))
+  define limit(nom, avar) (nom + ((sgauss(0) >= 0) ? avar : -avar))
+*
+*
+  dowhile run < mc_runs    $ loop starts here
+*
+*    alter c1 = unif(1e-09, 0.1)
+*    alter c1 = aunif(1e-09, 100e-12)
+*    alter c1 = gauss(1e-09, 0.1, 3)
+*    alter c1 = agauss(1e-09, 100e-12, 3)
 *
     alter c1 = unif(1e-09, 0.1)
     alter l1 = unif(10e-06, 0.1)
     alter c2 = unif(1e-09, 0.1)
     alter l2 = unif(10e-06, 0.1)
     alter l3 = unif(40e-06, 0.1)
-    alter c3 = unif(250e-12, 0.1)
+    alter c3 = limit(250e-12, 25e-12)
 *
     ac oct 100 250K 10Meg
 *	

src/frontend/control.c

@@ -56,7 +56,8 @@ int stackp = 0;
  * destroyed safely
  */
 
-static char *noredirect[] = { "stop", NULL } ;  /* Only one?? */
+/* no redirection after the following commands (we may need more to add here!) */
+static char *noredirect[] = { "stop", "define", NULL } ;
 
 
 static struct control *

src/frontend/numparam/xpressn.c

@@ -17,7 +17,8 @@
 #include "compatmode.h"
 
 /* random numbers in /maths/misc/randnumb.c */
-extern double gauss(void);
+extern double gauss0(void);
+extern double drand(void);
 
 void debugwarn(tdico *d, char *s);
 
@@ -46,11 +47,37 @@ ternary_fcn (int conditional, double if_value, double else_value)
 
 
 static double
-agauss (double nominal_val, double variation, double sigma)
+agauss (double nominal_val, double abs_variation, double sigma)
 {
   double stdvar;
-  stdvar=variation/sigma;
-  return (nominal_val+stdvar*gauss());
+  stdvar=abs_variation/sigma;
+  return (nominal_val+stdvar*gauss0());
+}
+
+static double
+gauss (double nominal_val, double rel_variation, double sigma)
+{
+  double stdvar;
+  stdvar=nominal_val * rel_variation / sigma;
+  return (nominal_val + stdvar * gauss0());
+}
+
+static double
+unif (double nominal_val, double rel_variation)
+{
+  return (nominal_val + nominal_val * rel_variation * drand());
+}
+
+static double
+aunif (double nominal_val, double abs_variation)
+{
+  return (nominal_val + abs_variation * drand());
+}
+
+static double
+limit (double nominal_val, double abs_variation)
+{
+    return (nominal_val + (drand() > 0 ? abs_variation : -1. * abs_variation));
 }
 
 static void
@@ -62,7 +89,8 @@ initkeys (void)
             "and or not div mod if else end while macro funct defined"
             " include for to downto is var");
   scopy_up (&fmathS,
-            "sqr sqrt sin cos exp ln arctan abs pow pwr max min int log sinh cosh tanh ternary_fcn v agauss sgn");
+            "sqr sqrt sin cos exp ln arctan abs pow pwr max min int log sinh cosh"
+            " tanh ternary_fcn v agauss sgn gauss unif aunif limit");
 }
 
 static double
@@ -1253,9 +1281,16 @@ formula (tdico * dico, char *s, bool *perror)
                       u = ternary_fcn ((int) v, w, u);
                   else if ((fu == 20))
                       u = agauss (v, w, u);
+                  else if ((fu == 22))
+                      u = gauss (v, w, u);
+                  else if ((fu == 23))
+                      u = unif (v, u);
+                  else if ((fu == 24))
+                      u = aunif (v, u);
+                  else if ((fu == 25))
+                      u = limit (v, u);
                   else
                       u = mathfunction (fu, v, u);
-
                 }
             }
           i = k;

src/frontend/parse.c

@@ -164,6 +164,8 @@ struct func ft_funcs[] = {
         { "norm",   cx_norm } ,
         { "rnd",    cx_rnd } ,
         { "sunif",  cx_sunif } ,
+        { "poisson",cx_poisson } ,
+        { "exponential", cx_exponential } ,
         { "sgauss", cx_sgauss } ,
         { "pos",    cx_pos } ,
         { "mean",   cx_mean } ,

src/frontend/spiceif.c

@@ -1513,8 +1513,10 @@ do {\
     	_foo(lname,TRANan,1);
     }
     ft_curckt->ci_curTask->jobs->JOBname = NULL;
-    ckt->CKTcurJob = (&(ft_curckt->ci_curTask->taskOptions)) -> jobs;
-    
+//    ckt->CKTcurJob = (&(ft_curckt->ci_curTask->taskOptions)) -> jobs;
+    ckt->CKTcurJob=ft_curckt->ci_curTask->jobs;
+
+
     _foo(ft_curckt->ci_curTask->jobs->JOBname, char, -1);
 
     ft_curckt->ci_curTask->jobs->JOBnextJob = NULL;

src/maths/cmaths/cmath2.c

@@ -32,7 +32,9 @@ extern void srandom (unsigned int seed);
 */
 extern void checkseed(void); /* seed random or set by 'set rndseed=value'*/
 extern double drand(void);  /* from randnumb.c */
-extern double gauss(void);  /* from randnumb.c */
+extern double gauss0(void);  /* from randnumb.c */
+extern int poisson(double);  /* from randnumb.c */
+extern double exprand(double);  /* from randnumb.c */
 
 static double *
 d_tan(double *dd, int length)
@@ -209,43 +211,52 @@ cx_uminus(void *data, short int type, int length, int *newlength, short int *new
     }
 }
 
+
+/* random integers drawn from a uniform distribution 
+   *data in: integer numbers, their absolut values are used, 
+       maximum is RAND_MAX (32767) 
+   *data out: random integers in interval [0, data[i][ 
+       standard library function rand() is used
+*/
 void *
 cx_rnd(void *data, short int type, int length, int *newlength, short int *newtype)
 {
     *newlength = length;
     checkseed();
     if (type == VF_COMPLEX) {
-      ngcomplex_t *c;
-      ngcomplex_t *cc = (ngcomplex_t *) data;
-      int i;
+        ngcomplex_t *c;
+        ngcomplex_t *cc = (ngcomplex_t *) data;
+        int i;
 
-      c = alloc_c(length);
-      *newtype = VF_COMPLEX;
-      for (i = 0; i < length; i++) {
-        int j, k;
-        j = (int)floor(realpart(&cc[i]));
-        k = (int)floor(imagpart(&cc[i]));
-        realpart(&c[i]) = j ? rand() % j : 0; //random() % j : 0;
-        imagpart(&c[i]) = k ? rand() % k : 0; //random() % k : 0;
-      }
-      return ((void *) c);
+        c = alloc_c(length);
+        *newtype = VF_COMPLEX;
+        for (i = 0; i < length; i++) {
+            int j, k;
+            j = (int)floor(realpart(&cc[i]));
+            k = (int)floor(imagpart(&cc[i]));
+            realpart(&c[i]) = j ? rand() % j : 0;
+            imagpart(&c[i]) = k ? rand() % k : 0;
+        }
+        return ((void *) c);
     } else {
-      double *d;
-      double *dd = (double *) data;
-      int i;
+        double *d;
+        double *dd = (double *) data;
+        int i;
 
-      d = alloc_d(length);
-      *newtype = VF_REAL;
-      for (i = 0; i < length; i++) {
-      int j;
-
-      j = (int)floor(dd[i]);
-      d[i] = j ? rand() % j : 0; //random() % j : 0;
+        d = alloc_d(length);
+        *newtype = VF_REAL;
+        for (i = 0; i < length; i++) {
+            int j;
+            j = (int)floor(dd[i]);
+            d[i] = j ? rand() % j : 0;
+        }
+        return ((void *) d);
     }
-    return ((void *) d);
-  }
 }
 
+/* random numbers drawn from a uniform distribution 
+   *data out: random numbers in interval [-1, 1[ 
+*/
 void *
 cx_sunif(void *data, short int type, int length, int *newlength, short int *newtype)
 {
@@ -254,29 +265,108 @@ cx_sunif(void *data, short int type, int length, int *newlength, short int *newt
     *newlength = length;
     checkseed();
     if (type == VF_COMPLEX) {
-      ngcomplex_t *c;
-      int i;
+        ngcomplex_t *c;
+        int i;
 
-      c = alloc_c(length);
-      *newtype = VF_COMPLEX;
-      for (i = 0; i < length; i++) {
-        realpart(&c[i]) = drand();
-        imagpart(&c[i]) = drand();
-      }
-      return ((void *) c);
+        c = alloc_c(length);
+        *newtype = VF_COMPLEX;
+        for (i = 0; i < length; i++) {
+            realpart(&c[i]) = drand();
+            imagpart(&c[i]) = drand();
+        }
+        return ((void *) c);
     } else {
-      double *d;
-      int i;
+        double *d;
+        int i;
 
-      d = alloc_d(length);
-      *newtype = VF_REAL;
-      for (i = 0; i < length; i++) {
-        d[i] = drand();
+        d = alloc_d(length);
+        *newtype = VF_REAL;
+        for (i = 0; i < length; i++) {
+            d[i] = drand();
+        }
+        return ((void *) d);
     }
-    return ((void *) d);
-  }
 }
 
+/* random numbers drawn from a poisson distribution 
+   *data in:  lambda
+   *data out: random integers according to poisson distribution, 
+       with lambda given by each vector element 
+*/
+void *
+cx_poisson(void *data, short int type, int length, int *newlength, short int *newtype)
+{
+    NG_IGNORE(data);
+
+    *newlength = length;
+    checkseed();
+    if (type == VF_COMPLEX) {
+        ngcomplex_t *c;
+        ngcomplex_t *cc = (ngcomplex_t *) data;
+        int i;
+
+        c = alloc_c(length);
+        *newtype = VF_COMPLEX;
+        for (i = 0; i < length; i++) {
+            realpart(&c[i]) = poisson (realpart(&cc[i]));
+            imagpart(&c[i]) = poisson (imagpart(&cc[i]));
+        }
+        return ((void *) c);
+    } else {
+        double *d;
+        double *dd = (double *) data;
+        int i;
+
+        d = alloc_d(length);
+        *newtype = VF_REAL;
+        for (i = 0; i < length; i++) {
+            d[i] = poisson(dd[i]);
+        }
+        return ((void *) d);
+    }
+}
+
+/* random numbers drawn from an exponential distribution 
+   *data in:  Mean values
+   *data out: exponentially distributed random numbers, 
+       with mean given by each vector element 
+*/
+void *
+cx_exponential(void *data, short int type, int length, int *newlength, short int *newtype)
+{
+    NG_IGNORE(data);
+
+    *newlength = length;
+    checkseed();
+    if (type == VF_COMPLEX) {
+        ngcomplex_t *c;
+        ngcomplex_t *cc = (ngcomplex_t *) data;
+        int i;
+
+        c = alloc_c(length);
+        *newtype = VF_COMPLEX;
+        for (i = 0; i < length; i++) {
+            realpart(&c[i]) = exprand(realpart(&cc[i]));
+            imagpart(&c[i]) = exprand(imagpart(&cc[i]));
+        }
+        return ((void *) c);
+    } else {
+        double *d;
+        double *dd = (double *) data;
+        int i;
+
+        d = alloc_d(length);
+        *newtype = VF_REAL;
+        for (i = 0; i < length; i++) {
+            d[i] = exprand(dd[i]);
+        }
+        return ((void *) d);
+    }
+}
+
+/* random numbers drawn from a Gaussian distribution 
+   mean 0, std dev 1
+*/
 void *
 cx_sgauss(void *data, short int type, int length, int *newlength, short int *newtype)
 {
@@ -285,29 +375,32 @@ cx_sgauss(void *data, short int type, int length, int *newlength, short int *new
     *newlength = length;
     checkseed();
     if (type == VF_COMPLEX) {
-      ngcomplex_t *c;
-      int i;
+        ngcomplex_t *c;
+        int i;
 
-      c = alloc_c(length);
-      *newtype = VF_COMPLEX;
-      for (i = 0; i < length; i++) {
-        realpart(&c[i]) = gauss();
-        imagpart(&c[i]) = gauss();
-      }
-      return ((void *) c);
+        c = alloc_c(length);
+        *newtype = VF_COMPLEX;
+        for (i = 0; i < length; i++) {
+            realpart(&c[i]) = gauss0();
+            imagpart(&c[i]) = gauss0();
+        }
+        return ((void *) c);
     } else {
-      double *d;
-      int i;
+        double *d;
+        int i;
 
-      d = alloc_d(length);
-      *newtype = VF_REAL;
-      for (i = 0; i < length; i++) {
-        d[i] = gauss();
+        d = alloc_d(length);
+        *newtype = VF_REAL;
+        for (i = 0; i < length; i++) {
+            d[i] = gauss0();
+        }
+        return ((void *) d);
     }
-    return ((void *) d);
-  }
 }
 
+
+
+
 /* Compute the avg of a vector.
    Created by A.M.Roldan 2005-05-21  */
    

src/maths/cmaths/cmath2.h

@@ -14,6 +14,8 @@ void * cx_uminus(void *data, short int type, int length, int *newlength, short i
 void * cx_rnd(void *data, short int type, int length, int *newlength, short int *newtype);
 void * cx_sunif(void *data, short int type, int length, int *newlength, short int *newtype);
 void * cx_sgauss(void *data, short int type, int length, int *newlength, short int *newtype);
+void * cx_poisson(void *data, short int type, int length, int *newlength, short int *newtype);
+void * cx_exponential(void *data, short int type, int length, int *newlength, short int *newtype);
 void * cx_mean(void *data, short int type, int length, int *newlength, short int *newtype);
 void * cx_length(void *data, short int type, int length, int *newlength, short int *newtype);
 void * cx_vector(void *data, short int type, int length, int *newlength, short int *newtype);

src/maths/misc/randnumb.c

@@ -69,7 +69,8 @@ double exprand(double);
 
 void checkseed(void);
 double drand(void);
-double gauss(void);
+double gauss0(void);
+int poisson(double);
 
 
 /* Check if a seed has been set by the command 'set rndseed=value'
@@ -197,7 +198,7 @@ unsigned int CombLCGTausInt2(void)
 
 /***  gauss  ***/
 
-double gauss(void)
+double gauss0(void)
 {
   static bool gliset = TRUE;
   static double glgset = 0.0;
@@ -241,6 +242,26 @@ void rgauss(double* py1, double* py2)
 }	
 
 
+
+/** Code by: Inexpensive
+    http://everything2.com/title/Generating+random+numbers+with+a+Poisson+distribution **/
+int poisson(double lambda)
+{
+  int k=0;                          //Counter
+  const int max_k = 1000;           //k upper limit
+  double p = CombLCGTaus();         //uniform random number
+  double P = exp(-lambda);        //probability
+  double sum=P;                     //cumulant
+  if (sum>=p) return 0;             //done allready
+  for (k=1; k<max_k; ++k) {         //Loop over all k:s
+    P*=lambda/(double)k;           //Calc next prob
+    sum+=P;                         //Increase cumulant
+    if (sum>=p) break;              //Leave loop
+  }
+  return k;                         //return random number
+}
+
+
 /* return an exponentially distributed random number */
 double exprand(double mean)
 {