using MathNet.Numerics.LinearAlgebra;
using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using SMAT_CE.Matlab;
using MathNet.Numerics;
using System.Diagnostics;
using System.IO;

namespace SMAT_CE.DownscaleR.Ozone
{
    /// <summary>
    /// Ozone downscaler analysis
    /// </summary>
    public class OzoneAnalysis
    {
        private CommonSettings _settings;
        private string _innermsg;
        // note: please store the "id", "lat", "lon", "ozonebase" in order
        //private Tuple<List<string>, List<double>, List<double>, List<double>> _monitordatasfromsmat;
        private Tuple<string[], double[], double[], double[]> _monitordatasfromsmat;
        // note: please store the "id", "lat", "lon", "ozone" in order
        //private Tuple<List<string>, List<double>, List<double>, List<double>> _modeldatasfromsmat;
        private Tuple<string[], double[], double[], double[]> _modeldatasfromsmat;
        private List<MCMC> _mcmcoutputlist;
        private List<OzoneModel> _resultoutputlist;

        private MCMC _mcmcoutputitem;
        private OzoneMonitor _monitoritem;
        private OzoneModel _modelitem;
        private int count;
        private List<OzoneMonitor> _monitorlistfromsmat;
        private Matrix<double> _matrixmonitorlatlon;
        private Matrix<double> _matrixmonitorozone;
        private Matrix<double> _matrixmodellatlon;
        private Matrix<double> _matrixmodelozone;
        private bool ok;

        private OzoneAnalysis(CommonSettings settings)
        {
            _settings = settings;
            _mcmcoutputlist = new List<MCMC>();
            _resultoutputlist = new List<OzoneModel>();
            _monitorlistfromsmat = new List<OzoneMonitor>();
        }
        public OzoneAnalysis(
            CommonSettings settings,
            Tuple<string[], double[], double[], double[]> monitordatas,
            Tuple<string[], double[], double[], double[]> modeldatas
            )
            : this(settings)
        {
            _monitordatasfromsmat = monitordatas;
            _modeldatasfromsmat = modeldatas;
        }

        /// <summary>
        /// Downscaler settings
        /// </summary>
        public CommonSettings Settings
        {
            get { return _settings; }
        }
        /// <summary>
        /// Running msg
        /// </summary>
        public string CurrentMsg
        {
            get { return _innermsg; }
        }
        /// <summary>
        /// MCMC output datas
        /// </summary>
        public List<MCMC> MCMCOutputData
        {
            get { return _mcmcoutputlist; }
        }
        /// <summary>
        /// Result output datas
        /// </summary>
        public List<OzoneModel> ResultOutputData
        {
            get { return _resultoutputlist; }
        }

        /// <summary>
        /// Handle SMAT input data, include matrixing
        /// </summary>
        private bool HandleSMATInputData()
        {
            try
            {
                if (_monitordatasfromsmat == null || _modeldatasfromsmat == null)
                {
                    _innermsg = "The input data is null";
                    return false;
                }
                string[] lstid;
                double[] lstlat, lstlon, lstconc;
                string id;
                double lat, lon, conc;
                List<double> lstlatlon = new List<double>();
                List<double> lstozone = new List<double>();
                // handle the monitor data from smat
                lstid = _monitordatasfromsmat.Item1;
                lstlat = _monitordatasfromsmat.Item2;
                lstlon = _monitordatasfromsmat.Item3;
                lstconc = _monitordatasfromsmat.Item4;
                count = lstid.Length;
                for (int i = 0; i < count; i++)
                {
                    id = lstid[i];
                    lat = lstlat[i];
                    lon = lstlon[i];
                    conc = lstconc[i];
                    _monitoritem = new OzoneMonitor(id, lat, lon, conc);
                    _monitorlistfromsmat.Add(_monitoritem);
                    // filter data
                    if (conc < 0.0)
                        continue;
                    if (conc == 0.0 && conc + 1.0 == 1.0)
                        conc = 0.1;
                    lstlatlon.Add(lat);
                    lstlatlon.Add(lon);
                    lstozone.Add(conc);
                }
                if (lstconc.Length <= 0)
                    return false;
                _matrixmonitorlatlon = Matrix<double>.Build.DenseOfColumnMajor(2, lstlatlon.Count / 2, lstlatlon.ToArray());
                _matrixmonitorlatlon = _matrixmonitorlatlon.Transpose();
                _matrixmonitorozone = Matrix<double>.Build.DenseOfColumnMajor(lstozone.Count, 1, lstozone.ToArray());
                // handle the model data from smat
                lstid = _modeldatasfromsmat.Item1;
                lstlat = _modeldatasfromsmat.Item2;
                lstlon = _modeldatasfromsmat.Item3;
                lstconc = _modeldatasfromsmat.Item4;
                count = lstid.Length;
                lstlatlon.Clear();
                lstozone.Clear();
                for (int i = 0; i < count; i++)
                {
                    id = lstid[i];
                    lat = lstlat[i];
                    lon = lstlon[i];
                    conc = lstconc[i];
                    _modelitem = new OzoneModel(id, lat, lon, conc);
                    _resultoutputlist.Add(_modelitem);
                    // filter data
                    if (conc < 0.0)
                        continue;
                    if (conc == 0.0 && conc + 1.0 == 1.0)
                        conc = 0.1;
                    lstlatlon.Add(lat);
                    lstlatlon.Add(lon);
                    lstozone.Add(conc);
                }
                if (lstconc.Length <= 0)
                    return false;
                _matrixmodellatlon = Matrix<double>.Build.DenseOfColumnMajor(2, lstlatlon.Count / 2, lstlatlon.ToArray());
                _matrixmodellatlon = _matrixmodellatlon.Transpose();
                _matrixmodelozone = Matrix<double>.Build.DenseOfColumnMajor(lstozone.Count, 1, lstozone.ToArray());
                // clear list
                lstlatlon.Clear(); lstozone.Clear();
                return true;
            }
            catch (Exception ex)
            {
                _innermsg = "Handle SMAT input data throws an error: " + Environment.NewLine + ex.Message;
            }
            return false;
        }

        /// <summary>
        /// Start analysis
        /// </summary>
        public bool Start()
        {
            try
            {
                ok = HandleSMATInputData();
                if (!ok) return false;

                #region declare fixed params
                double krig_tol = 0.00001;
                double phi_K = 1.0 / _settings.Cmaqres;
                double phi_Q = 1.5 / _settings.Cmaqres;
                #endregion

                double[] data_d, data_d_clone, data_c1, data_c2, data_c3;
                int k = 0, len, count;
                List<Vector<double>> lstConcRow = new List<Vector<double>>();
                List<Vector<double>> lstLatlonRow = new List<Vector<double>>();

                #region spatialize model matrix
                double vct2 = MathNet.Numerics.Constants.Pi / 180;
                data_d = _matrixmodellatlon.ToColumnWiseArray();
                data_d_clone = new double[data_d.Length];
                MathNet.Numerics.Control.LinearAlgebraProvider.ScaleArray(vct2, data_d, data_d_clone);
                var s1 = Matrix<double>.Build.DenseOfColumnMajor(_matrixmodellatlon.RowCount, _matrixmodellatlon.ColumnCount, data_d_clone);

                double[] low = MatrixCompute.MinColumnWise(s1).Item1;
                double[] up = MatrixCompute.MaxColumnWise(s1).Item1;
                data_c1 = new double[20];
                data_c2 = new double[20];
                for (int i = 0; i < 20; i++)
                {
                    data_c1[i] = low[0] + ((up[0] - low[0]) / 19) * i;
                    data_c2[i] = low[1] + ((up[1] - low[1]) / 19) * i;
                }
                var qlat = Matrix<double>.Build.DenseOfColumnMajor(1, 20, data_c1);
                qlat = qlat.Transpose();
                var qlon = Matrix<double>.Build.DenseOfColumnMajor(1, 20, data_c2);
                qlon = qlon.Transpose();

                SpecialFunction.SpatializeMatrix(ref s1);

                var s2 = MatrixCompute.UniqueRow(_matrixmonitorlatlon);
                data_d = s2.ToColumnWiseArray();
                data_d_clone = new double[data_d.Length];
                MathNet.Numerics.Control.LinearAlgebraProvider.ScaleArray(vct2, data_d, data_d_clone);
                s2 = Matrix<double>.Build.DenseOfColumnMajor(s2.RowCount, s2.ColumnCount, data_d_clone);
                SpecialFunction.SpatializeMatrix(ref s2);

                data_d = _matrixmonitorlatlon.ToColumnWiseArray();
                data_d_clone = new double[data_d.Length];
                MathNet.Numerics.Control.LinearAlgebraProvider.ScaleArray(vct2, data_d, data_d_clone);
                _matrixmonitorlatlon = Matrix<double>.Build.DenseOfColumnMajor(_matrixmonitorlatlon.RowCount, _matrixmonitorlatlon.ColumnCount, data_d_clone);
                SpecialFunction.SpatializeMatrix(ref _matrixmonitorlatlon);

                data_c1 = MatrixCompute.SortColumnWise(MatrixCompute.Repmat(qlat, 20, 1)).Item1[0];
                data_c2 = MatrixCompute.Repmat(qlon, 20, 1).ToColumnWiseArray();
                len = data_c1.Length;
                data_d = new double[len * 2];
                Array.Copy(data_c1, 0, data_d, 0, len);
                Array.Copy(data_c2, 0, data_d, len, len);
                var s4 = Matrix<double>.Build.DenseOfColumnMajor(len, 2, data_d);
                SpecialFunction.SpatializeMatrix(ref s4);

                //GC.Collect();
                #endregion

                #region al least meet the condition: K(s,u)>0.1
                var s3 = SpecialFunction.ComputeMatrix_s3(s2, s1, phi_K, krig_tol);
                data_d = MatrixCompute.MaxRowWise(s3).Item1;
                len = data_d.Length;
                //lstConcRow.Clear(); lstLatlonRow.Clear();
                for (int i = 0; i < len; i++)
                {
                    if (data_d[i] > 0.1)
                    {
                        lstConcRow.Add(s2.Row(i));
                        lstLatlonRow.Add(s3.Row(i));
                    }
                }
                //len = lstConcRow.Count;
                s2 = Matrix<double>.Build.DenseOfRowVectors(lstConcRow);
                s3 = Matrix<double>.Build.SparseOfRowVectors(lstLatlonRow);

                GC.Collect();
                #endregion

                #region get the matrix "mdta"
                data_d = _matrixmonitorozone.ToColumnWiseArray();
                var tuple = MatrixCompute.MatrixRowIsmemberMatrixRow(_matrixmonitorlatlon, s2);
                count = tuple.Item1.Where((bool a) => a == true).Count();
                Debug.WriteLine("count = {0}", count);
                k = 0;
                data_c1 = new double[count]; data_c2 = new double[count]; data_d_clone = new double[count * 2];
                for (int i = 0; i < tuple.Item1.Length; i++)
                {
                    if (tuple.Item1[i])
                    {
                        data_c1[k] = data_d[i];
                        data_c2[k] = tuple.Item2[i];
                        k++;
                    }
                }
                Array.Copy(data_c1, 0, data_d_clone, 0, count);
                Array.Copy(data_c2, 0, data_d_clone, count, count);
                Matrix<double> mdta = Matrix<double>.Build.DenseOfColumnMajor(count, 2, data_d_clone);
                if (count > 1)
                {
                    data_c2 = mdta.Column(1).ToArray();
                    count = data_c2.Length;
                    int[] index_arr = MatrixCompute.SortArray(data_c2).Item2;
                    Matrix<double> item2 = Matrix<double>.Build.Dense(count, 2);
                    for (int i = 0; i < count; i++)
                    {
                        item2.SetRow(i, mdta.Row(index_arr[i] - 1));
                    }
                    data_c1 = item2.Column(0).ToArray();
                    data_c2 = item2.Column(1).ToArray();
                    double[] match = MatrixCompute.UniqueArray(data_c2);
                    int vct3 = match.Length, vct4 = item2.RowCount;
                    if (vct3 < vct4)
                    {
                        double[] temp = new double[vct3], wgt = new double[vct3];
                        List<double> mean = new List<double>();
                        for (int i = 0; i < vct3; i++)
                        {
                            mean.Clear();
                            for (int j = 0; j < vct4; j++)
                            {
                                if (match[i] == data_c2[j])
                                    mean.Add(data_c1[j]);
                            }
                            temp[i] = mean.ToArray().Average();
                            wgt[i] = 1.0;
                        }
                        data_d = new double[vct3 * 3];
                        Array.Copy(temp, 0, data_d, 0, vct3);
                        Array.Copy(match, 0, data_d, vct3, vct3);
                        Array.Copy(wgt, 0, data_d, vct3 * 2, vct3);
                        mdta = Matrix<double>.Build.DenseOfColumnMajor(vct3, 3, data_d);
                        mean.Clear();
                    }
                    else
                    {
                        data_c3 = new double[vct4];
                        for (int i = 0; i < vct4; i++)
                        {
                            data_c3[i] = 1.0;
                        }
                        data_d = new double[vct4 * 3];
                        Array.Copy(data_c1, 0, data_d, 0, vct4);
                        Array.Copy(data_c2, 0, data_d, vct4, vct4);
                        Array.Copy(data_c3, 0, data_d, vct4 * 2, vct4);
                        mdta = Matrix<double>.Build.DenseOfColumnMajor(vct4, 3, data_d);
                    }
                }

                //GC.Collect();
                #endregion

                #region compute the distance
                var tup_s4_s5 = SpecialFunction.ComputeMatrix_s4_s5(s1, s4, phi_Q, krig_tol);
                s4 = tup_s4_s5.Item1;
                var s5 = tup_s4_s5.Item2;

                GC.Collect();
                #endregion

                #region compute the parameter "phi_b0s"
                double[] probs = new double[9];
                data_c1 = mdta.Column(0).ToArray();
                data_c2 = mdta.Column(1).ToArray();
                data_c3 = mdta.Column(2).ToArray();
                lstLatlonRow.Clear();
                for (int i = 0; i < data_c2.Length; i++)
                {
                    lstLatlonRow.Add(s2.Row((int)data_c2[i] - 1));
                }
                var s2_row = Matrix<double>.Build.DenseOfRowVectors(lstLatlonRow);
                var ds = MatrixCompute.Pdist2(s2_row, s2_row);
                double[] tmp = new double[mdta.RowCount];
                double max;
                data_d = _matrixmodelozone.Column(0).ToArray();
                for (int i = 0; i < data_c2.Length; i++)
                {
                    double[] data_d2 = s3.Row((int)data_c2[i] - 1).ToArray();
                    max = data_d2.Max();
                    double v = data_d[Array.IndexOf(data_d2, max)];
                    tmp[i] = v;
                }
                LinearRegression sl = MatrixCompute.Regstats(tmp, data_c1);
                if (sl.ResidualVariance == 0.0)
                    sl.ResidualVariance = 0.01;
                double[] b0sinvmat_d;
                max = MatrixCompute.MaxColumnWise(ds).Item1.Max();
                data_d = ds.ToColumnWiseArray();
                for (int i = 0; i < 9; i++)
                {
                    double p = i * 0.1 + 0.1;
                    b0sinvmat_d = new double[data_d.Length];
                    for (int j = 0; j < data_d.Length; j++)
                    {
                        double d = -3.0 / (p * max);
                        d *= data_d[j];
                        d = Math.Pow(MathNet.Numerics.Constants.E, d);
                        b0sinvmat_d[j] = d;
                    }
                    double[] matrixEv = new double[b0sinvmat_d.Length];
                    System.Numerics.Complex[] vectorEv = new System.Numerics.Complex[ds.RowCount];
                    double[] matrixD = new double[b0sinvmat_d.Length];
                    MathNet.Numerics.Control.LinearAlgebraProvider.EigenDecomp(true, ds.RowCount, b0sinvmat_d, matrixEv, vectorEv, matrixD);
                    double v = 0.0;
                    for (int j = 0; j < ds.RowCount; j++)
                    {
                        double d = vectorEv[j].Real;
                        d = Math.Log(d, MathNet.Numerics.Constants.E);
                        v += d;
                    }
                    v *= -0.5;
                    double[] tmp1 = new double[data_c1.Length];
                    for (int j = 0; j < data_c1.Length; j++)
                    {
                        tmp1[j] = data_c1[j] - sl.Coefficients.Item1 - sl.Coefficients.Item2 * tmp[j];
                    }
                    double[] tmp2 = new double[b0sinvmat_d.Length];
                    MathNet.Numerics.Control.LinearAlgebraProvider.ScaleArray(sl.ResidualVariance, b0sinvmat_d, tmp2);
                    MathNet.Numerics.Control.LinearAlgebraProvider.LUInverse(tmp2, ds.RowCount);
                    double[] tmp3 = new double[ds.ColumnCount];
                    MathNet.Numerics.Control.LinearAlgebraProvider.MatrixMultiply(tmp1, 1, tmp1.Length, tmp2, ds.RowCount, ds.ColumnCount, tmp3);
                    double[] tmp4 = new double[1];
                    MathNet.Numerics.Control.LinearAlgebraProvider.MatrixMultiply(tmp3, 1, tmp3.Length, tmp1, tmp1.Length, 1, tmp4);
                    v -= 0.5 * tmp4[0];
                    probs[i] = v;
                }
                double phi_b0s = 3.0 / ((Array.IndexOf(probs, probs.Max()) * 0.1 + 0.1) * max);
                Debug.WriteLine("phi_b0s = {0}", phi_b0s);

                GC.Collect();
                #endregion

                #region get the matrixes for MCMC
                int order_b0sinvmat, row_vct8, column_vct8, row_krigmat, column_krigmat;
                double[] data_b0sinvmat = SpecialFunction.Compute_b0sinvmat(phi_b0s, ds, out order_b0sinvmat);
                double[] data_vct8 = SpecialFunction.Compute_vct8(phi_b0s, s1, s2_row, out row_vct8, out column_vct8);
                double[] data_krigmat = SpecialFunction.Compute_krigmat(data_vct8, row_vct8, column_vct8, data_b0sinvmat, order_b0sinvmat, out  row_krigmat, out  column_krigmat);
                double[] data_mvars = SpecialFunction.Compute_mvars(data_krigmat, row_krigmat, column_krigmat, data_vct8);

                GC.Collect();
                #endregion

                double a_s2 = 2.5;
                double b_s2 = 1.5;
                double tY2_a = a_s2;
                double tY2_b = b_s2;

                #region get the parameters for MCMC
                int n = mdta.RowCount;
                double[] b = new double[] { 0.0, 1.0 };
                lstLatlonRow.Clear();
                for (int i = 0; i < data_c2.Length; i++)
                {
                    lstLatlonRow.Add(s3.Row((int)data_c2[i] - 1));
                }
                var s3_row = Matrix<double>.Build.DenseOfRowVectors(lstLatlonRow);
                data_d = s3_row.ToColumnWiseArray();
                for (int i = 0; i < data_d.Length; i++)
                {
                    if (data_d[i] > 0.0) data_d[i] = 1.0;
                    else data_d[i] = 0.0;
                }
                var s3_row2 = Matrix<double>.Build.DenseOfColumnMajor(s3_row.RowCount, s3_row.ColumnCount, data_d);
                data_d = MatrixCompute.SumColumnWise(s3_row2);
                bool[] ind = new bool[data_d.Length];
                for (int i = 0; i < data_d.Length; i++)
                {
                    if (data_d[i] != 0.0 && (data_d[i] + 1.0 != 1.0)) ind[i] = true;
                }
                count = ind.Where((bool a) => a == true).Count();
                int row_s3_ind, column_s3_ind, row_s5_ind, column_s5_ind;
                double[] data_s3_ind = SpecialFunction.GetParamsForXt_s3_ind(ind, count, s3_row, out row_s3_ind, out column_s3_ind);
                double[] data_s5_ind = SpecialFunction.GetParamsForXt_s5_ind(ind, count, s5, out  row_s5_ind, out  column_s5_ind);
                double[] data_cmaq2_ind = SpecialFunction.GetParamsForXt_cmaq2_ind(ind, count, _matrixmodelozone);
                double[] Q = new double[s4.RowCount];
                var Xt = SpecialFunction.ComputeMatrix_XtOrXtp(data_s3_ind, row_s3_ind, column_s3_ind, data_s5_ind, row_s5_ind, column_s5_ind, data_cmaq2_ind, Q, n);
                double s2b = 1.0;
                double tY2 = SpecialFunction.Compute_tY2(data_c1, n, Xt);
                Debug.WriteLine("tY2 = {0}", tY2);
                double[] W_data = (double[])data_c3.Clone();
                double[] W = MathNet.Numerics.LinearAlgebra.Matrix<double>.Build.Diagonal(W_data).ToColumnWiseArray();

                GC.Collect();
                #endregion

                double[] Ybar = new double[row_krigmat];
                double[] Y2bar = new double[row_krigmat];
                int kp = 0;
                var bdraws = Matrix<double>.Build.Dense(_settings.Numit, 2);
                var vardraws = Matrix<double>.Build.Dense(_settings.Numit, 2);

                double vctt1 = a_s2 + (double)n / 2.0;
                double vctt2 = b_s2 * 2.0;
                double vctt3 = (tY2_a + (double)n - 2.0) / 2.0;
                double[] vctt4 = new double[] { 0.002, 0.0, 0.0, 0.002 };
                count = _settings.Numit * _settings.Thin + _settings.Burn;
                List<int> lst2 = new List<int>();
                int m = _settings.Burn + 1;
                while (m <= count)
                {
                    lst2.Add(m);
                    m += _settings.Thin;
                }
                Matrix<double> Lt, Xtp;
                double[] prop, pmean, pvar, kY, kY2;
                double alph;

                #region compute the result
                double[] data_Xt, data_Xt_trans, data_XtXt;
                int row_Xt, column_Xt, order_XtXt;
                double[] data_b0s, data_c1_b0s;
                double[] data_s4 = s4.ToColumnWiseArray();
                int row_s4 = s4.RowCount, column_s4 = s4.ColumnCount;
                double[] data_matrixModelConc_c1 = _matrixmodelozone.Column(0).ToArray();
                for (int mcmcit = 1; mcmcit <= count; mcmcit++)
                {
                    if (mcmcit % 100 == 0) Debug.WriteLine("time: {0}, mcmcit = {1}", new object[] { DateTime.Now.ToString("yyyy-MM-dd HH:mm:ss.fff"), mcmcit });
                    if (mcmcit == count) Debug.WriteLine("time: {0}, mcmcit = {1}, finished", new object[] { DateTime.Now.ToString("yyyy-MM-dd HH:mm:ss.fff"), mcmcit });

                    Lt = SpecialFunction.ComputeMatrix_Lt_1(W, tY2, data_b0sinvmat, order_b0sinvmat, s2b);

                    data_Xt = Xt.ToColumnWiseArray();
                    data_Xt_trans = Xt.ToRowWiseArray();
                    row_Xt = Xt.RowCount;
                    column_Xt = Xt.ColumnCount;
                    order_XtXt = column_Xt;
                    data_XtXt = new double[column_Xt * column_Xt];
                    MathNet.Numerics.Control.LinearAlgebraProvider.MatrixMultiply(data_Xt_trans, column_Xt, row_Xt, data_Xt, row_Xt, column_Xt, data_XtXt);

                    data_b0s = SpecialFunction.Compute_b0s(Lt, data_c1, data_Xt, row_Xt, column_Xt, b, tY2, n);
                    data_c1_b0s = new double[data_b0s.Length];
                    MathNet.Numerics.Control.LinearAlgebraProvider.SubtractArrays(data_c1, data_b0s, data_c1_b0s);

                    s2b = SpecialFunction.Compute_s2b(vctt1, data_b0s, data_b0sinvmat, order_b0sinvmat, vctt2);
                    tY2 = SpecialFunction.Compute_tY2(vctt3, tY2_b, data_c1_b0s, W, data_Xt, row_Xt, column_Xt, data_XtXt, order_XtXt, data_Xt_trans);

                    Lt = SpecialFunction.ComputeMatrix_Lt_2(data_XtXt, order_XtXt, tY2, vctt4);
                    b = SpecialFunction.Compute_b(Lt, data_Xt_trans, row_Xt, column_Xt, data_c1_b0s, tY2);

                    prop = SpecialFunction.Compute_prop(data_s4, row_s4, column_s4, Q);
                    Xtp = SpecialFunction.ComputeMatrix_XtOrXtp(data_s3_ind, row_s3_ind, column_s3_ind, data_s5_ind, row_s5_ind, column_s5_ind, data_cmaq2_ind, prop, n);
                    alph = SpecialFunction.Compute_alph(data_c1_b0s, Xtp, b, data_Xt, row_Xt, column_Xt);
                    if (Math.Log(MathNet.Numerics.Generate.Uniform(1)[0], MathNet.Numerics.Constants.E) < alph)
                    {
                        if (SpecialFunction.Judge_tY2_Xtp(tY2, Xtp))
                        {
                            Q = (double[])prop.Clone();
                            Xt = Xtp.Clone();
                        }
                    }
                    //GC.Collect();

                    if (lst2.Contains(mcmcit))
                    {
                        kp += 1;
                        pmean = SpecialFunction.Compute_pmean(b, data_matrixModelConc_c1, data_krigmat, row_krigmat, column_krigmat, data_b0s);
                        pvar = SpecialFunction.Compute_pvar(s2b, data_mvars, tY2);
                        var tup_kY = SpecialFunction.Compute_kY(pmean, pvar);
                        kY = tup_kY.Item1;
                        kY2 = tup_kY.Item2;
                        Ybar = SpecialFunction.Compute_Ybar(Ybar, kp, kY);
                        Y2bar = SpecialFunction.Compute_Y2bar(Y2bar, kp, kY2);
                        vardraws.SetRow(kp - 1, new double[] { s2b, tY2 });
                        bdraws.SetRow(kp - 1, b);
                        //GC.Collect();
                    }
                }
                lst2.Clear();
                #endregion

                #region prepare the output data for smat
                // MCMC
                count = bdraws.RowCount;
                for (int i = 0; i < count; i++)
                {
                    _mcmcoutputitem = new MCMC(bdraws.At(i, 0), bdraws.At(i, 1), vardraws.At(i, 0), vardraws.At(i, 1));
                    _mcmcoutputlist.Add(_mcmcoutputitem);
                }
                // result
                double[] result = SpecialFunction.Compute_Result(kp, Y2bar, Ybar);
                count = result.Length;
                List<OzoneModel> dsresult = new List<OzoneModel>(count);
                for (int i = 0; i < count; i++)
                {
                    _modelitem = new OzoneModel("", _matrixmodellatlon.At(i, 0), _matrixmodellatlon.At(i, 1), 0.0, Ybar[i], result[i]);
                    dsresult.Add(_modelitem);
                }
                //Output(dsresult);
                count = _resultoutputlist.Count;
                for (int i = 0; i < count; i++)
                {
                    _modelitem = _resultoutputlist[i];
                    for (int j = dsresult.Count - 1; j >= 0; j--)
                    {
                        if (dsresult[j].Lat == _modelitem.Lat && dsresult[j].Lon == _modelitem.Lon)
                        {
                            _modelitem.DS_Prediction = dsresult[j].DS_Prediction;
                            _modelitem.DS_SEpred = dsresult[j].DS_SEpred;
                            dsresult.RemoveAt(j);
                            break;
                        }
                    }
                }
                #endregion
                lstConcRow.Clear(); lstLatlonRow.Clear();
                return true;
            }
            catch (Exception ex)
            {
                _innermsg = "The ozone downscaler analysis throws an error: " + Environment.NewLine + ex.Message;
            }
            return false;
        }
        [Obsolete("DS暂时输出数据，方便与MATLAB结果对比")]
        private void Output(List<OzoneModel> result)
        {
            try
            {
                string filename, filepath;
                filename = DateTime.Now.ToString("yyyyMMddHHmmssfff") + ".csv";
                filepath = AppDomain.CurrentDomain.BaseDirectory + @"\Result\";
                if (!Directory.Exists(filepath))
                    Directory.CreateDirectory(filepath);
                filepath += filename;
                FileStream fs = new FileStream(filepath, FileMode.Create, FileAccess.Write);
                StreamWriter streamWriter = new StreamWriter(fs, Encoding.Default);
                StringBuilder sb = new StringBuilder();
                sb.Append("Latitude,Longitude,Prediction,SEpred");
                streamWriter.WriteLine(sb.ToString());
                count = result.Count;
                for (int i = 0; i < count; i++)
                {
                    _modelitem = result[i];
                    sb.Clear();
                    sb.AppendFormat("{0},{1},{2},{3}", new object[] { _modelitem.Lat, _modelitem.Lon, _modelitem.DS_Prediction, _modelitem.DS_SEpred });
                    streamWriter.WriteLine(sb.ToString());
                }
                streamWriter.Close();
                fs.Close();
                sb.Clear();
            }
            catch
            {
                
            }
        }
    }
}