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