X-Git-Url: http://git.osdn.jp/view?a=blobdiff_plain;f=dev4%2Fpsychlops%2Fextention%2Fmath%2Fsolver.cs;fp=dev4%2Fpsychlops%2Fextention%2Fmath%2Fsolver.cs;h=e26365fcd9d4a5d5f457a23e1313471d144bf30c;hb=397aea915b85fc9457e2131db7ca78d8269b30b0;hp=0000000000000000000000000000000000000000;hpb=a18f323a3bf447b6f2ce12d0d1a4f39a07814e93;p=psychlops%2Fsilverlight.git diff --git a/dev4/psychlops/extention/math/solver.cs b/dev4/psychlops/extention/math/solver.cs new file mode 100644 index 0000000..e26365f --- /dev/null +++ b/dev4/psychlops/extention/math/solver.cs @@ -0,0 +1,384 @@ +using System.Threading; +using System; +using BigNum; + +namespace Psychlops +{ + + namespace Solver + { + + + static public class Constants + { + public static readonly int LIMIT = 5; + public delegate double Func0(double x); + public delegate double Func1(double x, double a); + public delegate double Func2(double x, double a, double b); + public delegate double Func3(double x, double a, double b, double c); + public delegate double Func4(double x, double a, double b, double c, double d); + public delegate double Func5(double x, double a, double b, double c, double d, double e); + } + + + static public class Binomial + { + public static double factorial(int fact) + { + double x = 1; + for (int i = 1; i <= fact; i++) + { + x *= i; + } + return x; + } + public static double combination(int nn, int kk) + { + //return (double)(Int64)(factorial(n) / (factorial(kk) * factorial(n - kk))); + BigInt n = new BigInt((Int64)nn, new PrecisionSpec()); + BigInt k = new BigInt((Int64)kk, new PrecisionSpec()); + BigInt n_k = n - k; + n.Factorial(); + k.Factorial(); + n_k.Factorial(); + return (double)(Int64)(n / (k * n_k) ); + } + public static double likelihood(double pp, int yes, int no) + { + double p; + if (pp < 0) p = 0; else if (pp > 1) p = 1; else p = pp; + return combination(yes + no, yes) * System.Math.Pow(p, yes) * System.Math.Pow(1 - p, no); + } + + public static double likelihood(double pp, int yes, int no, double comb) + { + double p; + if (pp < 0) p = 0; else if (pp > 1) p = 1; else p = pp; + return comb * System.Math.Pow(p, yes) * System.Math.Pow(1 - p, no); + } + + } + + + public class BernoulliProcess + { + public struct Data + { + public double x; + public int pos, neg, comb; + public void set(double cond, int posit, int negat) + { + x = cond; + pos = posit; + neg = negat; + } + + }; + public Data[] elems; + public int length; + public void set(double[][] num) + { + elems = new Data[num.GetLength(0)]; + for (int i = 0; i < elems.GetLength(0); i++) + { + elems[i].set(num[i][0], (int)num[i][1], (int)num[i][2]); + } + } + }; + + + + class BinomialLikelihoodThread + { + public bool finished; + public static bool started; + internal Thread th; + public static BinomialLikelihoodThread current; + + public Interval[] itvl; + public double[] step; + public double[] champ; + public double champ_like; + public BernoulliProcess data; + + + public Constants.Func0 func0; + public Constants.Func1 func1; + public Constants.Func2 func2; + public Constants.Func3 func3; + public Constants.Func4 func4; + public Constants.Func5 func5; + + + public BinomialLikelihoodThread() + { + } + public void waitLoop() + { + finished = false; + for(int i=0; i<10; i++) { + champ[i] = 0; + } + current = this; + started = false; + } + + public void loop1() { waitLoop(); th = new Thread(new ThreadStart(loop1_)); } + public void loop2() { waitLoop(); th = new Thread(new ThreadStart(loop2_)); } + public void loop3() { waitLoop(); th = new Thread(new ThreadStart(loop3_)); } + public void loop4() { waitLoop(); th = new Thread(new ThreadStart(loop4_)); } + public void loop5() { waitLoop(); th = new Thread(new ThreadStart(loop5_)); } + void loop1_() + { + started = true; + double p; + double like = 1.0; + champ_like=0.0; + int L = data.length; + for (double a = itvl[0].begin.val; a < itvl[0].end.val; a += step[0]) + { + like = 1.0; + for(int i=0; i champ_like) { + champ_like = like; + champ[0] = a; + } + } + finished = true; + } + void loop2_() + { + started = true; + double p; + double like = 1.0; + champ_like = 0.0; + int L = data.length; + for (double a = itvl[0].begin.val; a < itvl[0].end.val; a += step[0]) + { + for (double b = itvl[1].begin.val; b < itvl[1].end.val; b += step[1]) + { + like = 1.0; + for (int i = 0; i < L; i++) + { + p = func2(data.elems[i].x, a, b); + like *= Binomial.likelihood(p, data.elems[i].pos, data.elems[i].neg, data.elems[i].comb); + } + if (like > champ_like) + { + champ_like = like; + champ[0] = a; + } + } + } + finished = true; + } + void loop3_() + { + started = true; + double p; + double like = 1.0; + champ_like = 0.0; + int L = data.length; + for (double a = itvl[0].begin.val; a < itvl[0].end.val; a += step[0]) + { + for (double b = itvl[1].begin.val; b < itvl[1].end.val; b += step[1]) + { + for (double c = itvl[2].begin.val; c < itvl[2].end.val; c += step[2]) + { + like = 1.0; + for (int i = 0; i < L; i++) + { + p = func3(data.elems[i].x, a, b, c); + like *= Binomial.likelihood(p, data.elems[i].pos, data.elems[i].neg, data.elems[i].comb); + } + if (like > champ_like) + { + champ_like = like; + champ[0] = a; + } + } + } + } + finished = true; + } + void loop4_() + { + started = true; + finished = true; + } + void loop5_() + { + started = true; + finished = true; + } + }; + + + + class BinomialLikelihood + { + public int iteration; + + public Interval[] itvl; + public double[] step; + public double[] champ; + public double champ_like; + public BernoulliProcess data; + + public Constants.Func0 func0; + public Constants.Func1 func1; + public Constants.Func2 func2; + public Constants.Func3 func3; + public Constants.Func4 func4; + public Constants.Func5 func5; + + public BinomialLikelihood() + { + iteration = 2; + } + + public void begin(Constants.Func0 d_func) + { } + public void begin(Constants.Func1 d_func) + { + func1 = d_func; + + BinomialLikelihoodThread[] l = new BinomialLikelihoodThread[4]; + + double r, low, high; + for (int k = 0; k < iteration; k++) + { + champ_like = 0; + + begin_base(l); + for (int i = 0; i < 4; i++) + { + l[i].data = data; + l[i].func1 = d_func; + l[i].loop1(); + } + end_base(l); + + for (int j = 0; j < Constants.LIMIT; j++) + { + r = itvl[j].end.val - itvl[j].begin.val; + low = champ[j] - r / 8.0 < itvl[j].begin.val ? itvl[j].begin.val : champ[j] - r / 8.0; + high = champ[j] + r / 8.0 > itvl[j].end.val ? itvl[j].end.val : champ[j] + r / 8.0; + itvl[j] = new Interval(low, high); + } + } + } + public void begin(Constants.Func2 d_func) + { + func2 = d_func; + + BinomialLikelihoodThread[] l = new BinomialLikelihoodThread[4]; + + double r, low, high; + for (int k = 0; k < iteration; k++) + { + champ_like = 0; + + begin_base(l); + for (int i = 0; i < 4; i++) + { + l[i].data = data; + l[i].func2 = d_func; + l[i].loop2(); + } + end_base(l); + + for (int j = 0; j < Constants.LIMIT; j++) + { + r = itvl[j].end.val - itvl[j].begin.val; + low = champ[j] - r / 8.0 < itvl[j].begin.val ? itvl[j].begin.val : champ[j] - r / 8.0; + high = champ[j] + r / 8.0 > itvl[j].end.val ? itvl[j].end.val : champ[j] + r / 8.0; + itvl[j] = new Interval(low, high); + } + } + } + public void begin(Constants.Func3 d_func) + { + func3 = d_func; + + BinomialLikelihoodThread[] l = new BinomialLikelihoodThread[4]; + + double r, low, high; + for (int k = 0; k < iteration; k++) + { + champ_like = 0; + + begin_base(l); + for(int i=0; i<4; i++) { + l[i].data = data; + l[i].func3 = d_func; + l[i].loop3(); + } + end_base(l); + + for (int j = 0; j < Constants.LIMIT; j++) + { + r = itvl[j].end.val - itvl[j].begin.val; + low = champ[j]-r/8.0 < itvl[j].begin.val ? itvl[j].begin.val : champ[j]-r/8.0; + high = champ[j]+r/8.0 > itvl[j].end.val ? itvl[j].end.val : champ[j]+r/8.0; + itvl[j] = new Interval(low, high); + } + } + } + + public void begin(Constants.Func4 d_func) + { + } + public void begin(Constants.Func5 d_func) + { + } + + void begin_base(BinomialLikelihoodThread[] l) + { + data.length = data.elems.GetLength(0); + for (int i = 0; i < data.elems.GetLength(0); i++) + { + data.elems[i].comb = (int)Binomial.combination(data.elems[i].pos + data.elems[i].neg, data.elems[i].pos); + } + + for (int j = 0; j < Constants.LIMIT; j++) { step[j] = (itvl[j].end.val - itvl[j].begin.val) / 256.0; } + for (int i = 0; i < 4; i++) + { + l[i].itvl[0] = new Interval((itvl[0].begin.val) + (i * step[0] * 64), (itvl[0].begin.val) + ((i + 1) * step[0] * 64)); + l[i].step[0] = step[0]; + for (int j = 1; j < Constants.LIMIT; j++) + { + l[i].itvl[j] = itvl[j]; + l[i].step[j] = step[j]; + } + l[i].data = data; + } + } + + void end_base(BinomialLikelihoodThread[] l) + { + for (int i = 0; i < 4; i++) + { + l[i].th.Join(); + } + + for(int j=0; j