psychlops/core/math/psychlops_m_util.cpp

   1 /*
   2  *  psychlops_m_util.cpp
   3  *  Psychlops Standard Library (Universal)
   4  *
   5  *  Last Modified 2006/01/04 by Kenchi HOSOKAWA
   6  *  (C) 2006 Kenchi HOSOKAWA, Kazushi MARUYA and Takao SATO
   7  */
   8
   9
  10 #include <stdlib.h>
  11 #include <Math.h>
  12 #include <time.h>
  13
  14
  15 #if defined(__GNUC__)
  16 #  include "dSFMT2.0/dSFMT.c"
  17 #elif defined(_MSC_VER)
  18 #  include "dSFMT2.0/dSFMT.c"
  19 #elif defined(__BORLANDC__)
  20 #  include "dSFMT2.0/dSFMT.c"
  21 #endif
  22
  23 #include "dSFMT2.0/dSFMT.h"
  24 #include "okumura/igamma.h"
  25
  26
  27
  28 #include "../ApplicationInterfaces/psychlops_code.h"
  29 #include "psychlops_m_util.h"
  30 #include "psychlops_m_matrix.h"
  31
  32
  33 namespace Psychlops {
  34
  35
  36         const double PI = 3.14159265358979;\r
  37         const double Math::E = 2.718281828459045;\r
  38         const double Math::LOG2E = 1.44269504088896340736;
  39
  40
  41         /*              Functions               */
  42
  43         double random() {
  44                 return dsfmt_gv_genrand_close_open();
  45         }
  46         double random(double factor) {
  47                 return dsfmt_gv_genrand_close_open()*factor;
  48         }
  49         int random(int factor) {
  50                 return (int)(dsfmt_gv_genrand_close_open()*factor);
  51         }
  52         double random(double min, double max) {
  53                 double dif = (max-min);
  54                 if(dif<0) throw Exception("min exceeds max in random(min, max)");
  55                 return random(dif)+min;
  56         }
  57         void random(double array[], int size) {
  58                 for(int i=0; i<size; i++) array[i]=random();
  59         }
  60         void random(double array[], int size, double factor) {
  61                 for(int i=0; i<size; i++) array[i]=random(factor);
  62         }
  63         void random(double array[], int size, double min, double max) {
  64                 double factor = (max-min);
  65                 if(factor<0) throw Exception("min exceeds max in random(min, max)");
  66                 for(int i=0; i<size; i++) array[i]=random(factor)+min;
  67         }
  68         void random(Matrix &mtx) {
  69                 double *array, *str; int step; mtx.get_element_ptr(array, str, step);
  70                 int size = mtx.getRows()*mtx.getCols();
  71                 for(int i=0; i<size; i++) array[i]=random();
  72         }
  73         void random(Matrix &mtx, double factor) {
  74                 double *array, *str; int step; mtx.get_element_ptr(array, str, step);
  75                 int size = mtx.getRows()*mtx.getCols();
  76                 for(int i=0; i<size; i++) array[i]=random(factor);
  77         }
  78         void random(Matrix &mtx, double min, double max) {
  79                 double *array, *str; int step; mtx.get_element_ptr(array, str, step);
  80                 int size = mtx.getRows()*mtx.getCols();
  81                 double factor = (max-min);
  82                 if(factor<0) throw Exception("min exceeds max in random(min, max)");
  83                 for(int i=0; i<size; i++) array[i]=random(factor)+min;
  84         }
  85         void randomize(unsigned long seed) {
  86                 srand( (int)seed );
  87                 dsfmt_gv_init_gen_rand( seed );
  88 //              for(int i=0; i<1000000; i++) random();
  89         }
  90         void randomize() {
  91                 srand( (int)time(NULL) );
  92                 dsfmt_gv_init_gen_rand( (unsigned long)time(NULL) );
  93 //              for(int i=0; i<1000000; i++) random();
  94         }
  95
  96         void Math::setscramble(/*unsigned*/ int array[], unsigned int length) {
  97                 unsigned int rnd;
  98                 for(unsigned int i=0;i<length;i++) { array[i] = -1; }
  99                 for(unsigned int i=0;i<length;i++) {
 100                         rnd = rand()%length;
 101                         while(array[rnd]!=-1) {
 102                                 rnd++;
 103                                 if(rnd>=length) { rnd=0; }
 104                         }
 105                         array[rnd] = i;
 106                 }
 107         }
 108
 109
 110         int Math::round(double val) {
 111                 double integer_part, particle = modf(val, &integer_part);
 112                 return ( ( particle<0.5 | (particle==0.5 && (int)integer_part%2==0) ) ? (int)integer_part : (int)integer_part+1 );
 113         }\r
 114         double Math::log2(double val) {\r
 115                 return log(val)*LOG2E;\r
 116         }\r
 117         /*\r
 118         double Math::max(double val1, double val2) {\r
 119                 return val1>val2 ? val1 : val2;\r
 120         }\r
 121         double Math::min(double val1, double val2) {\r
 122                 return val1<val2 ? val1 : val2;\r
 123         }\r
 124         */\r
 125
 126 \r
 127         double Math::gaussian(double x, double sigma) {\r
 128                 return exp( -(x*x) / (2.0*sigma*sigma) );\r
 129         }\r
 130
 131         double Math::normalDistibution(double x, double mu, double sigma) {
 132                 return exp( -( (x-mu)*(x-mu) / (2*sigma*sigma) ) ) / sqrt(2*PI*sigma*sigma);
 133         }
 134
 135         double Math::cumulativeNormalDistibution(double x, double mu, double sigma) {
 136                 return .5 + .5*erf( (x-mu)/(sigma*sqrt(2.0) ) );
 137         }
 138
 139 /*
 140         double Math::mod(double lhs, double rhs) {
 141                 return lhs - floor(lhs/rhs)*rhs;
 142         }
 143
 144         double Math::abs(double val) {
 145                 return (val<0 ? -val : val);
 146         }
 147 */
 148
 149
 150 }       /*      <- namespace Psycholops         */