using System;
using System.Collections.Generic;
using System.ComponentModel;
using System.Data;
using System.Drawing;
using System.Drawing.Imaging;
using System.Linq;
using System.IO;
using System.Threading;
using System.Text;
using System.Runtime.InteropServices;
using System.Windows.Forms;
using System.Management;
using System.Diagnostics;
using System.Threading.Tasks;

namespace Fst
{
    public partial class Form1 : Form
    {
        #region DLL

        [DllImport("user32.dll", EntryPoint = "GetWindowText")]
        public static extern int GetWindowText(IntPtr hwnd, StringBuilder text, int len);
        #endregion

        #region STATISTICS
        private class SORTFDR
        {
            public double p;
            public int id;
            public int rank;
            public SORTFDR(double p1, int id1)
            {
                p = p1;
                id = id1;
            }
        }

        public static void Add(double[] a, double[] b, double[] c)
        {
            int m = c.GetLength(0);
            for (int i = 0; i < m; ++i)
                a[i] = b[i] + c[i];
        }

        public static void Add(double[,] a, double[,] b, double[,] c)
        {
            int m = c.GetLength(0);
            int n = c.GetLength(1);
            for (int i = 0; i < m; ++i)
                for (int j = 0; j < n; ++j)
                    a[i, j] = b[i, j] + c[i, j];
        }
        
        public static void Mul(double[] a, double[] b, double val)
        {
            int n = a.GetLength(0);
            for (int i = 0; i < n; ++i)
                a[i] = b[i] * val;
        }

        public static void Mul(double[,] a, double[,] b, double val)
        {
            int m = a.GetLength(0);
            int n = a.GetLength(1);
            for (int i = 0; i < m; ++i)
                for (int j = 0; j < n; ++j)
                    a[i, j] = b[i, j] * val;
        }

        public static void SetVal<T>(T[,] a, T val)
        {
            int m = a.GetLength(0);
            int n = a.GetLength(1);
            for (int i = 0; i < m; ++i)
                for (int j = 0; j < n; ++j)
                    a[i, j] = val;
        }

        public static void QuickSort(int[] a, int left, int right)
        {
            if (left > right) return;
            int i = left, j = right, temp = a[left];
            while (i != j)
            {
                while (a[j] >= temp && j > i) j--;
                if (j > i) a[i++] = a[j];

                while (a[i] <= temp && j > i) i++;
                if (j > i) a[j--] = a[i];
            }
            a[i] = temp;
            QuickSort(a, left, i - 1);
            QuickSort(a, i + 1, right);
        }

        public static void GetRandSeq(int[] td, Random rnd)
        {
            int nn = td.Length;
            for (int i = 0; i < nn; ++i)
                td[i] = ((rnd.Next() << 16) | i);
            QuickSort(td, 0, nn - 1);
            for (int i = 0; i < nn; ++i)
                td[i] &= 0xFFFF;
        }

        public static void Swap<T>(ref T val1, ref T val2)
        {
            T t = val1;
            val1 = val2;
            val2 = t;
        }

        public static void Permute<T>(T[] val, Random rnd, int n = -1)
        {
            for (int i = n == -1 ? val.Length : n; i > 1; --i)
                Swap(ref val[i - 1], ref val[rnd.Next(i)]);
        }

        public static double[,] BINOMIAL;//small big

        public static double Factorial(int n)
        {
            switch (n)
            {
                case 0: return 1;
                case 1: return 1;
                case 2: return 2;
                case 3: return 6;
                case 4: return 24;
                case 5: return 120;
                case 6: return 720;
                case 7: return 5040;
                case 8: return 40320;
                case 9: return 362880;
                case 10: return 3628800;
                case 11: return 39916800;
                case 12: return 479001600;
                case 13: return 6227020800;
                case 14: return 87178291200;
                case 15: return 1307674368000;
                case 16: return 20922789888000;
                case 17: return 355687428096000;
                case 18: return 6402373705728000;
                case 19: return 121645100408832000;
                case 20: return 2432902008176640000;
                case 21: return 51090942171709440000.0;
                case 22: return 1124000727777607680000.0;
                case 23: return 25852016738884976640000.0;
                case 24: return 620448401733239439360000.0;
                default:
                    double re = 620448401733239439360000.0;
                    for (int i = 25; i <= n; ++i)
                        re *= i;
                    return re;
            }
        }

        public static double Binomial(int n, int r)
        {
            return Factorial(n) / Factorial(r) / Factorial(n - r);
        }

        public static double[] GAMMA_P = new double[] { 0.99999999999980993, 676.5203681218851, -1259.1392167224028, 771.32342877765313, -176.61502916214059, 12.507343278686905, -0.13857109526572012, 9.9843695780195716e-6, 1.5056327351493116e-7 };

        public static double[] GAMMA_P2 = new double[] { 1.000000000190015, 76.18009172947146, -86.50532032941677, 24.01409824083091, -1.231739572450155, 0.0012086509738662, -0.000005395239385 };

        public static double Gamma(double x)
        {
            int g = 7;
            if (x < 0.5) return Math.PI / (Math.Sin(Math.PI * x) * Gamma(1 - x));
            x -= 1;
            double a = GAMMA_P[0];
            double t = x + g + 0.5;
            for (int i = 1; i < 9; ++i)
                a += GAMMA_P[i] / (x + i);
            return 2.5066282746310005 * Math.Pow(t, x + 0.5) * Math.Exp(-t) * a;
        }

        public static double GammaLn(double x)
        {
            if (x < 3.72008E-44) return GammaLn(3.72008E-44);
            double y = x, ser = GAMMA_P2[0];
            for (int i = 1; i < 7; ++i)
                ser += (GAMMA_P2[i] / ++y);
            return Math.Log(2.5066282746310005 * ser / x) - (x + 5.5 - (x + 0.5) * Math.Log(x + 5.5));
        }

        public static double GetRandNorm(Random rnd)
        {
            double U1 = Math.Sqrt(-2 * Math.Log(rnd.NextDouble()));
            double U2 = 2 * Math.PI * rnd.NextDouble();
            return U1 * Math.Cos(U2);
        }
        
        public static double GetRandGam(Random rnd, double alpha, double beta = 1)
        {
            if (alpha < 1) return GetRandGam(rnd, 1.0 + alpha, beta) * Math.Pow(rnd.NextDouble(), 1.0 / alpha);
            double x, v, u;
            double d = alpha - 1.0 / 3.0;
            double c = (1.0 / 3.0) / Math.Sqrt(d);

            while (true)
            {
                do
                {
                    x = GetRandNorm(rnd); ;
                    v = 1.0 + c * x;
                }
                while (v <= 0);

                v = v * v * v;
                u = rnd.NextDouble();

                if (u < 1 - 0.0331 * x * x * x * x) break;
                if (Math.Log(u) < 0.5 * x * x + d * (1 - v + Math.Log(v))) break;
            }
            return beta * d * v;
        }
        
        public static void GetRandDirichlet(Random rnd, double[,] q, double[,] p, double F)
        {
            if (F == 0)
            {
                Array.Copy(p, q, p.Length);
                return;
            }

            double f = 1.0 / F - 1;
            if (f < 0.01) f = 0.01;
            if (f > 1e30) f = 1e30;

            int L = q.GetLength(0), K = q.GetLength(1);
            Mul(q, p, f);
            for (int l = 0; l < L; ++l)
            {
                double sum = 0;
                for (int k = 0; k < K; ++k)
                {
                    q[l, k] = GetRandGam(rnd, q[l, k]) + 1e-300;
                    sum += q[l, k];
                }
                for (int k = 0; k < K; ++k)
                    q[l, k] /= sum;
            }
        }
               
        #endregion

        #region Simulation

        //v	L	k	C	Sc	Fis	Fst	Nf	Nm	No	nrep
        public static int v = 2;//ploidy
        public static int L = 32; //# ssr
        public static int K = 4; //# alleles at ssr
        public static int C = 1, Sc = 1;
        public static int S = 1; //# pops
        public static double Fis = 0.5; //diff among inds at ssr
        public static double Fst = 0.0; //diff among pops at ssr
        public static double Fct = 0.0; //diff among pops at ssr
        public static int NF = 1, NM = 1, NO = 1;
        public static int Ns = 0; //# inds in each pop
        public static int Nt = 0; //# inds in each pop
        public static int Nrep = 0;
        public static double nrfreq = 0.75;

        public static Dictionary<ulong, Genotype> genotypes;
        public static Genotype[] genotypeid;

        public void PrepareGenotypes()
        {
            int[] allele = new int[v];
            int ngssr = (int)Binomial(v + K - 1, v);
            genotypes = new Dictionary<ulong, Genotype>();
            genotypeid = new Genotype[ngssr];

            if (v == 2)
            {
                if (L > 0) for (int a1 = 0, id = 0; a1 < K; ++a1)
                {
                    allele[0] = a1;
                    for (int a2 = a1; a2 < K; ++a2, ++id)
                    {
                        allele[1] = a2;
                        Genotype t = new Genotype(allele, id);
                        genotypes[t.hash] = t;
                        genotypeid[id] = t;
                    }
                }
            }
            if (v == 4)
            {
                if (L > 0) for (int a1 = 0, id = 0; a1 < K; ++a1)
                    {
                        allele[0] = a1;
                        for (int a2 = a1; a2 < K; ++a2)
                        {
                            allele[1] = a2;
                            for (int a3 = a2; a3 < K; ++a3)
                            {
                                allele[2] = a3;
                                for (int a4 = a3; a4 < K; ++a4, ++id)
                                {
                                    allele[3] = a4;
                                    Genotype t = new Genotype(allele, id);
                                    genotypes[t.hash] = t;
                                    genotypeid[id] = t;
                                }
                            }
                        }
                    }
            }
            if (v == 6)
            {
                if (L > 0) for (int a1 = 0, id = 0; a1 < K; ++a1)
                    {
                        allele[0] = a1;
                        for (int a2 = a1; a2 < K; ++a2)
                        {
                            allele[1] = a2;
                            for (int a3 = a2; a3 < K; ++a3)
                            {
                                allele[2] = a3;
                                for (int a4 = a3; a4 < K; ++a4)
                                {
                                    allele[3] = a4;
                                    for (int a5 = a4; a5 < K; ++a5)
                                    {
                                        allele[4] = a5;
                                        for (int a6 = a5; a6 < K; ++a6, ++id)
                                        {
                                            allele[5] = a6;
                                            Genotype t = new Genotype(allele, id);
                                            genotypes[t.hash] = t;
                                            genotypeid[id] = t;
                                        }
                                    }
                                }
                            }
                        }
                    }
            }
        }

        static Dictionary<ulong, double> pxy = new Dictionary<ulong, double>();

        public class Genotype
        {
            public int[] allele;
            Dictionary<int, int> allelecount;
            public double hwecoef = -1;
            public ulong hash;
            public int allelepattern = 0;
            public int id;
            public double ho = -1, px2 = -1, sswi = -1;

            public Genotype(int[] _allele, int _id)
            {
                allele = new int[v];
                Array.Copy(_allele, allele, v);
                GetHweCoef();
                id = _id;
            }

            public Genotype(Random rnd, Genotype f, Genotype m, int[] buf)
            {
                allele = new int[v];
                GetRandSeq(buf, rnd);
                for (int i = 0; i < v >> 1; ++i)
                    allele[i           ] = f.allele[buf[i]];
                GetRandSeq(buf, rnd);
                for (int i = 0; i < v >> 1; ++i)
                    allele[(v >> 1) + i] = m.allele[buf[i]];
                GetHweCoef();
            }

            public double Pxy(Genotype y)
            {
                ulong ha = 0;
                if (hash < y.hash)
                    ha = ((ulong)y.hash << 32) | (ulong)hash;
                else
                    ha = ((ulong)hash << 32) | (ulong)y.hash;
                if (!pxy.ContainsKey(ha))
                {
                    double re = 0;
                    foreach (int k in allelecount.Keys)
                    {
                        if (y.allelecount.ContainsKey(k))
                            re += allelecount[k] * y.allelecount[k];
                    }
                    pxy[ha] = re / (double)(v * v); 
                }
                return pxy[ha];
            }

            public double SSwi()
            {
                if (sswi < 0)
                {
                    sswi = 0;
                    for (int i = 0; i < v; ++i)
                        for (int j = i + 1; j < v; ++j)
                            if (allele[i] != allele[j])
                                sswi++;
                    sswi /= v;
                }
                return sswi;
            }

            public double GetHo()
            {
                if (ho < 0)
                {
                    double re = 0;
                    for (int i = 0; i < v; ++i)
                        for (int j = i + 1; j < v; ++j)
                            if (allele[i] != allele[j])
                                re++;
                    ho = re / (v * (v - 1) / 2);
                }
                return ho;
            }

            public double GetPx2()
            {
                if (px2 < 0)
                {
                    px2 = 0;
                    foreach(int k in allelecount.Keys)
                    {
                        double f = allelecount[k] / (double)v;
                        px2 += f * f;
                    }
                }
                return px2;
            }

            public double GetHweCoef()
            {
                if (hwecoef >= 0) return hwecoef;

                allelepattern = 0;
                allelecount = new Dictionary<int, int>();
                foreach (int a in allele)
                {
                    if (!allelecount.ContainsKey(a)) allelecount[a] = 0;
                    allelecount[a]++;
                }
                for (int a = 0; a < K; ++a)
                    if (!allelecount.ContainsKey(a))
                        hash <<= 4;
                    else
                        hash = (hash << 4) | (ulong)allelecount[a];

                if (genotypes.ContainsKey(hash))
                    id = genotypes[hash].id;

                allelepattern = 0;

                foreach (var a in from o in allelecount orderby o.Value descending select o)
                    allelepattern = (allelepattern << 4) | a.Value;

                switch (allelepattern)
                {
                    default: hwecoef = 1; break;
                    case 0x1: hwecoef = 1.0; break;
                    case 0x2: hwecoef = 1.0; break;
                    case 0x11: hwecoef = 2.0; break;
                    case 0x3: hwecoef = 1.0; break;
                    case 0x21: hwecoef = 3.0; break;
                    case 0x111: hwecoef = 6.0; break;
                    case 0x4: hwecoef = 1.0; break;
                    case 0x31: hwecoef = 4.0; break;
                    case 0x22: hwecoef = 6.0; break;
                    case 0x211: hwecoef = 12.0; break;
                    case 0x1111: hwecoef = 24.0; break;
                    case 0x5: hwecoef = 1.0; break;
                    case 0x41: hwecoef = 5.0; break;
                    case 0x32: hwecoef = 10.0; break;
                    case 0x311: hwecoef = 20.0; break;
                    case 0x221: hwecoef = 30.0; break;
                    case 0x2111: hwecoef = 60.0; break;
                    case 0x11111: hwecoef = 120.0; break;
                    case 0x6: hwecoef = 1.0; break;
                    case 0x51: hwecoef = 6.0; break;
                    case 0x42: hwecoef = 15.0; break;
                    case 0x411: hwecoef = 30.0; break;
                    case 0x33: hwecoef = 20.0; break;
                    case 0x321: hwecoef = 60.0; break;
                    case 0x3111: hwecoef = 120.0; break;
                    case 0x222: hwecoef = 90.0; break;
                    case 0x2211: hwecoef = 180.0; break;
                    case 0x21111: hwecoef = 360.0; break;
                    case 0x111111: hwecoef = 720.0; break;
                    case 0x7: hwecoef = 1.0; break;
                    case 0x61: hwecoef = 7.0; break;
                    case 0x52: hwecoef = 21.0; break;
                    case 0x511: hwecoef = 42.0; break;
                    case 0x43: hwecoef = 35.0; break;
                    case 0x421: hwecoef = 105.0; break;
                    case 0x4111: hwecoef = 210.0; break;
                    case 0x331: hwecoef = 140.0; break;
                    case 0x322: hwecoef = 210.0; break;
                    case 0x3211: hwecoef = 420.0; break;
                    case 0x31111: hwecoef = 840.0; break;
                    case 0x2221: hwecoef = 630.0; break;
                    case 0x22111: hwecoef = 1260.0; break;
                    case 0x211111: hwecoef = 2520.0; break;
                    case 0x1111111: hwecoef = 5040.0; break;
                    case 0x8: hwecoef = 1.0; break;
                    case 0x71: hwecoef = 8.0; break;
                    case 0x62: hwecoef = 28.0; break;
                    case 0x611: hwecoef = 56.0; break;
                    case 0x53: hwecoef = 56.0; break;
                    case 0x521: hwecoef = 168.0; break;
                    case 0x5111: hwecoef = 336.0; break;
                    case 0x44: hwecoef = 70.0; break;
                    case 0x431: hwecoef = 280.0; break;
                    case 0x422: hwecoef = 420.0; break;
                    case 0x4211: hwecoef = 840.0; break;
                    case 0x41111: hwecoef = 1680.0; break;
                    case 0x332: hwecoef = 560.0; break;
                    case 0x3311: hwecoef = 1120.0; break;
                    case 0x3221: hwecoef = 1680.0; break;
                    case 0x32111: hwecoef = 3360.0; break;
                    case 0x311111: hwecoef = 6720.0; break;
                    case 0x2222: hwecoef = 2520.0; break;
                    case 0x22211: hwecoef = 5040.0; break;
                    case 0x221111: hwecoef = 10080.0; break;
                    case 0x2111111: hwecoef = 20160.0; break;
                    case 0x11111111: hwecoef = 40320.0; break;
                }
                return hwecoef;
            }

            private double CondDeltaSub(Dictionary<int, int> x, Dictionary<int, int> y, double coef, int k, int rd, double[,] fre, int l, double g)
            {
                int xk = x.ContainsKey(k) ? x[k] : 0, yk = y.ContainsKey(k) ? y[k] : 0;
                int mink = Math.Min(xk, yk); mink = Math.Min(rd, mink);

                double re = 0;
                if (k == K - 1)
                {
                    if (rd <= mink)
                    {
                        int zk = rd;
                        double gk = fre[l, k] * g;
                        re = coef / (Factorial(xk - zk) * Factorial(zk) * Factorial(yk - zk));
                        for (int j = 0; j < xk + yk - zk; ++j)
                            re *= gk + j;
                    }
                }
                else for (int zk = 0; zk <= mink; ++zk)
                {
                    double c2 = 1.0 / (Factorial(xk - zk) * Factorial(zk) * Factorial(yk - zk));
                    double gk = fre[l, k] * g;
                    for (int j = 0; j < xk + yk - zk; ++j)
                        c2 *= gk + j;
                    re += CondDeltaSub(x, y, coef * c2, k + 1, rd - zk, fre, l, g);
                }
                return re;
            }

            public double CondDelta(Genotype y, int d, double[,] fre, int l, double f)
            {
                if (allelecount.Count == 0) return 0;
                if (f <= 1e-6) f = 1e-6;
                if (f >= 0.999999) f = 0.999999;

                double g = 1.0 / f - 1;
                double coef = 1.0;
                for (int j = 0; j < v + v - d; ++j)
                    coef *= g + j;
                coef = Factorial(v - d) * Factorial(d) * Factorial(v - d) / coef;

                return CondDeltaSub(allelecount, y.allelecount, coef, 0, d, fre, l, g);
            }

            public double CondMother(Genotype m, double[,] fre, int l, double f)
            {
                if (allelecount.Count == 0 || m.allelecount.Count == 0) return 0;
                if (f <= 1e-6) f = 1e-6;
                if (f >= 0.999999) f = 0.999999;

                int K = fre.GetLength(1);
                double re = m.hwecoef, re2 = 1, g = 1.0 / f - 1;
                for (int k = 0; k < K; ++k)
                {
                    int nfk = allelecount.ContainsKey(k) ? allelecount[k] : 0;
                    int nc = nfk + (m.allelecount.ContainsKey(k) ? m.allelecount[k] : 0);
                    double gk = fre[l, k] * g;
                    for (int j = nfk; j < nc; ++j)
                        re *= gk + j;
                }
                for (int j = v; j < v + v; ++j)
                    re2 *= g + j;
                return re / re2;
            }

            public double GFZ(double[,] fre, int l, double f)
            {
                if (allelecount.Count == 0) return 0;
                if (f <= 1e-6) f = 1e-6;
                if (f >= 0.999999) f = 0.999999;

                double re = hwecoef;
                double re2 = 1, g = 1.0 / f - 1;
                foreach (var b in allelecount)
                {
                    double gk = fre[l, b.Key] * g;
                    for (int j = 0; j < b.Value; ++j)
                        re *= gk + j;
                }
                for (int j = 0; j < v; ++j)
                    re2 *= g + j;
                return re / re2;
            }
        };

        public class Individual
        {
            public int id;
            public int fid;
            public int mid;
            public int popid;
            public Genotype[] ssr;

            public Individual(Random rnd, int _id, Population pop)
            {
                id = _id;
                popid = pop.id;
                ssr = new Genotype[L];
                for (int l = 0; l < L; ++l)
                    ssr[l] = pop.DrawGenotype(rnd, l);
            }

            public Individual(Random rnd, int _id, Individual f, Population pop)
            {
                id = _id;
                popid = pop.id;
                ssr = new Genotype[L];
                for (int l = 0; l < L; ++l)
                    ssr[l] = pop.DrawCondGenotype(rnd, f.ssr[l], l);
            }

            public Individual(Random rnd, int _id, Population pop, Individual f, Individual m, int[] buf)
            {
                id = _id;
                popid = pop.id;
                fid = f.id;
                mid = m.id;
                ssr = new Genotype[L];

                for (int l = 0; l < L; ++l)
                    ssr[l] = genotypeid[new Genotype(rnd, f.ssr[l], m.ssr[l], buf).id];
            }

        };

        public class Population
        {
            public int id;
            public int ninds;
            public Population[] subpop;
            public Individual[] inds;

            public double[,] f;
            public double[,] tf;
            public double[,] f_obs;
            public double[] hs_obs;
            public double[] ht_obs;
            public double[] hi_obs;
            public double[] px2_obs;
            public double[] ps2_obs;
            public double[] pt2_obs;
            public double[] msl;
            public double[] mtl;
            public double[][,] cond; //L*Ng*Ng
            public double[][] gfz;

            public double[] SSwi, SSwp, SStot;

            public Population(Random rnd, bool preparecond, bool isregion, bool randfreq)
            {
                //root
                id        = -1;
                ninds     = Ns * S;

                inds      = new Individual[ninds * S];
                subpop    = new Population[S];

                f     = new double[L, K];
                tf    = new double[L, K];
                f_obs = new double[L, K];

                px2_obs = new double[L];
                ps2_obs = new double[L];
                pt2_obs = new double[L];

                hi_obs = new double[L];
                hs_obs = new double[L];
                ht_obs = new double[L];

                SSwi = new double[L];
                SSwp = new double[L];
                SStot= new double[L];
                msl  = new double[L];
                mtl  = new double[L];

                //Initial freq
                if (!randfreq)
                {
                    SetVal(f, 1.0 / K);
                }
                else
                {
                    //Dirichlet(1,1,...,1)
                    SetVal(tf, 1.0 / K);
                    GetRandDirichlet(rnd, f, tf, 0.5);
                }

                //Subpop freq
                Fct = 1 - Math.Sqrt(1 - Fst);
                for (int c = 0, pid = 0; c < C; ++c)
                {
                    if (isregion)
                    {
                        GetRandDirichlet(rnd, tf, f, Fct);
                        for (int s = 0; s < Sc; ++s)
                            subpop[pid++] = new Population(rnd, pid, Ns, tf, Fct);
                    }
                    else for (int s = 0; s < Sc; ++s)
                        subpop[pid++] = new Population(rnd, pid, Ns, f, Fst);
                }

                if (preparecond) PrepareCondGenotype();
            }

            public Population(Random rnd, int _id, int _ninds, double[,] _f, double Fst)
            {
                //subpop
                id     = _id;
                ninds  = _ninds;
                inds      = new Individual[ninds];

                f     = new double[L, K];
                tf    = new double[L, K];
                f_obs = new double[L, K];

                px2_obs= new double[L];
                ps2_obs= new double[L];

                hs_obs= new double[L];
                hi_obs= new double[L];

                SSwi = new double[L];
                SSwp = new double[L];
                msl  = new double[L];

                GetRandDirichlet(rnd, f, _f, Fst);
                PrepareCondGenotype();
            }

            public Genotype DrawGenotype(Random rnd, int l)
            {
                double m = rnd.NextDouble();
                double[] gfz = this.gfz[l];
                int low = 0, high = gfz.Length - 1;

                while (low < high - 1)
                {
                    int mid = (low + high) >> 1;
                    if (gfz[mid] >= m)
                        high = mid;
                    else if (gfz[mid] < m)
                        low = mid;
                }
                while (high >= 1 && gfz[high - 1] >= m) high--;
                return genotypeid[high];
            }

            public Genotype DrawCondGenotype(Random rnd, Genotype f, int l)
            {
                double m = rnd.NextDouble();
                double[,] gfz = cond[l];
                int low = 0, high = gfz.GetLength(1) - 1, fid = f.id;

                while (low < high - 1)
                {
                    int mid = (low + high) >> 1;
                    if (gfz[fid, mid] >= m)
                        high = mid;
                    else if (gfz[fid, mid] < m)
                        low = mid;
                }
                while (high >= 1 && gfz[fid, high - 1] >= m) high--;
                return genotypeid[high];
            }

            public void PrepareCondGenotype()
            {
                int ngssr = (int)Binomial(v + K - 1, v);
                gfz = new double[L][];
                cond = new double[L][,];
                for (int l = 0; l < L; ++l)
                {
                    gfz[l] = new double[ngssr];
                    cond[l] = new double[ngssr, ngssr];

                    double cu1 = 0;
                    for (int i = 0; i < ngssr; ++i)
                    {
                        double cu2 = 0;
                        cu1 += genotypeid[i].GFZ(f, l, Fis);
                        gfz[l][i] = cu1;
                        for (int j = 0; j < ngssr; ++j)
                        {
                            cu2 += genotypeid[i].CondMother(genotypeid[j], f, l, Fis);
                            cond[l][i, j] = cu2;
                        }
                    }
                }
            }
            
            public void GenerateSubpopInds(Random rnd)
            {
                int oid = 0, pid = 0;
                int[] buf = new int[v];
                for (int i = 0; i < NF; ++i)
                {
                    Individual f = new Individual(rnd, pid++, this);
                    for (int j = 0; j < NM; ++j)
                    {
                        Individual m = new Individual(rnd, pid++, f, this);
                        for (int k = 0; k < NO; ++k)
                            inds[oid++] = new Individual(rnd, oid, this, f, m, buf);
                    }
                }

                //estimating obs allele frequency
                SetVal(f_obs, 0);
                for (int l = 0; l < L; ++l)
                {
                    for (int i = 0; i < ninds; ++i)
                    {
                        foreach (int a in inds[i].ssr[l].allele)
                            f_obs[l, a]++;
                        SSwi[l] += inds[i].ssr[l].SSwi();
                        hi_obs[l] += inds[i].ssr[l].GetHo();
                        px2_obs[l] += inds[i].ssr[l].GetPx2();
                    }
                }

                Mul(f_obs, f_obs, 1.0 / (ninds * v));
                Mul(hi_obs, hi_obs, 1.0 / ninds);
                Mul(px2_obs, px2_obs, 1.0 / ninds);

                for (int l = 0; l < L; ++l)
                {
                    double sswpl = Ns * v * Ns * v;
                    ps2_obs[l] = 0;
                    for (int k = 0; k < K; ++k)
                    {
                        ps2_obs[l] += f_obs[l, k] * f_obs[l, k];
                        double na = f_obs[l, k] * Ns * v;
                        sswpl -= na * na;
                    }
                    SSwp[l] = sswpl / (2 * Ns * v);
                    hs_obs[l] = 1.0 - ps2_obs[l];
                }
            }

            public void Inbreeding(Random rnd, 
                            ref double fis_Wright, ref double fis_Nei, ref double fis_WC, ref double fis_Corr,
                            ref double fst_Wright, ref double fst_Nei, ref double fst_WC, ref double fst_Corr,
                            ref double tfis_Wright, ref double tfis_Nei, ref double tfis_WC, 
                            ref double tfst_Wright, ref double tfst_Nei, ref double tfst_WC,
                            ref double fst_Hedrick, ref double tfst_Hedrick,
                            ref double fis_Huang, ref double fst_Huang)
            {
                //create individuals
                for (int ss = 0; ss < S; ++ss)
                {
                    Population s = subpop[ss];
                    s.GenerateSubpopInds(rnd);

                    //Copy inds
                    Array.Copy(s.inds, 0, inds, ss * Ns, Ns);

                    //Copy freq
                    Add(f_obs, f_obs, s.f_obs);
                    Add(hi_obs, hi_obs, s.hi_obs);
                    Add(hs_obs, hs_obs, s.hs_obs);
                    Add(px2_obs, px2_obs, s.px2_obs);
                    Add(ps2_obs, ps2_obs, s.ps2_obs);

                    Add(SSwp, SSwp, s.SSwp);
                    Add(SSwi, SSwi, s.SSwi);
                }

                //Unify freq
                Mul(f_obs, f_obs, 1.0 / S);
                Mul(hi_obs, hi_obs, 1.0 / S);
                Mul(hs_obs, hs_obs, 1.0 / S);
                Mul(px2_obs, px2_obs, 1.0 / S);
                Mul(ps2_obs, ps2_obs, 1.0 / S);

                for (int l = 0; l < L; ++l)
                {
                    double sstotl = Nt * v * Nt * v;
                    pt2_obs[l] = 0;
                    for (int k = 0; k < K; ++k)
                    {
                        pt2_obs[l] += f_obs[l, k] * f_obs[l, k];
                        double na = Nt * v * f_obs[l, k];
                        sstotl -= na * na;
                    }
                    SStot[l] = sstotl / (2 * Nt * v);
                    ht_obs[l] = 1.0 - pt2_obs[l];
                }

                //Corr
                double Fit = 1 - (1 - Fis) * (1 - Fst);
                double tx = 1.0 / v + (v - 1.0) / v * Fis;
                double tFS = (tx + Fis) / 2, tHS = (tx + 2 * Fis + Fis * Fis / tx) / 4;

                double pfs = Binomial(NO, 2) * NF * NM;
                double phs = Binomial(NO * NM, 2) * NF - pfs;
                pfs /= Binomial(Ns, 2);
                phs /= Binomial(Ns, 2);
                double txx = tFS * pfs + tHS * phs;
                double ts = 1.0 / Ns * tx + (Ns - 1.0) / Ns * txx;

                double vtxx = txx + (1 - txx) * Fst;
                double vtx = 1.0 / v + (v - 1.0) / v * Fit;
                double vts = ts + (1 - ts) * Fst;
                double vtss = C * Binomial(Sc, 2) / Binomial(S, 2) * Fct;
                double vtt = 1.0 / S * vts + (S - 1.0) / S * vtss;

                double txxinf = txx;
                double vtssinf = Fct / C;

                //Estimating Nes
                double tso = 0, txxo = 0;
                foreach (Population s in subpop)
                {
                    double ts_s = 0, tx_s = 0, txx_s = 0;
                    s.SubpopNe(ref ts_s, ref tx_s, ref txx_s);
                    tso += ts_s;
                    txxo += txx_s;
                    Add(msl, msl, s.msl);
                }
                Mul(msl, msl, 1.0 / S);
                tso /= S; txxo /= S; 

                //Estimating Net
                double vtto = 0, vtso = 0, vtsso = 0;
                TotpopNe(ref vtto, ref vtso, ref vtsso);

                //Theoritical bias
                tfis_Wright = (tx - ts) / (1 - ts);
                tfis_Wright = (vtx - vts) / (1 - vts);
                tfst_Wright = (vts - vtt) / (1 - vtt);

                tfis_Nei = (Fis - ts) / (1 - ts);
                tfis_Nei = (Fit - vts) / (1 - vts);
                tfst_Nei = (vts - vtt) / (1 - vtt);

                tfis_WC = (Fis - txx) / (1 - txx);
                tfis_WC = (Fit - vtxx) / (1 - vtxx);
                tfst_WC = (vtxx - vtss) / (1 - vtss);

                //estimate f
                int cl = 0;
                double Nei_hi = 0, Nei_hs = 0, Nei_ht = 0;
                double Wright1924f1 = 0, Wright1924f2 = 0, Wright1924f3 = 0, Wright1924f4 = 0;
                double Weir1984vwi = 0, Weir1984vai = 0, Weir1984vap = 0;
                double Huang2019f1 = 0, Huang2019f2 = 0, Huang2019f3 = 0, Huang2019f4 = 0;
                double Hedrick2005f1 = 0, Hedrick2005f2 = 0, Hedrick2005f3 = 0, Hedrick2005f4 = 0;

                double wrfis = 0, wrfst = 0;
                double neifis = 0, neifst = 0;
                double wcfis = 0, wcfst = 0;
                double huangfis = 0, huangfst = 0;
                double hedrickfst = 0;
                for (int l = 0; l < L; ++l)
                {
                    if (ht_obs[l] == 0) continue;

                    int obsk = 0;
                    for (int k = 0; k < K; ++k)
                        if (f_obs[l, k] > 0) obsk++;

                    //Nei
                    Nei_hi += hi_obs[l];
                    Nei_hs += hs_obs[l];
                    Nei_ht += ht_obs[l];

                    //Hedrick
                    Hedrick2005f1 += (ht_obs[l] - hs_obs[l]) * (S - 1 + hs_obs[l]);
                    Hedrick2005f2 += ht_obs[l] * (S - 1) * (1 - hs_obs[l]);
                    Hedrick2005f3 += (S - 1 + hs_obs[l]);
                    Hedrick2005f4 += (S - 1) * (1 - hs_obs[l]);

                    //Wright
                    Wright1924f1 += px2_obs[l] - ps2_obs[l];
                    Wright1924f2 += 1.0 - ps2_obs[l];
                    Wright1924f3 += ps2_obs[l] - pt2_obs[l];
                    Wright1924f4 += 1.0 - pt2_obs[l];

                    //Huang
                    Huang2019f1 += hi_obs[l];//hi 
                    Huang2019f2 += 1.0 - msl[l]; // hs
                    Huang2019f3 += 1.0 - msl[l];
                    Huang2019f4 += 1.0 - mtl[l];

                    //Weir
                    double vwi = SSwi[l] / (Nt * (v - 1));
                    double vai = (SSwp[l] - vwi * (Nt * v - S)) / (Nt * v - S * v);
                    double vap = (SStot[l] - vwi * (Nt * v - 1) - vai * (Nt * v - v)) / (Nt * v - Ns * v);
                    Weir1984vwi += vwi;
                    Weir1984vai += vai;
                    Weir1984vap += vap;
                    cl++;

                    if (true)
                    {
                        hedrickfst += (1 - hs_obs[l] / ht_obs[l]) * (S - 1 + hs_obs[l]) / ((S - 1) * (1 - hs_obs[l]));

                        neifis += 1 - hi_obs[l] / hs_obs[l];
                        neifst += 1 - hs_obs[l] / ht_obs[l];

                        wrfis += (px2_obs[l] - ps2_obs[l]) / (1.0 - ps2_obs[l]);
                        wrfst += (ps2_obs[l] - pt2_obs[l]) / (1.0 - pt2_obs[l]);

                        wcfis += vai / (vai + vwi);
                        wcfst += vap / (vai + vwi + vap);

                        huangfis += (px2_obs[l] - msl[l]) / (1.0 - msl[l]);
                        huangfst += (msl[l] - mtl[l]) / (1.0 - mtl[l]);
                    }
                }
                wrfis /= cl; wrfst /= cl; ;
                neifis /= cl; neifst /= cl;
                wcfis /= cl; wcfst /= cl;
                huangfis /= cl; huangfst /= cl;
                hedrickfst /= cl;
                tfst_Hedrick = hedrickfst;

                //Nei
                fis_Nei = 1.0 - Nei_hi / Nei_hs;
                fst_Nei = 1.0 - Nei_hs / Nei_ht;

                //Hedrick
                fst_Hedrick = Hedrick2005f1 / Hedrick2005f2;

                //Wright
                fis_Wright = Wright1924f1 / Wright1924f2;
                fst_Wright = Wright1924f3 / Wright1924f4;

                //Weir
                fis_WC = Weir1984vai / (Weir1984vwi + Weir1984vai);
                fst_WC = Weir1984vap / (Weir1984vwi + Weir1984vai + Weir1984vap);

                //Huang
                fis_Huang = 1 - Huang2019f1 / Huang2019f2;
                fst_Huang = 1 - Huang2019f3 / Huang2019f4;

                //Corr
                double fis_Nei_corr = fis_Nei * (1 - tso) + tso;
                double fst_Nei_corr = (fst_Nei * (1 - vtto) + vtto - tso) / (1 - tso);

                double fis_Wright_corr = (v * (fis_Wright * (1 - tso) + tso) - 1) / (v - 1);
                double fst_Wright_corr = (fst_Wright - tso + vtto - fst_Wright * vtto) / (1 - tso);

                double fis_WC_corr = fis_WC * (1 - txxo) + txxo;
                double fst_WC_corr = (fst_WC * (1 - vtsso) + vtsso - txxo) / (1 - txxo);

                //Corr2, not used
                fis_Corr = fis_Wright_corr;
                fst_Corr = fst_Wright_corr;

                //fis_Corr = fis_Huang;
                //fst_Corr = fst_Huang;

                if (S == 1)
                {
                    fst_Wright = 0;
                    fst_Nei = 0;
                    fst_WC = 0;
                    fst_Corr = 0;
                    fst_Huang = 0;

                    tfst_Nei = 0;
                    tfst_Wright = 0;
                    tfst_WC = 0;
                }
            }

            public struct d2
            {
                public double a;
                public double b;
            }

            private void SubpopNe(ref double ts, ref double tx, ref double txx)
            {
                double[,] indxy = new double[Ns, Ns];
                for (int x = 0; x < Ns; ++x)
                {
                    Individual xi = inds[x];
                    for (int y = x + 1; y < Ns; ++y)
                    {
                        Individual yi = inds[y];
                        double sum = 0;
                        for (int l = 0; l < L; ++l)
                        {
                            if (hs_obs[l] == 0) continue;
                            sum += xi.ssr[l].Pxy(yi.ssr[l]);
                        }
                        indxy[x, y] = indxy[y, x] = sum;
                    }
                }

                int nnr = (int)(Ns * (Ns - 1) / 2 * nrfreq);
                double Msl = 0, Hs = 0, Hi = 0;
                int cl = 0;

                for (int l = 0; l < L; ++l)
                {
                    if (hs_obs[l] == 0) { msl[l] = 1; continue; }
                    cl++;

                    Hs += hs_obs[l];
                    Hi += hi_obs[l];
                    List<d2> pxyl = new List<d2>();
                    for (int x = 0; x < Ns; ++x)
                    {
                        Individual xi = inds[x];
                        for (int y = x + 1; y < Ns; ++y)
                        {
                            Individual yi = inds[y];
                            d2 c;
                            c.b = xi.ssr[l].Pxy(yi.ssr[l]);
                            c.a = indxy[x, y] - c.b;
                            pxyl.Add(c);
                        }
                    }

                    pxyl = pxyl.OrderBy(o => o.a).ToList();
                    pxyl.RemoveRange(nnr, pxyl.Count - nnr);
                    msl[l] = pxyl.Sum(o => o.b) / nnr;
                    Msl += msl[l];
                }

                ts = 1 - Hs / (cl - Msl);
                tx = 1.0 / v + (v - 1.0) / v * (1 - Hi / (cl - Msl));
                txx = (Ns * ts - tx) / (Ns - 1);
            }

            public double Pxy(Population y, int l)
            {
                double s = 0;
                for (int k = 0; k < K; ++k)
                    s += f_obs[l, k] * y.f_obs[l, k];
                return s;
            }

            private void TotpopNe(ref double tt, ref double ts, ref double tss)
            {
                if (S == 1) return;
                double[,] sxy = new double[S, S];
                double[,,] sxyl = new double[S, S, L];

                for (int x = 0; x < S; ++x)
                {
                    Population xi = subpop[x];
                    for (int y = x + 1; y < S; ++y)
                    {
                        Population yi = subpop[y];
                        double sum = 0;
                        for (int l = 0; l < L; ++l)
                        {
                            if (ht_obs[l] == 0) continue;
                            double sl = xi.Pxy(yi, l);
                            sxyl[x, y, l] = sxyl[y, x, l] = sl;
                            sum += sl;
                        }
                        sxy[x, y] = sxy[y, x] = sum;
                    }
                }

                int nnr = (int)(S * (S - 1) / 2 * nrfreq);
                double Mtl = 0, Hs = 0, Ht = 0;
                int cl = 0;

                for (int l = 0; l < L; ++l)
                {
                    if (ht_obs[l] == 0) { mtl[l] = 1; continue; }
                    cl++;

                    Hs += hs_obs[l];
                    Ht += ht_obs[l];
                    List<d2> pxyl = new List<d2>();
                    for (int x = 0; x < S; ++x)
                    {
                        Population xi = subpop[x];
                        for (int y = x + 1; y < S; ++y)
                        {
                            Population yi = subpop[y];
                            d2 c;
                            c.b = sxyl[x, y, l];
                            c.a = sxy[x, y] - sxyl[x, y, l];
                            pxyl.Add(c);
                        }
                    }

                    pxyl = pxyl.OrderBy(o => o.a).ToList();
                    pxyl.RemoveRange(nnr, pxyl.Count - nnr);
                    mtl[l] = pxyl.Sum(o => o.b) / nnr;
                    Mtl += mtl[l];
                }

                tt = 1 - Ht / (cl - Mtl);
                ts = 1 - Hs / (cl - Mtl);
                tss = (Nt * Nt * tt - S * Ns * Ns * ts) / (Nt * Nt - S * Ns * Ns);
            }
        };
        
        #endregion

        string[] ARGS;
        Process[] PROC;
        private void RunCommands(object sender, EventArgs e)
        {
            //this.ShowIcon = false;
            this.WindowState = FormWindowState.Minimized;
            int ncore = (int)(Environment.ProcessorCount * 1.2);
            if (Environment.ProcessorCount < 13) ncore = Environment.ProcessorCount / 2;
            Random rnd = new Random();
            PROC = new Process[ncore];
            string[] name = new string[ncore];
            DateTime[] time = new DateTime[ncore];

            for (int i = 0; i < ncore; ++i)
            {
                time[i] = DateTime.Now;
                name[i] = "";
            }

            string APP = "";
            for (; ; )
            {
                string[] cmd2 = textBox1.Text.Split(new char[] { '\r', '\n' }, StringSplitOptions.RemoveEmptyEntries);
                List<string> cmd = new List<string>();

                for (int i = 0; i < cmd2.Length; ++i)
                {
                    string[] con = cmd2[i % cmd2.Length].Split(new char[] { ' ', '\t' }, StringSplitOptions.RemoveEmptyEntries);
                    APP = con[0] + "\\";
                    if (!File.Exists(APP + con[1] + ".txt"))
                    {
                        cmd.Add(cmd2[i % cmd2.Length].Replace("\t", " "));
                        continue;
                    }

                    int ncalced = File.ReadAllLines(APP + con[1] + ".txt").Length + 1;

                    if (con[0].ToLower() == "app1" && int.Parse(con[13]) > ncalced)
                        cmd.Add(cmd2[i % cmd2.Length].Replace("\t", " "));

                    if (con[0].ToLower() == "app2" && int.Parse(con[9]) > ncalced)
                        cmd.Add(cmd2[i % cmd2.Length].Replace("\t", " "));
                }

                if (cmd.Count == 0) break;

                int ll = cmd.Count;
                while (cmd.Count < ncore)
                    cmd.Add(cmd[rnd.Next(ll)]);

                int[] seq = new int[cmd.Count];
                for (int i = 0; i < cmd.Count; ++i)
                    seq[i] = i;
                Permute(seq, rnd);

                int cmdp = 0;
                string path = Environment.GetCommandLineArgs()[0];

                while (cmdp < cmd.Count)
                {
                    for (int i = 0; i < ncore; ++i)
                    {
                        this.Text = cmdp + " / " + cmd.Count;
                        if (cmdp >= seq.Length) continue;

                        if (PROC[i] == null || PROC[i].HasExited)
                        {
                            PROC[i] = Process.Start("fst.exe", cmd[seq[cmdp++]]);
                            //7PROC[i] = Process.Start(cmd[0][3] == '1' ? path.Replace(".vshost", "") : "PloidyInferAPP2.exe", cmd[seq[cmdp++]]);
                            time[i] = DateTime.Now;
                            Thread.Sleep(300);
                        }

                        if ((DateTime.Now - time[i]).Minutes >= 6)
                        {
                            time[i] = DateTime.Now;

                            StringBuilder sb = new StringBuilder(5120);
                            GetWindowText(PROC[i].MainWindowHandle, sb, sb.Capacity);

                            if (name[i] != sb.ToString())
                                name[i] = sb.ToString();
                            else
                            {
                                while (!PROC[i].HasExited)
                                {
                                    Process.Start("ntsd.exe", "-c q -p " + PROC[i].Id);
                                    Thread.Sleep(1000);
                                    Process.Start("taskkill.exe", "/IM ntsd.exe /f ");
                                    Thread.Sleep(1000);
                                }
                            }
                        }

                    }
                    Thread.Sleep(1000);
                }
            }
        }

        public Form1(string[] args)
        {
            ARGS = args;
            InitializeComponent();
        }

        private void toolStripButton2_Click(object sender, EventArgs e)
        {
            //app1
            if (Directory.Exists("app1"))
            {
                int nst = 13;
                DirectoryInfo di = new DirectoryInfo("app1");
                StringBuilder sb = new StringBuilder();
                Parallel.ForEach(di.GetFiles("*.txt"), new ParallelOptions() { MaxDegreeOfParallelism = 12 }, f =>
                //foreach (FileInfo f in di.GetFiles("*.txt"))
                {
                    if (f.Name[0] >= '0' && f.Name[0] <= '9')
                    {
                        string[] cont = File.ReadAllLines(f.FullName);

                        int npar = 18;
                        double[] ex = new double[npar];
                        double[] ex2 = new double[npar]; 
                        int[] nl = new int[npar];
                        foreach (string line in cont)
                        {
                            string[] grids = line.Split(new char[] { '\t' }, StringSplitOptions.None);
                            for (int i = 0; i < npar; ++i)
                            {
                                double v = double.Parse(grids[nst + i]);
                                if (double.IsNaN(v) || double.IsInfinity(v)) continue;
                                nl[i]++;
                                ex[i] += v;
                                ex2[i] += v * v;
                            }
                        }

                        for (int i = 0; i < npar; ++i)
                        {
                            ex[i] /= nl[i];
                            ex2[i] /= nl[i];
                        }

                        string[] grids2 = cont[0].Split(new char[] { '\t' }, StringSplitOptions.None);

                        lock (this)
                        {
                            sb.Append(grids2[0]);
                            for (int i = 1; i < nst; ++i)
                                sb.Append("\t" + grids2[i]);
                            for (int i = 0; i < npar; ++i)
                            {
                                if (double.IsInfinity(ex[i]))
                                    sb.Append("\tinf\tinf");
                                else if (double.IsNaN(ex[i]))
                                    sb.Append("\tnan\tnan");
                                else
                                {
                                    sb.Append("\t" + ex[i]);
                                    sb.Append("\t" + (ex2[i] - ex[i] * ex[i]));
                                }
                            }
                            sb.Append("\r\n");
                        }
                    }
                });
                File.WriteAllText("app1.txt", sb.ToString());
            }
        }

        private void APP1()
        {
            this.BeginInvoke(new Action(() => { this.WindowState = FormWindowState.Minimized; }));

            //app	id	v	L	k	C	Sc	Fis	Fst	Nf	Nm	No	nrep
            int c = 1;
            int id = int.Parse(ARGS[c++]);
            v = int.Parse(ARGS[c++]);
            L = int.Parse(ARGS[c++]);
            K = int.Parse(ARGS[c++]);
            C = int.Parse(ARGS[c++]);
            Sc = int.Parse(ARGS[c++]);
            Fis = double.Parse(ARGS[c++]);
            Fst = double.Parse(ARGS[c++]);
            NF = int.Parse(ARGS[c++]);
            NM = int.Parse(ARGS[c++]);
            NO = int.Parse(ARGS[c++]);
            nrfreq = double.Parse(ARGS[c++]);
            Nrep = int.Parse(ARGS[c++]);
            S = C * Sc;
            Ns = NF * NM * NO;
            Nt = S * Ns;

            int nc = 0;

            double Fis_Wright = 0, Fis_Nei = 0, Fis_WC = 0, Fis_Corr = 0;
            double Fis2_Wright = 0, Fis2_Nei = 0, Fis2_WC = 0, Fis2_Corr = 0;
            double Fst_Wright = 0, Fst_Nei = 0, Fst_WC = 0, Fst_Corr = 0;
            double Fst2_Wright = 0, Fst2_Nei = 0, Fst2_WC = 0, Fst2_Corr = 0;

            double tFis_Wright = 0, tFis_Nei = 0, tFis_WC = 0;
            double tFst_Wright = 0, tFst_Nei = 0, tFst_WC = 0;
            double Fst_Hedrick = 0, tFst_Hedrick = 0;
            double Fst_Huang = 0, Fis_Huang = 0;

            PrepareGenotypes();

            string app = ARGS[0];

            if (!Directory.Exists(app)) Directory.CreateDirectory(app);
            string path = app + "\\" + id + ".txt";
            bool debug = false;
            if (debug && File.Exists(path))
                File.Delete(path);

            StringBuilder sb = new StringBuilder();
            int seed = DateTime.Now.DayOfYear;
            seed *= 1000;
            seed += DateTime.Now.Hour;
            seed *= 100;
            seed += DateTime.Now.Minute;
            seed *= 100;
            seed += DateTime.Now.Second;
            seed *= 1000;
            seed += DateTime.Now.Millisecond;

            Random rnd2 = new Random(seed);
            this.BeginInvoke(new Action(() => { this.Text = id + " 0/0"; }));

            Parallel.For(0, Nrep, new ParallelOptions() { MaxDegreeOfParallelism = debug ? 1 : 1 }, m => 
            //for (int m = 0; m < nrep; ++m)
            {
                Random rnd;

                double fis_Wright = 0, fis_Nei = 0, fis_WC = 0, fis_Corr = 0, fis_Huang = 0;
                double fst_Wright = 0, fst_Nei = 0, fst_WC = 0, fst_Corr = 0, fst_Huang = 0;
                double tfis_Wright = 0, tfis_Nei = 0, tfis_WC = 0;
                double tfst_Wright = 0, tfst_Nei = 0, tfst_WC = 0;
                double fst_Hedrick = 0, tfst_Hedrick = 0;

                if (debug) rnd = new Random(rnd2.Next());
                else       rnd = rnd2;
                Population tot = new Population(rnd, true, true, true);

                tot.Inbreeding(rnd, ref fis_Wright, ref fis_Nei, ref fis_WC, ref fis_Corr,
                                    ref fst_Wright, ref fst_Nei, ref fst_WC, ref fst_Corr, 
                                    ref tfis_Wright, ref tfis_Nei, ref tfis_WC, 
                                    ref tfst_Wright, ref tfst_Nei, ref tfst_WC, 
                                    ref fst_Hedrick, ref tfst_Hedrick,
                                    ref fis_Huang, ref fst_Huang);

                lock (this)
                {
                    Fis_Wright += fis_Wright; Fis2_Wright += fis_Wright * fis_Wright;
                    Fis_Nei += fis_Nei; Fis2_Nei += fis_Nei * fis_Nei;
                    Fis_WC += fis_WC; Fis2_WC += fis_WC * fis_WC;
                    Fis_Corr += fis_Corr; Fis2_Corr += fis_Corr * fis_Corr;

                    Fst_Wright += fst_Wright; Fst2_Wright += fst_Wright * fst_Wright;
                    Fst_Nei += fst_Nei; Fst2_Nei += fst_Nei * fst_Nei;
                    Fst_WC += fst_WC; Fst2_WC += fst_WC * fst_WC;

                    Fst_Corr += fst_Corr; Fst2_Corr += fst_Corr * fst_Corr;

                    tFis_Wright = tfis_Wright;
                    tFis_Nei = tfis_Nei;
                    tFis_WC = tfis_WC;

                    tFst_Wright = tfst_Wright;
                    tFst_Nei = tfst_Nei;
                    tFst_WC = tfst_WC;

                    Fst_Hedrick += fst_Hedrick;
                    tFst_Hedrick = tfst_Hedrick;

                    Fst_Huang += fst_Huang;
                    Fis_Huang += fis_Huang;
                    
                    nc++;
                }

                if (!debug)
                {
                    //app	id	v	h2	S	qst_	Qis	Lqtl	Kqtl	L	K	NF	NM	NO	nrep
                    sb.Append(
                        id + "\t" +
                        v + "\t" +
                        L + "\t" +
                        K + "\t" +
                        C + "\t" +
                        Sc + "\t" +
                        Fis + "\t" +
                        Fst + "\t" +
                        NF + "\t" +
                        NM + "\t" +
                        NO + "\t" +
                        nrfreq + "\t" +
                        Nrep + "\t" +
                        fis_Wright.ToString("F10") + "\t" +
                        fis_Nei.ToString("F10") + "\t" +
                        fis_WC.ToString("F10") + "\t" +
                        fis_Corr.ToString("F10") + "\t" +
                        fst_Wright.ToString("F10") + "\t" +
                        fst_Nei.ToString("F10") + "\t" +
                        fst_WC.ToString("F10") + "\t" +
                        fst_Corr.ToString("F10") + "\t" +
                        tfis_Wright.ToString("F10") + "\t" +
                        tfis_Nei.ToString("F10") + "\t" +
                        tfis_WC.ToString("F10") + "\t" +
                        tfst_Wright.ToString("F10") + "\t" +
                        tfst_Nei.ToString("F10") + "\t" +
                        tfst_WC.ToString("F10") + "\t" +
                        fis_Huang.ToString("F10") + "\t" +
                        fst_Huang.ToString("F10") + "\t" +
                        fst_Hedrick.ToString("F10") + "\t" +
                        tfst_Hedrick.ToString("F10") + "\r\n");

                    if (rnd.NextDouble() < 0.3)
                    {
                        for (; ; )
                        {
                            try
                            {
                                File.AppendAllText(path, sb.ToString());
                                sb.Clear();
                                break;
                            }
                            catch { }
                        }
                        int line = File.ReadAllLines(path).Length;
                        this.BeginInvoke(new Action(() => { this.Text = id + " " + line + "/" + Nrep + " "; }));
                        if (line >= Nrep) Environment.Exit(0);
                    }
                }
            });

            if (!debug) for (; ; )
            {
                try
                {
                    File.AppendAllText(path, sb.ToString());
                    sb.Clear();
                    break;
                }
                catch { }
            }

            Fis_Wright /= nc; Fis2_Wright /= nc; Fis2_Wright -= Fis_Wright * Fis_Wright;
            Fis_Nei /= nc; Fis2_Nei /= nc; Fis2_Nei -= Fis_Nei * Fis_Nei;
            Fis_WC /= nc; Fis2_WC /= nc; Fis2_WC -= Fis_WC * Fis_WC;
            Fis_Corr /= nc; Fis2_Corr /= nc; Fis2_Corr -= Fis_Corr * Fis_Corr;

            Fst_Wright /= nc; Fst2_Wright /= nc; Fst2_Wright -= Fst_Wright * Fst_Wright;
            Fst_Nei /= nc; Fst2_Nei /= nc; Fst2_Nei -= Fst_Nei * Fst_Nei;
            Fst_WC /= nc; Fst2_WC /= nc; Fst2_WC -= Fst_WC * Fst_WC;
            Fst_Corr /= nc; Fst2_Corr /= nc; Fst2_Corr -= Fst_Corr * Fst_Corr;

            Fis_Huang /= nc; Fst_Huang /= nc;

            string s =
                Fis_Wright + "\t" + Fis2_Wright + "\r\n" +
                Fis_Nei + "\t" + Fis2_Nei + "\r\n" +
                Fis_WC + "\t" + Fis2_WC + "\r\n" +
                Fis_Corr + "\t" + Fis2_Corr + "\r\n" +
                Fst_Wright + "\t" + Fst2_Wright + "\r\n" +
                Fst_Nei + "\t" + Fst2_Nei + "\r\n" +
                Fst_WC + "\t" + Fst2_WC + "\r\n" +
                Fst_Corr + "\t" + Fst2_Corr + "\r\n";

            Environment.Exit(0);
        }

        private void Form1_Load(object sender, EventArgs e)
        {
            if (ARGS != null && ARGS.Length > 0 && ARGS[0].ToLower().Substring(0,3) == "app")
            {
                Thread thread = new Thread(new ThreadStart(APP1));
                thread.Start();
            }
        }
    }
}
