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