using System;

namespace Science.Mathematics.VectorCalculus
{
	/// <summary>
	/// Adaptive multidimensional Monte Carlo integration.
	/// </summary>
	public class IntegrationMultiD
	{
        private Function.ToLastType<double[],double> del;
        public IntegrationMultiD()
        {
        }
        public Function.ToLastType<double[], double> Function
        {
            set{del = value;}
        }

        public IntegrationMultiD(Function.ToLastType<double[], double> f)
		{
			del = f;
		}

        public IntegrationMultiD(Function.ToLastType<double[], double> f, double[] from, double[] to)
        {
            del = f;
            this.from = from;
            this.to = to;
        }
	
		public void Compute()
		{
            input = new double[2 * to.Length];
            for (int i = 0; i < to.Length; i++)
            {
                input[i] = from[i];
                input[to.Length + i] = to[i];
            }
            Vegas obj = new Vegas();      
            obj.Dummy = idum;
            obj.vegas(input, input.Length / 2, del, 0, nc, 5, 0);

            tgral = obj.BestEstimation;
            sd = obj.StandardDeviation;
            chi2a = obj.ChiSquare;
            result = obj.Result;            	                
		}

		private double[] input, from, to;
   		public double[] From
		{
            set { from = value;}
		}
        public double[] To
        {
            set { to = value;}
        }

		private ulong nc=100000;
		public ulong NumberOfCall
		{
			set{nc=value;}
		}

        private ulong idum = 17;
        public ulong Seed
        {
            set { idum = value; }
        }

        private double tgral, sd, chi2a;
        private string result;
        public double BestEstimation
        {
            get { return tgral; }
        }
        public double StandardDeviation
        {
            get { return sd; }
        }
        public double ChiSQUARE
        {
            get { return chi2a; }
        }
        public string Result
        {
            get { return result; }
        }

        public override string ToString()
        {
            string res = "";
            res += this.BestEstimation.ToString() + " +/- " +
                this.StandardDeviation.ToString();
            res += "\r\n";
            return res;
        } 
	}


    class Ran
    {
        ulong u, v, w;
        public Ran(ulong j)
        {
            v = 4101842887655102017L;
            w = 1;
            u = j ^ v; int64();
            v = u; int64();
            w = v; int64();
        }
        public ulong int64()
        {
            u = u * 2862933555777941757L + 7046029254386353087L;
            v ^= v >> 17; v ^= v << 31; v ^= v >> 8;
            w = 4294957665U * (w & 0xffffffff) + (w >> 32);
            ulong x = u ^ (u << 21); x ^= x >> 35; x ^= x << 4;
            return (x + v) ^ w;
        }
        public double doub()
        {
            return 5.42101086242752217E-20 * int64();
        }
        public uint int32()
        {
            return (uint)int64();
        }
        ulong breg;
        int bc;
        public byte int8()
        {
            if ((bc--) != 0) return (byte)(breg >>= 8);
            breg = int64();
            bc = 7;
            return (byte)breg;
        }
    }

    class Vegas
	{
		public Vegas()
		{
			ia = new int[MXDIM];
			kg = new int[MXDIM];
			d = new double[NDMX,MXDIM];
			di = new double[NDMX,MXDIM];
			dt = new double[MXDIM];
            dx = new double[MXDIM];
			r = new double[NDMX];
			x = new double[MXDIM];
			xi= new double[MXDIM,NDMX];
			xin= new double[NDMX];
		}
        private double ALPH = 1.5;
        private int NDMX = 50;
        private int MXDIM = 10;
        private double TINY = 1.0e-30;

		private int i,it,j,k,mds,nd,ndo,ng,npg;
		private int[] ia, kg;
		private double calls,dv2g,dxg,f,f2,f2b,fb,rc,ti,tsi,wgt,xjac,xn,xnd,xo;
		private double[,] d, di, xi;
		private double[] dt, dx, r, x, xin;
		private double schi,si,swgt;
	
		private ulong idum = 0;
		public ulong Dummy
		{
			set{idum = value;}
		}
		private double tgral, sd, chi2a;
		private string result;
		public double BestEstimation
		{
			get{return tgral;}
		}
		public double StandardDeviation
		{
			get{return sd;}
		}
		public double ChiSquare
		{
			get{return chi2a;}
		}
		public string Result
		{
			get{return result;}
		}

        public void vegas(double[] regn, int ndim,
            Function.ToLastType<double[],double> fxn, int init, ulong ncall,
            int itmx, int nprn)
        {
            Ran obj = new Ran(idum);
   
            if (init <= 0)
            {
                mds = ndo = 1;
                for (j = 0; j < ndim; j++) xi[j, 0] = 1.0;
            }
            if (init <= 1) si = swgt = schi = 0.0;
            if (init <= 2)
            {
                nd = NDMX;
                ng = 1;
                if (mds != 0)
                {
                    ng = (int)Math.Pow(ncall / 2.0 + 0.25, 1.0 / ndim);
                    mds = 1;
                    if ((2 * ng - NDMX) >= 0)
                    {
                        mds = -1;
                        npg = ng / NDMX + 1;
                        nd = ng / npg;
                        ng = npg * nd;
                    }
                }
                for (k = 1, i = 0; i < ndim; i++) k *= ng;
                npg = Math.Max((int)ncall / k, 2);
                calls = npg * k;
                dxg = 1.0 / ng;
                for (dv2g = 1, i = 0; i < ndim; i++) dv2g *= dxg;
                dv2g = Math.Pow(calls * dv2g, 2.0) / npg / npg / (npg - 1.0);
                xnd = nd;
                dxg *= xnd;
                xjac = 1.0 / calls;
                for (j = 0; j < ndim; j++)
                {
                    dx[j] = regn[j + ndim] - regn[j];
                    xjac *= dx[j];
                }
                if (nd != ndo)
                {
                    for (i = 0; i < nd; i++) r[i] = 1.0;
                    for (j = 0; j < ndim; j++) rebin(ndo / xnd, nd, r, xin, xi, j);
                    ndo = nd;
                }
                if (nprn >= 0)
                {
                    result += "Input parameters for vegas:" + "ndim=" + Convert.ToString(ndim) + " ncall=" + Convert.ToString(calls) + "\r\n";
                    result += "   it=" + Convert.ToString(it) + " itmx=" + Convert.ToString(itmx) + "\r\n";
                    result += "   nprn=" + Convert.ToString(nprn) + " ALPH=" + Convert.ToString(ALPH) + "\r\n";
                    result += "   mds=" + Convert.ToString(mds) + " nd=" + Convert.ToString(nd) + "\r\n";
                    for (j = 0; j < ndim; j++)
                    {
                        result += "   xl[" + Convert.ToString(j) + "]=" + Convert.ToString(regn[j]) +
                                 "   xu[" + Convert.ToString(j) + "]=" + Convert.ToString(regn[j + ndim]) + "\r\n";
                    }
                }
            }
            for (it = 1; it <= itmx; it++)
            {
                ti = tsi = 0.0;
                for (j = 0; j < ndim; j++)
                {
                    kg[j] = 1;
                    for (i = 0; i < nd; i++) d[i, j] = di[i, j] = 0.0;
                }
                for (; ; )
                {
                    fb = f2b = 0.0;
                    for (k = 1; k <= npg; k++)
                    {
                        wgt = xjac;
                        for (j = 0; j < ndim; j++)
                        {
                            xn = (kg[j] - obj.doub()) * dxg + 1.0;
                            ia[j] = Math.Max(Math.Min((int)(xn), NDMX), 1);
                            if (ia[j] > 1)
                            {
                                xo = xi[j, ia[j] - 1] - xi[j, ia[j] - 2];
                                rc = xi[j, ia[j] - 2] + (xn - ia[j]) * xo;
                            }
                            else
                            {
                                xo = xi[j, ia[j] - 1];
                                rc = (xn - ia[j]) * xo;
                            }
                            x[j] = regn[j] + rc * dx[j];
                            wgt *= xo * xnd;
                        }
                        f = wgt * fxn(x); // f=wgt*fxn(x,wgt)
                        f2 = f * f;
                        fb += f;
                        f2b += f2;
                        for (j = 0; j < ndim; j++)
                        {
                            di[ia[j] - 1, j] += f;
                            if (mds >= 0) d[ia[j] - 1, j] += f2;
                        }
                    }
                    f2b = Math.Sqrt(f2b * npg);
                    f2b = (f2b - fb) * (f2b + fb);
                    if (f2b <= 0.0) f2b = TINY;
                    ti += fb;
                    tsi += f2b;
                    if (mds < 0)
                    {
                        for (j = 0; j < ndim; j++) d[ia[j] - 1, j] += f2b;
                    }
                    for (k = ndim; k >= 1; k--)
                    {
                        kg[k - 1] %= ng;
                        if (++kg[k - 1] != 1) break;
                    }
                    if (k < 1) break;
                }
                tsi *= dv2g;
                wgt = 1.0 / tsi;
                si += wgt * ti;
                schi += wgt * ti * ti;
                swgt += wgt;
                tgral = si / swgt;
                chi2a = (schi - si * tgral) / (it - 0.9999);
                if (chi2a < 0.0) chi2a = 0.0;
                sd = Math.Sqrt(1.0 / swgt);
                tsi = Math.Sqrt(tsi);
                if (nprn >= 0)
                {
                    result += " iteration no. " + Convert.ToString(it) + ": integral =" + Convert.ToString(ti) +
                        "   +/-  " + Convert.ToString(tsi) + "\r\n";
                    result += " all iterations: integral=" + Convert.ToString(tgral) + "+/-" + Convert.ToString(sd) +
                        "  chi**2/IT n = " + Convert.ToString(chi2a) + "\r\n";

                    if (nprn != 0)
                    {
                        for (j = 0; j < ndim; j++)
                        {
                            result += " DATA FOR axis " + Convert.ToString(j + 1) + "\r\n";
                            result += "     X     " + "    delta i   " + "     X     " + "    delta i   " +
                                "     X     " + "    delta i   " + "\r\n";
                            for (i = nprn / 2; i < nd; i += nprn + 2)
                            {
                                result += " " + Convert.ToString(xi[j, i]) + " " + Convert.ToString(di[i, j]) +
                                    " " + Convert.ToString(xi[j, i + 1]) + " " + Convert.ToString(di[i + 1, j]) +
                                    " " + Convert.ToString(xi[j, i + 2]) + " " + Convert.ToString(di[i + 2, j]) + "\r\n";
                            }
                        }
                    }
                }
                for (j = 0; j < ndim; j++)
                {
                    xo = d[0, j];
                    xn = d[1, j];
                    d[0, j] = (xo + xn) / 2.0;
                    dt[j] = d[0, j];
                    for (i = 1; i < nd - 1; i++)
                    {
                        rc = xo + xn;
                        xo = xn;
                        xn = d[i + 1, j];
                        d[i, j] = (rc + xn) / 3.0;
                        dt[j] += d[i, j];
                    }
                    d[nd - 1, j] = (xo + xn) / 2.0;
                    dt[j] += d[nd - 1, j];
                }
                for (j = 0; j < ndim; j++)
                {
                    rc = 0.0;
                    for (i = 0; i < nd; i++)
                    {
                        if (d[i, j] < TINY) d[i, j] = TINY;
                        r[i] = Math.Pow((1.0 - d[i, j] / dt[j]) /
                            (Math.Log(dt[j]) - Math.Log(d[i, j])), ALPH);
                        rc += r[i];
                    }
                    rebin(rc / xnd, nd, r, xin, xi, j);
                }
            }
        }       
	   
        private void rebin(double rc, int nd, double[] r, double[] xin, double[,] xi, int j)
        {
            int i, k = -1;
            double dr = 0.0, xn = 0.0, xo = 0.0;

            for (i = 0; i < nd - 1; i++)
            {
                while (rc > dr) dr += r[++k];
                if (k > 0) xo = xi[j, k - 1];
                xn = xi[j, k];
                dr -= rc;
                xin[i] = xn - (xn - xo) * dr / r[k];
            }
            for (i = 0; i < nd - 1; i++) xi[j, i] = xin[i];
            xi[j, nd - 1] = 1.0;
        }
    }
}