dev4/psychlops/extention/math/solver.cs

   1 using System.Threading;\r
   2 using System;\r
   3 using BigNum;\r
   4 \r
   5 namespace Psychlops\r
   6 {\r
   7 \r
   8         namespace Solver\r
   9         {\r
  10 \r
  11                 \r
  12                 static public class Constants\r
  13                 {\r
  14                         public static readonly int LIMIT = 5;\r
  15                         public delegate double Func0(double x);\r
  16                         public delegate double Func1(double x, double a);\r
  17                         public delegate double Func2(double x, double a, double b);\r
  18                         public delegate double Func3(double x, double a, double b, double c);\r
  19                         public delegate double Func4(double x, double a, double b, double c, double d);\r
  20                         public delegate double Func5(double x, double a, double b, double c, double d, double e);\r
  21                 }\r
  22 \r
  23 \r
  24                 static public class Binomial\r
  25                 {\r
  26                         public static double factorial(int fact)\r
  27                         {\r
  28                                 double x = 1;\r
  29                                 for (int i = 1; i <= fact; i++)\r
  30                                 {\r
  31                                         x *= i;\r
  32                                 }\r
  33                                 return x;\r
  34                         }\r
  35                         public static double combination(int nn, int kk)\r
  36                         {\r
  37                                 //return (double)(Int64)(factorial(n) / (factorial(kk) * factorial(n - kk)));\r
  38                                 BigInt n = new BigInt((Int64)nn, new PrecisionSpec());\r
  39                                 BigInt k = new BigInt((Int64)kk, new PrecisionSpec());\r
  40                                 BigInt n_k = n - k;\r
  41                                 n.Factorial();\r
  42                                 k.Factorial();\r
  43                                 n_k.Factorial();\r
  44                                 return (double)(Int64)(n / (k * n_k) );\r
  45                         }\r
  46                         public static double likelihood(double pp, int yes, int no)\r
  47                         {\r
  48                                 double p;\r
  49                                 if (pp < 0) p = 0; else if (pp > 1) p = 1; else p = pp;\r
  50                                 return combination(yes + no, yes) * System.Math.Pow(p, yes) * System.Math.Pow(1 - p, no);\r
  51                         }\r
  52 \r
  53                         public static double likelihood(double pp, int yes, int no, double comb)\r
  54                         {\r
  55                                 double p;\r
  56                                 if (pp < 0) p = 0; else if (pp > 1) p = 1; else p = pp;\r
  57                                 return comb * System.Math.Pow(p, yes) * System.Math.Pow(1 - p, no);\r
  58                         }\r
  59 \r
  60                 }\r
  61 \r
  62 \r
  63                 public class BernoulliProcess\r
  64                 {\r
  65                         public struct Data\r
  66                         {\r
  67                                 public double x;\r
  68                                 public int pos, neg, comb;\r
  69                                 public void set(double cond, int posit, int negat)\r
  70                                 {\r
  71                                         x = cond;\r
  72                                         pos = posit;\r
  73                                         neg = negat;\r
  74                                 }\r
  75 \r
  76                         };\r
  77                         public Data[] elems;\r
  78                         public int length;\r
  79                         public void set(double[][] num)\r
  80                         {\r
  81                                 elems = new Data[num.GetLength(0)];\r
  82                                 for (int i = 0; i < elems.GetLength(0); i++)\r
  83                                 {\r
  84                                         elems[i].set(num[i][0], (int)num[i][1], (int)num[i][2]);\r
  85                                 }\r
  86                         }\r
  87                 };\r
  88 \r
  89 \r
  90 \r
  91                 class BinomialLikelihoodThread\r
  92                 {\r
  93                         public bool finished;\r
  94                         public static bool started;\r
  95                         internal Thread th;\r
  96                         public static BinomialLikelihoodThread current;\r
  97 \r
  98                         public Interval[] itvl;\r
  99                         public double[] step;\r
 100                         public double[] champ;\r
 101                         public double champ_like;\r
 102                         public BernoulliProcess data;\r
 103 \r
 104 \r
 105                         public Constants.Func0 func0;\r
 106                         public Constants.Func1 func1;\r
 107                         public Constants.Func2 func2;\r
 108                         public Constants.Func3 func3;\r
 109                         public Constants.Func4 func4;\r
 110                         public Constants.Func5 func5;\r
 111 \r
 112 \r
 113                         public BinomialLikelihoodThread()\r
 114                         {\r
 115                         }\r
 116                         public void waitLoop()\r
 117                         {\r
 118                                 finished = false;\r
 119                                 for(int i=0; i<10; i++) {\r
 120                                         champ[i] = 0;\r
 121                                 }\r
 122                                 current = this;\r
 123                                 started = false;\r
 124                         }\r
 125 \r
 126                         public void loop1() { waitLoop(); th = new Thread(new ThreadStart(loop1_)); }\r
 127                         public void loop2() { waitLoop(); th = new Thread(new ThreadStart(loop2_)); }\r
 128                         public void loop3() { waitLoop(); th = new Thread(new ThreadStart(loop3_)); }\r
 129                         public void loop4() { waitLoop(); th = new Thread(new ThreadStart(loop4_)); }\r
 130                         public void loop5() { waitLoop(); th = new Thread(new ThreadStart(loop5_)); }\r
 131                         void loop1_()\r
 132                         {\r
 133                                 started = true;\r
 134                                 double p;\r
 135                                 double like = 1.0;\r
 136                                 champ_like=0.0;\r
 137                                 int L = data.length;\r
 138                                 for (double a = itvl[0].begin.val; a < itvl[0].end.val; a += step[0])\r
 139                                 {\r
 140                                         like = 1.0;\r
 141                                         for(int i=0; i<L; i++) {\r
 142                                                 p = func1(data.elems[i].x, a);\r
 143                                                 like *= Binomial.likelihood(p, data.elems[i].pos, data.elems[i].neg, data.elems[i].comb);\r
 144                                         }\r
 145                                         if(like > champ_like) {\r
 146                                                 champ_like = like;\r
 147                                                 champ[0] = a;\r
 148                                         }\r
 149                                 }\r
 150                                 finished = true;\r
 151                         }\r
 152                         void loop2_()\r
 153                         {\r
 154                                 started = true;\r
 155                                 double p;\r
 156                                 double like = 1.0;\r
 157                                 champ_like = 0.0;\r
 158                                 int L = data.length;\r
 159                                 for (double a = itvl[0].begin.val; a < itvl[0].end.val; a += step[0])\r
 160                                 {\r
 161                                         for (double b = itvl[1].begin.val; b < itvl[1].end.val; b += step[1])\r
 162                                         {\r
 163                                                 like = 1.0;\r
 164                                                 for (int i = 0; i < L; i++)\r
 165                                                 {\r
 166                                                         p = func2(data.elems[i].x, a, b);\r
 167                                                         like *= Binomial.likelihood(p, data.elems[i].pos, data.elems[i].neg, data.elems[i].comb);\r
 168                                                 }\r
 169                                                 if (like > champ_like)\r
 170                                                 {\r
 171                                                         champ_like = like;\r
 172                                                         champ[0] = a;\r
 173                                                 }\r
 174                                         }\r
 175                                 }\r
 176                                 finished = true;\r
 177                         }\r
 178                         void loop3_()\r
 179                         {\r
 180                                 started = true;\r
 181                                 double p;\r
 182                                 double like = 1.0;\r
 183                                 champ_like = 0.0;\r
 184                                 int L = data.length;\r
 185                                 for (double a = itvl[0].begin.val; a < itvl[0].end.val; a += step[0])\r
 186                                 {\r
 187                                         for (double b = itvl[1].begin.val; b < itvl[1].end.val; b += step[1])\r
 188                                         {\r
 189                                                 for (double c = itvl[2].begin.val; c < itvl[2].end.val; c += step[2])\r
 190                                                 {\r
 191                                                         like = 1.0;\r
 192                                                         for (int i = 0; i < L; i++)\r
 193                                                         {\r
 194                                                                 p = func3(data.elems[i].x, a, b, c);\r
 195                                                                 like *= Binomial.likelihood(p, data.elems[i].pos, data.elems[i].neg, data.elems[i].comb);\r
 196                                                         }\r
 197                                                         if (like > champ_like)\r
 198                                                         {\r
 199                                                                 champ_like = like;\r
 200                                                                 champ[0] = a;\r
 201                                                         }\r
 202                                                 }\r
 203                                         }\r
 204                                 }\r
 205                                 finished = true;\r
 206                         }\r
 207                         void loop4_()\r
 208                         {\r
 209                                 started = true;\r
 210                                 finished = true;\r
 211                         }\r
 212                         void loop5_()\r
 213                         {\r
 214                                 started = true;\r
 215                                 finished = true;\r
 216                         }\r
 217                 };\r
 218 \r
 219 \r
 220 \r
 221                 class BinomialLikelihood\r
 222                 {\r
 223                         public int iteration;\r
 224 \r
 225                         public Interval[] itvl;\r
 226                         public double[] step;\r
 227                         public double[] champ;\r
 228                         public double champ_like;\r
 229                         public BernoulliProcess data;\r
 230 \r
 231                         public Constants.Func0 func0;\r
 232                         public Constants.Func1 func1;\r
 233                         public Constants.Func2 func2;\r
 234                         public Constants.Func3 func3;\r
 235                         public Constants.Func4 func4;\r
 236                         public Constants.Func5 func5;\r
 237 \r
 238                         public BinomialLikelihood()\r
 239                         {\r
 240                                 iteration = 2;\r
 241                         }\r
 242 \r
 243                         public void begin(Constants.Func0 d_func)\r
 244                         { }\r
 245                         public void begin(Constants.Func1 d_func)\r
 246                         {\r
 247                                 func1 = d_func;\r
 248 \r
 249                                 BinomialLikelihoodThread[] l = new BinomialLikelihoodThread[4];\r
 250 \r
 251                                 double r, low, high;\r
 252                                 for (int k = 0; k < iteration; k++)\r
 253                                 {\r
 254                                         champ_like = 0;\r
 255 \r
 256                                         begin_base(l);\r
 257                                         for (int i = 0; i < 4; i++)\r
 258                                         {\r
 259                                                 l[i].data = data;\r
 260                                                 l[i].func1 = d_func;\r
 261                                                 l[i].loop1();\r
 262                                         }\r
 263                                         end_base(l);\r
 264 \r
 265                                         for (int j = 0; j < Constants.LIMIT; j++)\r
 266                                         {\r
 267                                                 r = itvl[j].end.val - itvl[j].begin.val;\r
 268                                                 low = champ[j] - r / 8.0 < itvl[j].begin.val ? itvl[j].begin.val : champ[j] - r / 8.0;\r
 269                                                 high = champ[j] + r / 8.0 > itvl[j].end.val ? itvl[j].end.val : champ[j] + r / 8.0;\r
 270                                                 itvl[j] = new Interval(low, high);\r
 271                                         }\r
 272                                 }\r
 273                         }\r
 274                         public void begin(Constants.Func2 d_func)\r
 275                         {\r
 276                                 func2 = d_func;\r
 277 \r
 278                                 BinomialLikelihoodThread[] l = new BinomialLikelihoodThread[4];\r
 279 \r
 280                                 double r, low, high;\r
 281                                 for (int k = 0; k < iteration; k++)\r
 282                                 {\r
 283                                         champ_like = 0;\r
 284 \r
 285                                         begin_base(l);\r
 286                                         for (int i = 0; i < 4; i++)\r
 287                                         {\r
 288                                                 l[i].data = data;\r
 289                                                 l[i].func2 = d_func;\r
 290                                                 l[i].loop2();\r
 291                                         }\r
 292                                         end_base(l);\r
 293 \r
 294                                         for (int j = 0; j < Constants.LIMIT; j++)\r
 295                                         {\r
 296                                                 r = itvl[j].end.val - itvl[j].begin.val;\r
 297                                                 low = champ[j] - r / 8.0 < itvl[j].begin.val ? itvl[j].begin.val : champ[j] - r / 8.0;\r
 298                                                 high = champ[j] + r / 8.0 > itvl[j].end.val ? itvl[j].end.val : champ[j] + r / 8.0;\r
 299                                                 itvl[j] = new Interval(low, high);\r
 300                                         }\r
 301                                 }\r
 302                         }\r
 303                         public void begin(Constants.Func3 d_func)\r
 304                         {\r
 305                                 func3 = d_func;\r
 306 \r
 307                                 BinomialLikelihoodThread[] l = new BinomialLikelihoodThread[4];\r
 308 \r
 309                                 double r, low, high;\r
 310                                 for (int k = 0; k < iteration; k++)\r
 311                                 {\r
 312                                         champ_like = 0;\r
 313 \r
 314                                         begin_base(l);\r
 315                                         for(int i=0; i<4; i++) {\r
 316                                                 l[i].data = data;\r
 317                                                 l[i].func3 = d_func;\r
 318                                                 l[i].loop3();\r
 319                                         }\r
 320                                         end_base(l);\r
 321 \r
 322                                         for (int j = 0; j < Constants.LIMIT; j++)\r
 323                                         {\r
 324                                                 r = itvl[j].end.val - itvl[j].begin.val;\r
 325                                                 low  = champ[j]-r/8.0 < itvl[j].begin.val ? itvl[j].begin.val : champ[j]-r/8.0;\r
 326                                                 high = champ[j]+r/8.0 > itvl[j].end.val   ? itvl[j].end.val   : champ[j]+r/8.0;\r
 327                                                 itvl[j] = new Interval(low, high);\r
 328                                         }\r
 329                                 }\r
 330                         }\r
 331 \r
 332                         public void begin(Constants.Func4 d_func)\r
 333                         {\r
 334                         }\r
 335                         public void begin(Constants.Func5 d_func)\r
 336                         {\r
 337                         }\r
 338 \r
 339                         void begin_base(BinomialLikelihoodThread[] l)\r
 340                         {\r
 341                                 data.length = data.elems.GetLength(0);\r
 342                                 for (int i = 0; i < data.elems.GetLength(0); i++)\r
 343                                 {\r
 344                                         data.elems[i].comb = (int)Binomial.combination(data.elems[i].pos + data.elems[i].neg, data.elems[i].pos);\r
 345                                 }\r
 346 \r
 347                                 for (int j = 0; j < Constants.LIMIT; j++) { step[j] = (itvl[j].end.val - itvl[j].begin.val) / 256.0; }\r
 348                                 for (int i = 0; i < 4; i++)\r
 349                                 {\r
 350                                         l[i].itvl[0] = new Interval((itvl[0].begin.val) + (i * step[0] * 64), (itvl[0].begin.val) + ((i + 1) * step[0] * 64));\r
 351                                         l[i].step[0] = step[0];\r
 352                                         for (int j = 1; j < Constants.LIMIT; j++)\r
 353                                         {\r
 354                                                 l[i].itvl[j] = itvl[j];\r
 355                                                 l[i].step[j] = step[j];\r
 356                                         }\r
 357                                         l[i].data = data;\r
 358                                 }\r
 359                         }\r
 360 \r
 361                         void end_base(BinomialLikelihoodThread[] l)\r
 362                         {\r
 363                                 for (int i = 0; i < 4; i++)\r
 364                                 {\r
 365                                         l[i].th.Join();\r
 366                                 }\r
 367 \r
 368                                 for(int j=0; j<Constants.LIMIT; j++) { champ[j] = 0; }\r
 369                                 for(int i=0; i<4; i++) {\r
 370                                         if(champ_like < l[i].champ_like) {\r
 371                                                 champ_like = l[i].champ_like;\r
 372                                                 for(int j=0; j<Constants.LIMIT; j++) { champ[j] = l[i].champ[j]; }\r
 373                                         }\r
 374                                 }\r
 375 \r
 376                         }\r
 377 \r
 378 \r
 379                 }\r
 380 \r
 381 \r
 382         }\r
 383 \r
 384 }