-using System.Threading;\r
-using System;\r
-using BigNum;\r
-\r
-namespace Psychlops\r
-{\r
-\r
- namespace Solver\r
- {\r
- internal static class CONST\r
- {\r
- public static readonly uint MAX_ARG = 5;\r
- }\r
-\r
- \r
- static public class Constants\r
- {\r
- public static readonly int LIMIT = 5;\r
- public delegate double Func0(double x);\r
- public delegate double Func1(double x, double a);\r
- public delegate double Func2(double x, double a, double b);\r
- public delegate double Func3(double x, double a, double b, double c);\r
- public delegate double Func4(double x, double a, double b, double c, double d);\r
- public delegate double Func5(double x, double a, double b, double c, double d, double e);\r
- }\r
-\r
-\r
- static public class Binomial\r
- {\r
- public static double factorial(int fact)\r
- {\r
- double x = 1;\r
- for (int i = 1; i <= fact; i++)\r
- {\r
- x *= i;\r
- }\r
- return x;\r
- }\r
- public static double combination(int nn, int kk)\r
- {\r
- return (double)(Int64)(factorial(nn) / (factorial(kk) * factorial(nn - kk)));\r
- /*\r
- BigInt n = new BigInt((Int64)nn, new PrecisionSpec());\r
- BigInt k = new BigInt((Int64)kk, new PrecisionSpec());\r
- BigInt n_k = n - k;\r
- n.Factorial();\r
- k.Factorial();\r
- n_k.Factorial();\r
- return (double)(Int64)(n / (k * n_k) );\r
- * */\r
- }\r
- public static double likelihood(double pp, int yes, int no)\r
- {\r
- double p;\r
- if (pp < 0) p = 0; else if (pp > 1) p = 1; else p = pp;\r
- return combination(yes + no, yes) * System.Math.Pow(p, yes) * System.Math.Pow(1 - p, no);\r
- }\r
-\r
- public static double likelihood(double pp, int yes, int no, double comb)\r
- {\r
- double p;\r
- if (pp < 0) p = 0; else if (pp > 1) p = 1; else p = pp;\r
- return comb * System.Math.Pow(p, yes) * System.Math.Pow(1 - p, no);\r
- }\r
-\r
- }\r
-\r
-\r
- public class BernoulliProcess\r
- {\r
- public struct Data\r
- {\r
- public double x;\r
- public int pos, neg, comb;\r
- public void set(double cond, int posit, int negat)\r
- {\r
- x = cond;\r
- pos = posit;\r
- neg = negat;\r
- }\r
-\r
- };\r
- public Data[] elems;\r
- public int length;\r
- public void set(double[][] num)\r
- {\r
- elems = new Data[num.GetLength(0)];\r
- for (int i = 0; i < elems.GetLength(0); i++)\r
- {\r
- elems[i].set(num[i][0], (int)num[i][1], (int)num[i][2]);\r
- }\r
- }\r
- };\r
-\r
-\r
-\r
- class BinomialLikelihoodThread\r
- {\r
- public bool finished;\r
- public static bool started;\r
- internal Thread th;\r
- public static BinomialLikelihoodThread current;\r
-\r
- public Interval[] itvl;\r
- public double[] step;\r
- public double[] champ;\r
- public double champ_like;\r
- public BernoulliProcess data;\r
-\r
-\r
- public Constants.Func0 func0;\r
- public Constants.Func1 func1;\r
- public Constants.Func2 func2;\r
- public Constants.Func3 func3;\r
- public Constants.Func4 func4;\r
- public Constants.Func5 func5;\r
-\r
-\r
- public BinomialLikelihoodThread()\r
- {\r
- itvl = new Interval[CONST.MAX_ARG];\r
- step = new double[CONST.MAX_ARG];\r
- champ = new double[CONST.MAX_ARG];\r
- data = new BernoulliProcess();\r
- }\r
- public void waitLoop()\r
- {\r
- finished = false;\r
- for (int i = 0; i < CONST.MAX_ARG; i++)\r
- {\r
- champ[i] = 0;\r
- }\r
- current = this;\r
- started = false;\r
- }\r
-\r
- public void loop1() { waitLoop(); th = new Thread(new ThreadStart(loop1_)); th.Start(); }\r
- public void loop2() { waitLoop(); th = new Thread(new ThreadStart(loop2_)); th.Start(); }\r
- public void loop3() { waitLoop(); th = new Thread(new ThreadStart(loop3_)); th.Start(); }\r
- public void loop4() { waitLoop(); th = new Thread(new ThreadStart(loop4_)); th.Start(); }\r
- public void loop5() { waitLoop(); th = new Thread(new ThreadStart(loop5_)); th.Start(); }\r
- void loop1_()\r
- {\r
- started = true;\r
- double p;\r
- double like = 1.0;\r
- champ_like=0.0;\r
- int L = data.length;\r
- for (double a = itvl[0].begin.value; a < itvl[0].end.value; a += step[0])\r
- {\r
- like = 1.0;\r
- for(int i=0; i<L; i++) {\r
- p = func1(data.elems[i].x, a);\r
- like *= Binomial.likelihood(p, data.elems[i].pos, data.elems[i].neg, data.elems[i].comb);\r
- }\r
- if(like > champ_like) {\r
- champ_like = like;\r
- champ[0] = a;\r
- }\r
- }\r
- finished = true;\r
- }\r
- void loop2_()\r
- {\r
- started = true;\r
- double p = 0.0;\r
- double like = 1.0;\r
- champ_like = 0.0;\r
- int L = data.length;\r
- for (double a = itvl[0].begin.value; a < itvl[0].end.value; a += step[0])\r
- {\r
- for (double b = itvl[1].begin.value; b < itvl[1].end.value; b += step[1])\r
- {\r
- like = 1.0;\r
- for (int i = 0; i < L; i++)\r
- {\r
- p = func2(data.elems[i].x, a, b);\r
- like *= Binomial.likelihood(p, data.elems[i].pos, data.elems[i].neg, data.elems[i].comb);\r
- }\r
- if (like > champ_like)\r
- {\r
- champ_like = like;\r
- champ[0] = a;\r
- champ[1] = b;\r
- }\r
- }\r
- }\r
- finished = true;\r
- }\r
- void loop3_()\r
- {\r
- started = true;\r
- double p;\r
- double like = 1.0;\r
- champ_like = 0.0;\r
- int L = data.length;\r
- for (double a = itvl[0].begin.value; a < itvl[0].end.value; a += step[0])\r
- {\r
- for (double b = itvl[1].begin.value; b < itvl[1].end.value; b += step[1])\r
- {\r
- for (double c = itvl[2].begin.value; c < itvl[2].end.value; c += step[2])\r
- {\r
- like = 1.0;\r
- for (int i = 0; i < L; i++)\r
- {\r
- p = func3(data.elems[i].x, a, b, c);\r
- like *= Binomial.likelihood(p, data.elems[i].pos, data.elems[i].neg, data.elems[i].comb);\r
- }\r
- if (like > champ_like)\r
- {\r
- champ_like = like;\r
- champ[0] = a;\r
- champ[1] = b;\r
- champ[1] = c;\r
- }\r
- }\r
- }\r
- }\r
- finished = true;\r
- }\r
- void loop4_()\r
- {\r
- started = true;\r
- finished = true;\r
- }\r
- void loop5_()\r
- {\r
- started = true;\r
- finished = true;\r
- }\r
- };\r
-\r
-\r
-\r
- public class BinomialLikelihood\r
- {\r
- /*\r
- public static void showWindow(Constants.Func1 f)\r
- {\r
- Main.canvas.api_canvas.Dispatcher.BeginInvoke(new Action<Constants.Func1>(showWindow_), f);\r
- }\r
- internal static void showWindow_(Constants.Func1 f)\r
- {\r
- System.Windows.Controls.ChildWindow page = new PsychlopsSilverlight4.Pages.BinomialSolver(f);\r
- page.Show();\r
- }\r
- public static void showWindow(Constants.Func2 f)\r
- {\r
- Main.canvas.api_canvas.Dispatcher.BeginInvoke(new Action<Constants.Func2>(showWindow_), f);\r
- }\r
- internal static void showWindow_(Constants.Func2 f)\r
- {\r
- System.Windows.Controls.ChildWindow page = new PsychlopsSilverlight4.Pages.BinomialSolver(f);\r
- page.Show();\r
- }\r
- public static void showWindow(Constants.Func3 f)\r
- {\r
- Main.canvas.api_canvas.Dispatcher.BeginInvoke(new Action<Constants.Func3>(showWindow_), f);\r
- }\r
- internal static void showWindow_(Constants.Func3 f)\r
- {\r
- System.Windows.Controls.ChildWindow page = new PsychlopsSilverlight4.Pages.BinomialSolver(f);\r
- page.Show();\r
- }\r
- */\r
-\r
- public int iteration;\r
-\r
- public Interval[] itvl;\r
- public double[] step;\r
- public double[] champ;\r
- public double champ_like;\r
- public BernoulliProcess data;\r
-\r
- public Constants.Func0 func0;\r
- public Constants.Func1 func1;\r
- public Constants.Func2 func2;\r
- public Constants.Func3 func3;\r
- public Constants.Func4 func4;\r
- public Constants.Func5 func5;\r
-\r
- public BinomialLikelihood()\r
- {\r
- itvl = new Interval[CONST.MAX_ARG];\r
- step = new double[CONST.MAX_ARG];\r
- champ = new double[CONST.MAX_ARG];\r
- iteration = 2;\r
- }\r
-\r
- public void begin(Constants.Func0 d_func)\r
- { }\r
- public void begin(Constants.Func1 d_func)\r
- {\r
- func1 = d_func;\r
-\r
- BinomialLikelihoodThread[] l = new BinomialLikelihoodThread[4];\r
- for (int i = 0; i < 4; i++) { l[i] = new BinomialLikelihoodThread(); }\r
-\r
- for (int k = 0; k < iteration; k++)\r
- {\r
- begin_base(l);\r
- for (int i = 0; i < 4; i++)\r
- {\r
- l[i].data = data;\r
- l[i].func1 = d_func;\r
- l[i].loop1();\r
- }\r
- end_base(l);\r
- }\r
- }\r
- public void begin(Constants.Func2 d_func)\r
- {\r
- func2 = d_func;\r
-\r
- BinomialLikelihoodThread[] l = new BinomialLikelihoodThread[4];\r
- for (int i = 0; i < 4; i++) { l[i] = new BinomialLikelihoodThread(); }\r
-\r
- for (int k = 0; k < iteration; k++)\r
- {\r
- begin_base(l);\r
- for (int i = 0; i < 4; i++)\r
- {\r
- l[i].data = data;\r
- l[i].func2 = d_func;\r
- l[i].loop2();\r
- }\r
- end_base(l);\r
- }\r
- }\r
- public void begin(Constants.Func3 d_func)\r
- {\r
- func3 = d_func;\r
-\r
- BinomialLikelihoodThread[] l = new BinomialLikelihoodThread[4];\r
- for (int i = 0; i < 4; i++) { l[i] = new BinomialLikelihoodThread(); }\r
-\r
- for (int k = 0; k < iteration; k++)\r
- {\r
- begin_base(l);\r
- for(int i=0; i<4; i++) {\r
- l[i].data = data;\r
- l[i].func3 = d_func;\r
- l[i].loop3();\r
- }\r
- end_base(l);\r
- }\r
- }\r
-\r
- public void begin(Constants.Func4 d_func)\r
- {\r
- }\r
- public void begin(Constants.Func5 d_func)\r
- {\r
- }\r
-\r
- void begin_base(BinomialLikelihoodThread[] l)\r
- {\r
- champ_like = 0;\r
-\r
- data.length = data.elems.GetLength(0);\r
- for (int i = 0; i < data.elems.GetLength(0); i++)\r
- {\r
- data.elems[i].comb = (int)Binomial.combination(data.elems[i].pos + data.elems[i].neg, data.elems[i].pos);\r
- }\r
-\r
- for (int j = 0; j < Constants.LIMIT; j++) { step[j] = (itvl[j].end.value - itvl[j].begin.value) / 256.0; }\r
- for (int i = 0; i < 4; i++)\r
- {\r
- l[i].itvl[0] = new Interval((itvl[0].begin.value) + (i * step[0] * 64), (itvl[0].begin.value) + ((i + 1) * step[0] * 64));\r
- l[i].step[0] = step[0];\r
- for (int j = 1; j < Constants.LIMIT; j++)\r
- {\r
- l[i].itvl[j] = itvl[j];\r
- l[i].step[j] = step[j];\r
- }\r
- l[i].data = data;\r
- }\r
- }\r
-\r
- void end_base(BinomialLikelihoodThread[] l)\r
- {\r
- //for (int i = 0; i < 4; i++) { l[i].th.Join(); }\r
- while (!l[0].finished || !l[1].finished || !l[2].finished || !l[3].finished) { Thread.Sleep(100); } \r
-\r
- for(int j=0; j<Constants.LIMIT; j++) { champ[j] = 0; }\r
- champ_like = 0.0;\r
- for(int i=0; i<4; i++) {\r
- if(champ_like < l[i].champ_like) {\r
- champ_like = l[i].champ_like;\r
- for(int j=0; j<Constants.LIMIT; j++) { champ[j] = l[i].champ[j]; }\r
- }\r
- }\r
-\r
- double r, low, high;\r
- for (int j = 0; j < Constants.LIMIT; j++)\r
- {\r
- r = itvl[j].end.value - itvl[j].begin.value;\r
- low = champ[j] - r / 8.0 < itvl[j].begin.value ? itvl[j].begin.value : champ[j] - r / 8.0;\r
- high = champ[j] + r / 8.0 > itvl[j].end.value ? itvl[j].end.value : champ[j] + r / 8.0;\r
- itvl[j] = new Interval(low, high);\r
- }\r
-\r
- }\r
-\r
-\r
- }\r
-\r
-\r
- }\r
-\r
-}
\ No newline at end of file