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

namespace SMAT_CE.DownscaleR
{
    /// <summary>
    /// DS特殊方法
    /// </summary>
    public static class SpecialFunction
    {
        static int _maxThreadCount = MathNet.Numerics.Control.MaxDegreeOfParallelism;
        /// <summary>
        /// 矩阵空间化
        /// </summary>
        public static void SpatializeMatrix(ref Matrix<double> matrix)
        {
            double[] data_d, row, data_c1, data_c2, data_c3;
            int len = matrix.RowCount;
            data_c1 = new double[len]; data_c2 = new double[len]; data_c3 = new double[len];
            for (int i = 0; i < len; i++)
            {
                row = matrix.Row(i).ToArray();
                data_c1[i] = Math.Cos(row[0]) * Math.Cos(row[1]) * 6378.1;
                data_c2[i] = Math.Cos(row[0]) * Math.Sin(row[1]) * 6378.1;
                data_c3[i] = Math.Sin(row[0]) * 6378.1;
            }
            data_d = new double[len * 3];
            Array.Copy(data_c1, 0, data_d, 0, len);
            Array.Copy(data_c2, 0, data_d, len, len);
            Array.Copy(data_c3, 0, data_d, len * 2, len);
            matrix = Matrix<double>.Build.DenseOfColumnMajor(len, 3, data_d);
            data_d = new double[0]; row = new double[0]; data_c1 = new double[0]; data_c2 = new double[0]; data_c3 = new double[0];
        }
        /// <summary>
        /// 计算矩阵 s3；
        /// s{3} = exp(-phi_K .* pdist2(s{2}, s{1}))；
        /// s{3}(s{3} &lt; krig_tol) = 0；
        /// s{3} = sparse(s{3})；
        /// </summary>
        public static Matrix<double> ComputeMatrix_s3(Matrix<double> s2, Matrix<double> s1, double phi_K, double krig_tol, int threadCount = 0)
        {
            #region 多线程
            //threadCount = (threadCount > _maxThreadCount || threadCount <= 0) ? _maxThreadCount : threadCount;
            //if (s2.RowCount <= threadCount)
            //    threadCount = 1;

            //int row_base = s2.RowCount / threadCount;
            //int row_remain = s2.RowCount % threadCount;
            //List<Tuple<int, int>> indexList = new List<Tuple<int, int>>();
            //for (int i = 0; i < threadCount; i++)
            //{
            //    int row_count = row_base;
            //    if (i == threadCount - 1)
            //    {
            //        row_count = row_base + row_remain;
            //    }
            //    Tuple<int, int> item = new Tuple<int, int>(i * row_base, i * row_base + row_count);
            //    indexList.Add(item);
            //}

            //Task<Matrix<double>>[] tasks = new Task<Matrix<double>>[threadCount];
            //for (int i = 0; i < threadCount; i++)
            //{
            //    int from = indexList[i].Item1;
            //    int to = indexList[i].Item2;
            //    tasks[i] = Task<Matrix<double>>.Factory.StartNew(() =>
            //    {
            //        int len1 = to - from, len2 = s1.RowCount;
            //        double[] resultrow = new double[len1 * len2];
            //        for (int j = 0; j < len2; j++)
            //        {
            //            int n = 0;
            //            double[] datarow = new double[len1];
            //            double[] right = s1.Row(j).ToArray();
            //            for (int k = from; k < to; k++)
            //            {
            //                double[] left = s2.Row(k).ToArray();
            //                double d = MathNet.Numerics.Distance.Euclidean(left, right);
            //                d *= -phi_K;
            //                d = Math.Pow(Constants.E, d);
            //                if (d < krig_tol) d = 0.0;
            //                datarow[n] = d;
            //                n++;
            //            }
            //            Array.Copy(datarow, 0, resultrow, j * len1, len1);
            //        }
            //        return Matrix<double>.Build.SparseOfColumnMajor(len1, len2, resultrow);
            //    });
            //}
            //Task.WaitAll(tasks);

            //Matrix<double> result = Matrix<double>.Build.Sparse(s2.RowCount, s1.RowCount);
            //for (int i = 0; i < threadCount; i++)
            //{
            //    int from = indexList[i].Item1;
            //    var submatrix = tasks[i].Result;
            //    result.SetSubMatrix(from, 0, submatrix);
            //}
            //return result;
            #endregion
            #region 目前单线程
            //Matrix<double> result = Matrix<double>.Build.Sparse(s2.RowCount, s1.RowCount);
            //for (int i = 0; i < result.RowCount; i++)
            //{
            //    double[] left = s2.Row(i).ToArray();
            //    for (int j = 0; j < result.ColumnCount; j++)
            //    {
            //        double[] right = s1.Row(j).ToArray();
            //        double d = MathNet.Numerics.Distance.Euclidean(left, right);
            //        d *= -phi_K;
            //        d = Math.Pow(Constants.E, d);
            //        if (d < krig_tol) d = 0.0;
            //        result.At(i, j, d);
            //    }
            //}
            //return result;
            #endregion
            #region 目前单线程改进方式
            int len1 = s2.RowCount, len2 = s1.RowCount;
            double[] data = new double[len1 * len2], datarow, right, left;
            for (int j = 0; j < len2; j++)
            {
                datarow = new double[len1];
                right = s1.Row(j).ToArray();
                for (int i = 0; i < len1; i++)
                {
                    left = s2.Row(i).ToArray();
                    double d = MathNet.Numerics.Distance.Euclidean(left, right);
                    d *= -phi_K;
                    d = Math.Pow(MathNet.Numerics.Constants.E, d);
                    if (d < krig_tol) d = 0.0;
                    datarow[i] = d;
                }
                Array.Copy(datarow, 0, data, j * len1, len1);
            }
            datarow = new double[0]; right = new double[0]; left = new double[0];
            return Matrix<double>.Build.SparseOfColumnMajor(len1, len2, data);
            #endregion
        }
        /// <summary>
        /// 计算矩阵 s4 s5；
        /// s{5} = pdist2(s{1}, s{4})；
        /// s{4} = pdist2(s{4}, s{4})；
        /// s{5} = mrdivide(exp(-phi_Q .* s{5}), exp(-phi_Q .* s{4}))；
        /// s{5}(abs(s{5}) &lt; krig_tol) = 0；
        /// s{5} = sparse(s{5})；
        /// s{4} = sparse(chol(exp(-phi_Q .* s{4}), 'lower'))；
        /// </summary>
        public static Tuple<Matrix<double>, Matrix<double>> ComputeMatrix_s4_s5(Matrix<double> s1, Matrix<double> s4, double phi_Q, double krig_tol, int threadCount = 0)
        {
            #region 多线程
            //int threadCount2 = threadCount;
            //threadCount2 = (threadCount2 > _maxThreadCount || threadCount2 <= 0) ? _maxThreadCount : threadCount2;
            //if (s1.RowCount <= threadCount2)
            //    threadCount2 = 1;
            //int row_base = s1.RowCount / threadCount2;
            //int row_remain = s1.RowCount % threadCount2;

            //List<Tuple<int, int>> indexList = new List<Tuple<int, int>>();
            //for (int i = 0; i < threadCount2; i++)
            //{
            //    int row_count = row_base;
            //    if (i == threadCount2 - 1)
            //    {
            //        row_count = row_base + row_remain;
            //    }
            //    Tuple<int, int> item = new Tuple<int, int>(i * row_base, i * row_base + row_count);
            //    indexList.Add(item);
            //}

            //Task<Matrix<double>>[] tasks = new Task<Matrix<double>>[threadCount2];
            //for (int i = 0; i < threadCount2; i++)
            //{
            //    int from = indexList[i].Item1;
            //    int to = indexList[i].Item2;
            //    tasks[i] = Task<Matrix<double>>.Factory.StartNew(() =>
            //    {
            //        int len1 = to - from, len2 = s4.RowCount;
            //        double[] resultrow = new double[len1 * len2];
            //        for (int j = 0; j < len2; j++)
            //        {
            //            int n = 0;
            //            double[] datarow = new double[len1];
            //            double[] right = s4.Row(j).ToArray();
            //            for (int k = from; k < to; k++)
            //            {
            //                double[] left = s1.Row(k).ToArray();
            //                double d = MathNet.Numerics.Distance.Euclidean(left, right);
            //                d *= -phi_Q;
            //                d = Math.Pow(Constants.E, d);
            //                datarow[n] = d;
            //                n++;
            //            }
            //            Array.Copy(datarow, 0, resultrow, j * len1, len1);
            //        }
            //        return Matrix<double>.Build.SparseOfColumnMajor(len1, len2, resultrow);
            //    });
            //}
            //Task.WaitAll(tasks);

            //Matrix<double> s5_exp = Matrix<double>.Build.Sparse(s1.RowCount, s4.RowCount);
            //for (int i = 0; i < threadCount2; i++)
            //{
            //    int from = indexList[i].Item1;
            //    var submatrix = tasks[i].Result;
            //    s5_exp.SetSubMatrix(from, 0, submatrix);
            //}
            //double[] data_s5_exp = s5_exp.ToColumnWiseArray();

            //threadCount2 = threadCount;
            //threadCount2 = (threadCount2 > _maxThreadCount || threadCount2 <= 0) ? _maxThreadCount : threadCount2;
            //if (s4.RowCount <= threadCount2)
            //    threadCount2 = 1;
            //row_base = s4.RowCount / threadCount2;
            //row_remain = s4.RowCount % threadCount2;

            //indexList.Clear();
            //for (int i = 0; i < threadCount2; i++)
            //{
            //    int row_count = row_base;
            //    if (i == threadCount2 - 1)
            //    {
            //        row_count = row_base + row_remain;
            //    }
            //    Tuple<int, int> item = new Tuple<int, int>(i * row_base, i * row_base + row_count);
            //    indexList.Add(item);
            //}

            //tasks = new Task<Matrix<double>>[threadCount2];
            //for (int i = 0; i < threadCount2; i++)
            //{
            //    int from = indexList[i].Item1;
            //    int to = indexList[i].Item2;
            //    tasks[i] = Task<Matrix<double>>.Factory.StartNew(() =>
            //    {
            //        int len1 = to - from, len2 = s4.RowCount;
            //        double[] resultrow = new double[len1 * len2];
            //        for (int j = 0; j < len2; j++)
            //        {
            //            int n = 0;
            //            double[] datarow = new double[len1];
            //            double[] right = s4.Row(j).ToArray();
            //            for (int k = from; k < to; k++)
            //            {
            //                double[] left = s4.Row(k).ToArray();
            //                double d = MathNet.Numerics.Distance.Euclidean(left, right);
            //                d *= -phi_Q;
            //                d = Math.Pow(Constants.E, d);
            //                datarow[n] = d;
            //                n++;
            //            }
            //            Array.Copy(datarow, 0, resultrow, j * len1, len1);
            //        }
            //        return Matrix<double>.Build.SparseOfColumnMajor(len1, len2, resultrow);
            //    });
            //}
            //Task.WaitAll(tasks);

            //Matrix<double> s4_exp = Matrix<double>.Build.Sparse(s4.RowCount, s4.RowCount);
            //for (int i = 0; i < threadCount2; i++)
            //{
            //    int from = indexList[i].Item1;
            //    var submatrix = tasks[i].Result;
            //    s4_exp.SetSubMatrix(from, 0, submatrix);
            //}
            //double[] data_s4_exp = s4_exp.ToColumnWiseArray();

            //Task<Matrix<double>> t1 = Task<Matrix<double>>.Factory.StartNew(() =>
            //{
            //    double[] data = (double[])data_s4_exp.Clone();
            //    MathNet.Numerics.Control.LinearAlgebraProvider.LUInverse(data, s4_exp.RowCount);
            //    double[] result = new double[s5_exp.RowCount * s4_exp.RowCount];
            //    MathNet.Numerics.Control.LinearAlgebraProvider.MatrixMultiply(data_s5_exp, s5_exp.RowCount, s5_exp.ColumnCount, data, s4_exp.RowCount, s4_exp.ColumnCount, result);
            //    for (int i = 0; i < result.Length; i++)
            //    {
            //        if (Math.Abs(result[i]) < krig_tol)
            //            result[i] = 0.0;
            //    }
            //    return Matrix<double>.Build.SparseOfColumnMajor(s5_exp.RowCount, s4_exp.RowCount, result);
            //});
            //Task<Matrix<double>> t2 = Task<Matrix<double>>.Factory.StartNew(() =>
            //{
            //    double[] data = (double[])data_s4_exp.Clone();
            //    MathNet.Numerics.Control.LinearAlgebraProvider.CholeskyFactor(data, s4_exp.RowCount);
            //    return Matrix<double>.Build.SparseOfColumnMajor(s4_exp.RowCount, s4_exp.RowCount, data);
            //});
            //Task.WaitAll(t1, t2);
            //return Tuple.Create<Matrix<double>, Matrix<double>>(t2.Result, t1.Result);
            #endregion
            #region 目前单线程
            //Matrix<double> s5_exp = Matrix<double>.Build.Sparse(s1.RowCount, s4.RowCount);
            //for (int i = 0; i < s5_exp.RowCount; i++)
            //{
            //    double[] left = s1.Row(i).ToArray();
            //    for (int j = 0; j < s5_exp.ColumnCount; j++)
            //    {
            //        double[] right = s4.Row(j).ToArray();
            //        double d = MathNet.Numerics.Distance.Euclidean(left, right);
            //        d *= -phi_Q;
            //        d = Math.Pow(Constants.E, d);
            //        s5_exp.At(i, j, d);
            //    }
            //}
            //double[] data_s5_exp = s5_exp.ToColumnWiseArray();
            //Matrix<double> s4_exp = Matrix<double>.Build.Sparse(s4.RowCount, s4.RowCount);
            //for (int i = 0; i < s4_exp.RowCount; i++)
            //{
            //    double[] left = s4.Row(i).ToArray();
            //    for (int j = 0; j < s4_exp.ColumnCount; j++)
            //    {
            //        double[] right = s4.Row(j).ToArray();
            //        double d = MathNet.Numerics.Distance.Euclidean(left, right);
            //        d *= -phi_Q;
            //        d = Math.Pow(Constants.E, d);
            //        s4_exp.At(i, j, d);
            //    }
            //}
            //double[] data_s4_exp = s4_exp.ToColumnWiseArray();

            //double[] data = (double[])data_s4_exp.Clone();
            //MathNet.Numerics.Control.LinearAlgebraProvider.LUInverse(data, s4_exp.RowCount);
            //double[] result = new double[s5_exp.RowCount * s4_exp.RowCount];
            //MathNet.Numerics.Control.LinearAlgebraProvider.MatrixMultiply(data_s5_exp, s5_exp.RowCount, s5_exp.ColumnCount, data, s4_exp.RowCount, s4_exp.ColumnCount, result);
            //for (int i = 0; i < result.Length; i++)
            //{
            //    if (Math.Abs(result[i]) < krig_tol)
            //        result[i] = 0.0;
            //}
            //Matrix<double> s5 = Matrix<double>.Build.SparseOfColumnMajor(s5_exp.RowCount, s4_exp.RowCount, result);

            //data = (double[])data_s4_exp.Clone();
            //MathNet.Numerics.Control.LinearAlgebraProvider.CholeskyFactor(data, s4_exp.RowCount);
            //Matrix<double> s44 = Matrix<double>.Build.SparseOfColumnMajor(s4_exp.RowCount, s4_exp.RowCount, data);

            //return Tuple.Create<Matrix<double>, Matrix<double>>(s44, s5);
            #endregion
            #region 目前单线程改进方式
            int len1 = s1.RowCount, len2 = s4.RowCount;
            double[] data = new double[len1 * len2], datarow, right, left;
            for (int j = 0; j < len2; j++)
            {
                datarow = new double[len1];
                right = s4.Row(j).ToArray();
                for (int i = 0; i < len1; i++)
                {
                    left = s1.Row(i).ToArray();
                    double d = MathNet.Numerics.Distance.Euclidean(left, right);
                    d *= -phi_Q;
                    d = Math.Pow(MathNet.Numerics.Constants.E, d);
                    datarow[i] = d;
                }
                Array.Copy(datarow, 0, data, j * len1, len1);
            }
            Matrix<double> s5_exp = Matrix<double>.Build.SparseOfColumnMajor(len1, len2, data);
            data = new double[0]; datarow = new double[0]; left = new double[0]; right = new double[0];

            len1 = s4.RowCount;
            len2 = s4.RowCount;
            data = new double[len1 * len2];
            for (int j = 0; j < len2; j++)
            {
                datarow = new double[len1];
                right = s4.Row(j).ToArray();
                for (int i = 0; i < len1; i++)
                {
                    left = s4.Row(i).ToArray();
                    double d = MathNet.Numerics.Distance.Euclidean(left, right);
                    d *= -phi_Q;
                    d = Math.Pow(MathNet.Numerics.Constants.E, d);
                    datarow[i] = d;
                }
                Array.Copy(datarow, 0, data, j * len1, len1);
            }
            Matrix<double> s4_exp = Matrix<double>.Build.SparseOfColumnMajor(len1, len2, data);
            data = new double[0]; datarow = new double[0]; left = new double[0]; right = new double[0];

            double[] data_s5_exp = s5_exp.ToColumnWiseArray();
            double[] data_s4_exp = s4_exp.ToColumnWiseArray();
            data = (double[])data_s4_exp.Clone();
            MathNet.Numerics.Control.LinearAlgebraProvider.LUInverse(data, s4_exp.RowCount);
            double[] result = new double[s5_exp.RowCount * s4_exp.RowCount];
            MathNet.Numerics.Control.LinearAlgebraProvider.MatrixMultiply(data_s5_exp, s5_exp.RowCount, s5_exp.ColumnCount, data, s4_exp.RowCount, s4_exp.ColumnCount, result);
            for (int i = 0; i < result.Length; i++)
            {
                if (Math.Abs(result[i]) < krig_tol)
                    result[i] = 0.0;
            }
            Matrix<double> s5 = Matrix<double>.Build.SparseOfColumnMajor(s5_exp.RowCount, s4_exp.RowCount, result);
            data_s5_exp = new double[0]; data = new double[0]; result = new double[0];

            data = (double[])data_s4_exp.Clone();
            MathNet.Numerics.Control.LinearAlgebraProvider.CholeskyFactor(data, s4_exp.RowCount);
            Matrix<double> s44 = Matrix<double>.Build.SparseOfColumnMajor(s4_exp.RowCount, s4_exp.RowCount, data);
            data_s4_exp = new double[0]; data = new double[0];

            return Tuple.Create<Matrix<double>, Matrix<double>>(s44, s5);
            #endregion
        }
        /// <summary>
        /// 按顺序返回参数矩阵 b0sinvmat，krigmat，mvars；
        /// 计算公式：b0sinvmat = exp(-phi_b0s .* pdist2(s{2}(mdta{t}(:,2),:), s{2}(mdta{t}(:,2),:)))；(b0sinvmat = exp(-phi_b0s .* ds);)；
        /// b0sinvmat = mldivide(b0sinvmat, eye(size(b0sinvmat, 1)))；
        /// vct8 = exp(-phi_b0s .* pdist2(s{1}, s{2}(mdta{t}(:,2),:)))；(vct8 = exp(-phi_b0s .* pdist2(s{1}, s2_row))；
        /// krigmat = vct8；
        /// krigmat = krigmat * b0sinvmat；
        /// mvars = sum(krigmat .* vct8, 2)；
        /// </summary>
        [Obsolete("以前的版本")]
        public static Tuple<Matrix<double>, Matrix<double>, Matrix<double>> GetMatrixsFromPriors(double phi_b0s, Matrix<double> ds, Matrix<double> s1, Matrix<double> s2_row, int threadCount = 0)
        {
            int len1 = s1.RowCount, len2 = s2_row.RowCount;
            double[] data = new double[len1 * len2], datarow;
            for (int j = 0; j < len2; j++)
            {
                datarow = new double[len1];
                double[] right = s2_row.Row(j).ToArray();
                for (int i = 0; i < len1; i++)
                {
                    double[] left = s1.Row(i).ToArray();
                    double d = MathNet.Numerics.Distance.Euclidean(left, right);
                    d *= -phi_b0s;
                    d = Math.Pow(MathNet.Numerics.Constants.E, d);
                    datarow[i] = d;
                }
                Array.Copy(datarow, 0, data, j * len1, len1);
            }
            Matrix<double> vct8 = Matrix<double>.Build.DenseOfColumnMajor(len1, len2, data);
            double[] data_vct8 = vct8.ToColumnWiseArray();

            data = ds.ToColumnWiseArray();
            len1 = data.Length;
            double[] data_b0sinvmat = new double[len1];
            for (int i = 0; i < len1; i++)
            {
                double d = -phi_b0s * data[i];
                d = Math.Pow(MathNet.Numerics.Constants.E, d);
                data_b0sinvmat[i] = d;
            }
            MathNet.Numerics.Control.LinearAlgebraProvider.LUInverse(data_b0sinvmat, ds.RowCount);
            Matrix<double> b0sinvmat = Matrix<double>.Build.DenseOfColumnMajor(ds.RowCount, ds.ColumnCount, data_b0sinvmat);

            double[] data_krigmat = new double[vct8.RowCount * ds.ColumnCount];
            MathNet.Numerics.Control.LinearAlgebraProvider.MatrixMultiply(data_vct8, vct8.RowCount, vct8.ColumnCount, data_b0sinvmat, ds.RowCount, ds.ColumnCount, data_krigmat);
            Matrix<double> krigmat = Matrix<double>.Build.DenseOfColumnMajor(vct8.RowCount, ds.ColumnCount, data_krigmat);

            double[] data_mvars = new double[data_krigmat.Length];
            MathNet.Numerics.Control.LinearAlgebraProvider.PointWiseMultiplyArrays(data_krigmat, data_vct8, data_mvars);
            Matrix<double> mvars = Matrix<double>.Build.DenseOfColumnMajor(krigmat.RowCount, krigmat.ColumnCount, data_mvars);
            data_mvars = MatrixCompute.SumRowWise(mvars);
            mvars = Matrix<double>.Build.DenseOfColumnMajor(data_mvars.Length, 1, data_mvars);

            return Tuple.Create<Matrix<double>, Matrix<double>, Matrix<double>>(b0sinvmat, krigmat, mvars);
        }
        /// <summary>
        /// 计算参数 b0sinvmat；
        /// 计算公式：b0sinvmat = exp(-phi_b0s .* pdist2(s{2}(mdta{t}(:,2),:), s{2}(mdta{t}(:,2),:)))；(b0sinvmat = exp(-phi_b0s .* ds);)；
        /// b0sinvmat = mldivide(b0sinvmat, eye(size(b0sinvmat, 1)))；
        /// </summary>
        public static double[] Compute_b0sinvmat(double phi_b0s, Matrix<double> ds, out int order_b0sinvmat)
        {
            double[] data_1 = ds.ToColumnWiseArray();
            int len = data_1.Length;
            double[] result = new double[len];
            for (int i = 0; i < len; i++)
            {
                double d = -phi_b0s * data_1[i];
                d = Math.Pow(MathNet.Numerics.Constants.E, d);
                result[i] = d;
            }
            MathNet.Numerics.Control.LinearAlgebraProvider.LUInverse(result, ds.RowCount);
            order_b0sinvmat = ds.RowCount;
            data_1 = new double[0];
            return result;
        }
        /// <summary>
        /// 计算参数 vct8；
        /// 计算公式：vct8 = exp(-phi_b0s .* pdist2(s{1}, s{2}(mdta{t}(:,2),:)))；(vct8 = exp(-phi_b0s .* pdist2(s{1}, s2_row))；
        /// </summary>
        public static double[] Compute_vct8(double phi_b0s, Matrix<double> s1, Matrix<double> s2_row, out int row_vct8, out int column_vct8)
        {
            row_vct8 = s1.RowCount;
            column_vct8 = s2_row.RowCount;
            double[] data = new double[row_vct8 * column_vct8], datarow, right, left;
            for (int j = 0; j < column_vct8; j++)
            {
                datarow = new double[row_vct8];
                right = s2_row.Row(j).ToArray();
                for (int i = 0; i < row_vct8; i++)
                {
                    left = s1.Row(i).ToArray();
                    double d = MathNet.Numerics.Distance.Euclidean(left, right);
                    d *= -phi_b0s;
                    d = Math.Pow(MathNet.Numerics.Constants.E, d);
                    datarow[i] = d;
                }
                Array.Copy(datarow, 0, data, j * row_vct8, row_vct8);
            }
            datarow = new double[0]; left = new double[0]; right = new double[0];
            return data;
        }
        /// <summary>
        /// 计算参数 krigmat；
        /// 计算公式：krigmat = vct8；
        /// krigmat = krigmat * b0sinvmat；
        /// </summary>
        public static double[] Compute_krigmat(double[] data_vct8, int row_vct8, int column_vct8, double[] data_b0sinvmat, int order_b0sinvmat, out int row_krigmat, out int column_krigmat)
        {
            row_krigmat = row_vct8;
            column_krigmat = order_b0sinvmat;
            double[] result = new double[row_vct8 * order_b0sinvmat];
            MathNet.Numerics.Control.LinearAlgebraProvider.MatrixMultiply(data_vct8, row_vct8, column_vct8, data_b0sinvmat, order_b0sinvmat, order_b0sinvmat, result);
            return result;
        }
        /// <summary>
        /// 计算参数 mvars；
        /// 计算公式：mvars = sum(krigmat .* vct8, 2)；
        /// </summary>
        public static double[] Compute_mvars(double[] data_krigmat, int row_krigmat, int column_krigmat, double[] data_vct8)
        {
            double[] result = new double[data_krigmat.Length];
            MathNet.Numerics.Control.LinearAlgebraProvider.PointWiseMultiplyArrays(data_krigmat, data_vct8, result);
            Matrix<double> mvars = Matrix<double>.Build.DenseOfColumnMajor(row_krigmat, column_krigmat, result);
            result = MatrixCompute.SumRowWise(mvars);
            mvars = null;
            return result;
        }
        /// <summary>
        /// 按顺序返回计算 Xt 的矩阵参数 s3_ind，s5_ind，cmaq2_ind；
        /// 计算公式：s{3}(mdta{t}(:,2),ind)；
        /// s{5}(ind,:)；
        /// cmaq{2}(ind,1)；
        /// </summary>
        [Obsolete("以前的版本")]
        public static Tuple<Matrix<double>, Matrix<double>, Matrix<double>> GetParamsForXt(bool[] ind, Matrix<double> s3_row, Matrix<double> s5, Matrix<double> matrixModelConc)
        {
            //int count = ind.Where((bool a) => a == true).Count();
            int count = 0;
            for (int i = 0; i < ind.Length; i++)
                if (ind[i]) count++;

            double[] data = new double[s3_row.RowCount * count];
            for (int i = 0, j = 0; i < s3_row.ColumnCount; i++)
            {
                if (ind[i])
                {
                    double[] col = s3_row.Column(i).ToArray();
                    Array.Copy(col, 0, data, j * s3_row.RowCount, s3_row.RowCount);
                    j++;
                }
            }
            var s3_ind = Matrix<double>.Build.DenseOfColumnMajor(s3_row.RowCount, count, data);

            List<Vector<double>> lst = new List<Vector<double>>();
            for (int i = 0; i < s5.RowCount; i++)
            {
                if (ind[i]) lst.Add(s5.Row(i));
            }
            var s5_ind = Matrix<double>.Build.DenseOfRowVectors(lst);

            data = matrixModelConc.Column(0).ToArray();
            double[] result = new double[count];
            for (int i = 0, j = 0; i < data.Length; i++)
            {
                if (ind[i])
                {
                    result[j] = data[i];
                    j++;
                }
            }
            var cmaq2_ind = Matrix<double>.Build.DenseOfColumnMajor(count, 1, result);
            return Tuple.Create<Matrix<double>, Matrix<double>, Matrix<double>>(s3_ind, s5_ind, cmaq2_ind);
        }
        /// <summary>
        /// 返回计算 Xt 的矩阵参数 s3_ind；
        /// 计算公式：s{3}(mdta{t}(:,2),ind)；
        /// </summary>
        public static double[] GetParamsForXt_s3_ind(bool[] ind, int count, Matrix<double> s3_row, out int row_s3_ind, out int column_s3_ind)
        {
            double[] data = new double[s3_row.RowCount * count];
            double[] col;
            for (int i = 0, j = 0; i < s3_row.ColumnCount; i++)
            {
                if (ind[i])
                {
                    col = s3_row.Column(i).ToArray();
                    Array.Copy(col, 0, data, j * s3_row.RowCount, s3_row.RowCount);
                    j++;
                }
            }
            row_s3_ind = s3_row.RowCount;
            column_s3_ind = count;
            col = new double[0];
            return data;
        }
        /// <summary>
        /// 返回计算 Xt 的矩阵参数 s5_ind；
        /// 计算公式：s{5}(ind,:)；
        /// </summary>
        public static double[] GetParamsForXt_s5_ind(bool[] ind, int count, Matrix<double> s5, out int row_s5_ind, out int column_s5_ind)
        {
            List<Vector<double>> lst = new List<Vector<double>>();
            for (int i = 0; i < s5.RowCount; i++)
            {
                if (ind[i]) lst.Add(s5.Row(i));
            }
            row_s5_ind = count;
            column_s5_ind = s5.ColumnCount;
            var s5_ind = Matrix<double>.Build.DenseOfRowVectors(lst);
            lst.Clear();
            return s5_ind.ToColumnWiseArray();
        }
        /// <summary>
        /// 返回计算 Xt 的矩阵参数 cmaq2_ind；
        /// 计算公式：cmaq{2}(ind,1)；
        /// </summary>
        public static double[] GetParamsForXt_cmaq2_ind(bool[] ind, int count, Matrix<double> matrixModelConc)
        {
            double[] data_1 = matrixModelConc.Column(0).ToArray();
            double[] result = new double[count];
            for (int i = 0, j = 0; i < data_1.Length; i++)
            {
                if (ind[i])
                {
                    result[j] = data_1[i];
                    j++;
                }
            }
            data_1 = new double[0];
            return result;
        }
        /// <summary>
        /// 计算矩阵 Xt 或 Xtp；
        /// 计算公式：Xt = [ones(n,1) rdivide(s{3}(mdta{t}(:,2),ind) * (exp(s{5}(ind,:) * Q) .* cmaq{2}(ind,1)), s{3}(mdta{t}(:,2),ind) * exp(s{5}(ind,:) * Q))]；
        /// Xtp = [ones(n,1) rdivide(s{3}(mdta{t}(:,2),ind) * (exp(s{5}(ind,:) * prop) .* cmaq{2}(ind,1)), s{3}(mdta{t}(:,2),ind) * exp(s{5}(ind,:) * prop))]；
        /// </summary>
        [Obsolete("以前的版本")]
        public static Matrix<double> ComputeMatrix_XtOrXtp(Matrix<double> s3_ind, Matrix<double> s5_ind, Matrix<double> cmaq2_ind, double[] Q, int n)
        {
            // exp(s{5}(ind,:) * Q)
            double[] data_s5_ind = s5_ind.ToColumnWiseArray();
            double[] result1 = new double[s5_ind.RowCount];
            MathNet.Numerics.Control.LinearAlgebraProvider.MatrixMultiply(data_s5_ind, s5_ind.RowCount, s5_ind.ColumnCount, Q, Q.Length, 1, result1);
            for (int i = 0; i < result1.Length; i++)
            {
                result1[i] = Math.Pow(MathNet.Numerics.Constants.E, result1[i]);
            }
            // exp(s{5}(ind,:) * Q) .* cmaq{2}(ind,1)
            double[] data_cmaq2_ind = cmaq2_ind.ToColumnWiseArray();
            double[] result2 = new double[data_cmaq2_ind.Length];
            MathNet.Numerics.Control.LinearAlgebraProvider.PointWiseMultiplyArrays(result1, data_cmaq2_ind, result2);
            // s{3}(mdta{t}(:,2),ind) * exp(s{5}(ind,:) * Q)
            double[] data_s3_ind = s3_ind.ToColumnWiseArray();
            double[] result3 = new double[s3_ind.RowCount];
            MathNet.Numerics.Control.LinearAlgebraProvider.MatrixMultiply(data_s3_ind, s3_ind.RowCount, s3_ind.ColumnCount, result1, result1.Length, 1, result3);
            // s{3}(mdta{t}(:,2),ind) * (exp(s{5}(ind,:) * Q) .* cmaq{2}(ind,1))
            double[] result4 = new double[s3_ind.RowCount];
            MathNet.Numerics.Control.LinearAlgebraProvider.MatrixMultiply(data_s3_ind, s3_ind.RowCount, s3_ind.ColumnCount, result2, result2.Length, 1, result4);

            double[] data_xt_c2 = new double[s3_ind.RowCount];
            MathNet.Numerics.Control.LinearAlgebraProvider.PointWiseDivideArrays(result4, result3, data_xt_c2);
            double[] data_xt_c1 = new double[n];
            for (int i = 0; i < n; i++)
            {
                data_xt_c1[i] = 1.0;
            }
            double[] data_xt = new double[data_xt_c2.Length * 2];
            Array.Copy(data_xt_c1, 0, data_xt, 0, n);
            Array.Copy(data_xt_c2, 0, data_xt, data_xt_c2.Length, data_xt_c2.Length);
            return Matrix<double>.Build.DenseOfColumnMajor(data_xt_c2.Length, 2, data_xt);
        }
        /// <summary>
        /// 计算矩阵 Xt 或 Xtp；
        /// 计算公式：Xt = [ones(n,1) rdivide(s{3}(mdta{t}(:,2),ind) * (exp(s{5}(ind,:) * Q) .* cmaq{2}(ind,1)), s{3}(mdta{t}(:,2),ind) * exp(s{5}(ind,:) * Q))]；
        /// Xtp = [ones(n,1) rdivide(s{3}(mdta{t}(:,2),ind) * (exp(s{5}(ind,:) * prop) .* cmaq{2}(ind,1)), s{3}(mdta{t}(:,2),ind) * exp(s{5}(ind,:) * prop))]；
        /// </summary>
        public static Matrix<double> ComputeMatrix_XtOrXtp(double[] data_s3_ind, int row_s3_ind, int column_s3_ind, double[] data_s5_ind, int row_s5_ind, int column_s5_ind, double[] data_cmaq2_ind, double[] Q, int n)
        {
            // exp(s{5}(ind,:) * Q)
            double[] data_1 = new double[row_s5_ind];
            MathNet.Numerics.Control.LinearAlgebraProvider.MatrixMultiply(data_s5_ind, row_s5_ind, column_s5_ind, Q, Q.Length, 1, data_1);
            for (int i = 0; i < row_s5_ind; i++)
            {
                data_1[i] = Math.Pow(MathNet.Numerics.Constants.E, data_1[i]);
            }

            // exp(s{5}(ind,:) * Q) .* cmaq{2}(ind,1)
            double[] data_2 = new double[row_s5_ind];
            MathNet.Numerics.Control.LinearAlgebraProvider.PointWiseMultiplyArrays(data_1, data_cmaq2_ind, data_2);

            // s{3}(mdta{t}(:,2),ind) * exp(s{5}(ind,:) * Q)
            double[] data_3 = new double[row_s3_ind];
            MathNet.Numerics.Control.LinearAlgebraProvider.MatrixMultiply(data_s3_ind, row_s3_ind, column_s3_ind, data_1, data_1.Length, 1, data_3);

            // s{3}(mdta{t}(:,2),ind) * (exp(s{5}(ind,:) * Q) .* cmaq{2}(ind,1))
            double[] data_4 = new double[row_s3_ind];
            MathNet.Numerics.Control.LinearAlgebraProvider.MatrixMultiply(data_s3_ind, row_s3_ind, column_s3_ind, data_2, data_2.Length, 1, data_4);

            double[] data_xt_c2 = new double[row_s3_ind];
            MathNet.Numerics.Control.LinearAlgebraProvider.PointWiseDivideArrays(data_4, data_3, data_xt_c2);
            double[] data_xt_c1 = new double[n];
            for (int i = 0; i < n; i++)
                data_xt_c1[i] = 1.0;
            double[] data_xt = new double[row_s3_ind * 2];
            Array.Copy(data_xt_c1, 0, data_xt, 0, n);
            Array.Copy(data_xt_c2, 0, data_xt, row_s3_ind, row_s3_ind);

            data_1 = new double[0]; data_2 = new double[0]; data_3 = new double[0]; data_4 = new double[0]; data_xt_c1 = new double[0]; data_xt_c2 = new double[0];
            return Matrix<double>.Build.DenseOfColumnMajor(row_s3_ind, 2, data_xt);
        }
        /// <summary>
        /// 计算参数 tY2；
        /// 计算公式：tY2 = max(0.1, (mdta{t}(:,1)' * (eye(n) - mrdivide(Xt, Xt' * Xt) * Xt') * mdta{t}(:,1)) / (length(mdta{t}(:,1)) - 1))；
        /// </summary>
        public static double Compute_tY2(double[] data_c1, int n, Matrix<double> Xt)
        {
            // eye(n)
            double[] data_1 = Matrix<double>.Build.Diagonal(n, n, 1.0).ToColumnWiseArray();

            // Xt'
            double[] data_xt_trans = Xt.ToRowWiseArray();

            // (eye(n) - mrdivide(Xt, Xt' * Xt) * Xt')
            double[] data_xt = Xt.ToColumnWiseArray();
            double[] data_xt_1 = new double[Xt.ColumnCount * Xt.ColumnCount];
            MathNet.Numerics.Control.LinearAlgebraProvider.MatrixMultiply(data_xt_trans, Xt.ColumnCount, Xt.RowCount, data_xt, Xt.RowCount, Xt.ColumnCount, data_xt_1);
            MathNet.Numerics.Control.LinearAlgebraProvider.LUInverse(data_xt_1, Xt.ColumnCount);
            double[] data_xt_2 = new double[data_xt.Length];
            MathNet.Numerics.Control.LinearAlgebraProvider.MatrixMultiply(data_xt, Xt.RowCount, Xt.ColumnCount, data_xt_1, Xt.ColumnCount, Xt.ColumnCount, data_xt_2);
            double[] data_xt_3 = new double[Xt.RowCount * Xt.RowCount];
            MathNet.Numerics.Control.LinearAlgebraProvider.MatrixMultiply(data_xt_2, Xt.RowCount, Xt.ColumnCount, data_xt_trans, Xt.ColumnCount, Xt.RowCount, data_xt_3);
            double[] data_3 = new double[data_1.Length];
            MathNet.Numerics.Control.LinearAlgebraProvider.SubtractArrays(data_1, data_xt_3, data_3);

            // mdta{t}(:,1)' * (eye(n) - mrdivide(Xt, Xt' * Xt) * Xt') * mdta{t}(:,1)
            double[] data_4 = new double[n];
            MathNet.Numerics.Control.LinearAlgebraProvider.MatrixMultiply(data_c1, 1, data_c1.Length, data_3, n, n, data_4);
            double[] data_5 = new double[1];
            MathNet.Numerics.Control.LinearAlgebraProvider.MatrixMultiply(data_4, 1, n, data_c1, n, 1, data_5);

            double tY2 = data_5[0] / (data_c1.Length - 1);

            data_1 = new double[0]; data_xt_trans = new double[0]; data_xt = new double[0]; data_xt_1 = new double[0];
            data_xt_2 = new double[0]; data_xt_3 = new double[0]; data_3 = new double[0]; data_4 = new double[0];
            return Math.Max(0.1, tY2);
        }
        /// <summary>
        /// 计算矩阵 Lt；
        /// 计算公式：Lt = chol(W ./ tY2 + b0sinvmat ./ s2b)；
        /// </summary>
        [Obsolete("以前的版本")]
        public static Matrix<double> ComputeMatrix_Lt_1(double[] W, double tY2, Matrix<double> b0sinvmat, double s2b)
        {
            // W ./ tY2
            double[] result1 = new double[W.Length];
            MathNet.Numerics.Control.LinearAlgebraProvider.ScaleArray(1.0 / tY2, W, result1);

            // b0sinvmat ./ s2b
            double[] data_b0sinvmat = b0sinvmat.ToColumnWiseArray();
            double[] result2 = new double[data_b0sinvmat.Length];
            MathNet.Numerics.Control.LinearAlgebraProvider.ScaleArray(1.0 / s2b, data_b0sinvmat, result2);

            double[] result = new double[W.Length];
            MathNet.Numerics.Control.LinearAlgebraProvider.AddArrays(result1, result2, result);
            MathNet.Numerics.Control.LinearAlgebraProvider.CholeskyFactor(result, b0sinvmat.RowCount);
            var Lt = Matrix<double>.Build.DenseOfColumnMajor(b0sinvmat.RowCount, b0sinvmat.RowCount, result);
            return Lt.Transpose();
        }
        /// <summary>
        /// 计算矩阵 Lt；
        /// 计算公式：Lt = chol(W ./ tY2 + b0sinvmat ./ s2b)；
        /// </summary>
        public static Matrix<double> ComputeMatrix_Lt_1(double[] W, double tY2, double[] data_b0sinvmat, int order_b0sinvmat, double s2b)
        {
            // W ./ tY2
            double alpha = 1.0 / tY2;
            double[] data_1 = new double[W.Length];
            MathNet.Numerics.Control.LinearAlgebraProvider.ScaleArray(alpha, W, data_1);

            // b0sinvmat ./ s2b
            alpha = 1.0 / s2b;
            double[] data_2 = new double[data_b0sinvmat.Length];
            MathNet.Numerics.Control.LinearAlgebraProvider.ScaleArray(alpha, data_b0sinvmat, data_2);

            // chol(W ./ tY2 + b0sinvmat ./ s2b)
            double[] data_3 = new double[W.Length];
            MathNet.Numerics.Control.LinearAlgebraProvider.AddArrays(data_1, data_2, data_3);
            MathNet.Numerics.Control.LinearAlgebraProvider.CholeskyFactor(data_3, order_b0sinvmat);

            var Lt = Matrix<double>.Build.DenseOfColumnMajor(order_b0sinvmat, order_b0sinvmat, data_3);
            data_1 = new double[0]; data_2 = new double[0]; data_3 = new double[0];
            return Lt.Transpose();
        }
        /// <summary>
        /// 计算矩阵 b0s；
        /// 计算公式：b0s= mldivide(Lt, mldivide(Lt', (mdta{t}(:,1) - Xt * b) ./ tY2)) + mldivide(Lt, randn(n,1))；
        /// </summary>
        [Obsolete("以前的版本")]
        public static Matrix<double> ComputeMatrix_b0s(Matrix<double> Lt, double[] data_c1, Matrix<double> Xt, double[] b, double tY2, int n)
        {
            // inv(Lt)
            double[] result1 = Lt.ToColumnWiseArray();
            MathNet.Numerics.Control.LinearAlgebraProvider.LUInverse(result1, Lt.RowCount);
            // mldivide(Lt, randn(n,1))
            double[] data_st_random = MathNet.Numerics.Generate.Standard(n);
            double[] result2 = new double[Lt.RowCount];
            MathNet.Numerics.Control.LinearAlgebraProvider.MatrixMultiply(result1, Lt.RowCount, Lt.ColumnCount, data_st_random, data_st_random.Length, 1, result2);
            // inv(Lt')
            double[] result3 = Lt.ToRowWiseArray();
            MathNet.Numerics.Control.LinearAlgebraProvider.LUInverse(result3, Lt.ColumnCount);
            // mldivide(Lt', (mdta{t}(:,1) - Xt * b) ./ tY2)
            double[] data_Xt = Xt.ToColumnWiseArray();
            double[] data_Xt_b = new double[Xt.RowCount];
            MathNet.Numerics.Control.LinearAlgebraProvider.MatrixMultiply(data_Xt, Xt.RowCount, Xt.ColumnCount, b, b.Length, 1, data_Xt_b);
            double[] data_1 = new double[Xt.RowCount];
            MathNet.Numerics.Control.LinearAlgebraProvider.SubtractArrays(data_c1, data_Xt_b, data_1);
            double alpha = 1.0 / tY2;
            double[] data_2 = new double[Xt.RowCount];
            MathNet.Numerics.Control.LinearAlgebraProvider.ScaleArray(alpha, data_1, data_2);
            double[] result4 = new double[Lt.ColumnCount];
            MathNet.Numerics.Control.LinearAlgebraProvider.MatrixMultiply(result3, Lt.ColumnCount, Lt.RowCount, data_2, data_2.Length, 1, result4);
            // mldivide(Lt, mldivide(Lt', (mdta{t}(:,1) - Xt * b) ./ tY2))
            double[] result5 = new double[Lt.RowCount];
            MathNet.Numerics.Control.LinearAlgebraProvider.MatrixMultiply(result1, Lt.RowCount, Lt.ColumnCount, result4, result4.Length, 1, result5);

            double[] result = new double[result5.Length];
            MathNet.Numerics.Control.LinearAlgebraProvider.AddArrays(result5, result2, result);
            return Matrix<double>.Build.DenseOfColumnMajor(result.Length, 1, result);
        }
        /// <summary>
        /// 计算矩阵 b0s；
        /// 计算公式：b0s= mldivide(Lt, mldivide(Lt', (mdta{t}(:,1) - Xt * b) ./ tY2)) + mldivide(Lt, randn(n,1))；
        /// </summary>
        [Obsolete("以前的版本")]
        public static Matrix<double> ComputeMatrix_b0s(Matrix<double> Lt, double[] data_c1, double[] data_Xt, int row_Xt, int column_Xt, double[] b, double tY2, int n)
        {
            // inv(Lt)
            double[] result1 = Lt.ToColumnWiseArray();
            MathNet.Numerics.Control.LinearAlgebraProvider.LUInverse(result1, Lt.RowCount);
            // mldivide(Lt, randn(n,1))
            double[] data_st_random = MathNet.Numerics.Generate.Standard(n);
            double[] result2 = new double[Lt.RowCount];
            MathNet.Numerics.Control.LinearAlgebraProvider.MatrixMultiply(result1, Lt.RowCount, Lt.ColumnCount, data_st_random, data_st_random.Length, 1, result2);
            // inv(Lt')
            double[] result3 = Lt.ToRowWiseArray();
            MathNet.Numerics.Control.LinearAlgebraProvider.LUInverse(result3, Lt.ColumnCount);
            // mldivide(Lt', (mdta{t}(:,1) - Xt * b) ./ tY2)
            //double[] data_Xt = Xt.ToColumnWiseArray();
            double[] data_Xt_b = new double[row_Xt];
            MathNet.Numerics.Control.LinearAlgebraProvider.MatrixMultiply(data_Xt, row_Xt, column_Xt, b, b.Length, 1, data_Xt_b);
            double[] data_1 = new double[row_Xt];
            MathNet.Numerics.Control.LinearAlgebraProvider.SubtractArrays(data_c1, data_Xt_b, data_1);
            double alpha = 1.0 / tY2;
            double[] data_2 = new double[row_Xt];
            MathNet.Numerics.Control.LinearAlgebraProvider.ScaleArray(alpha, data_1, data_2);
            double[] result4 = new double[Lt.ColumnCount];
            MathNet.Numerics.Control.LinearAlgebraProvider.MatrixMultiply(result3, Lt.ColumnCount, Lt.RowCount, data_2, data_2.Length, 1, result4);
            // mldivide(Lt, mldivide(Lt', (mdta{t}(:,1) - Xt * b) ./ tY2))
            double[] result5 = new double[Lt.RowCount];
            MathNet.Numerics.Control.LinearAlgebraProvider.MatrixMultiply(result1, Lt.RowCount, Lt.ColumnCount, result4, result4.Length, 1, result5);

            double[] result = new double[result5.Length];
            MathNet.Numerics.Control.LinearAlgebraProvider.AddArrays(result5, result2, result);
            return Matrix<double>.Build.DenseOfColumnMajor(result.Length, 1, result);
        }
        /// <summary>
        /// 计算数组 b0s；
        /// 计算公式：b0s= mldivide(Lt, mldivide(Lt', (mdta{t}(:,1) - Xt * b) ./ tY2)) + mldivide(Lt, randn(n,1))；
        /// </summary>
        public static double[] Compute_b0s(Matrix<double> Lt, double[] data_c1, double[] data_Xt, int row_Xt, int column_Xt, double[] b, double tY2, int n)
        {
            // inv(Lt)
            double[] data_1 = Lt.ToColumnWiseArray();
            MathNet.Numerics.Control.LinearAlgebraProvider.LUInverse(data_1, Lt.RowCount);

            // mldivide(Lt, randn(n,1))
            double[] data_st_random = MathNet.Numerics.Generate.Standard(n);
            double[] data_2 = new double[Lt.RowCount];
            MathNet.Numerics.Control.LinearAlgebraProvider.MatrixMultiply(data_1, Lt.RowCount, Lt.ColumnCount, data_st_random, data_st_random.Length, 1, data_2);

            // inv(Lt')
            double[] data_3 = Lt.ToRowWiseArray();
            MathNet.Numerics.Control.LinearAlgebraProvider.LUInverse(data_3, Lt.ColumnCount);

            // Xt * b
            double[] data_4 = new double[row_Xt];
            MathNet.Numerics.Control.LinearAlgebraProvider.MatrixMultiply(data_Xt, row_Xt, column_Xt, b, b.Length, 1, data_4);

            // mdta{t}(:,1) - Xt * b
            double[] data_5 = new double[row_Xt];
            MathNet.Numerics.Control.LinearAlgebraProvider.SubtractArrays(data_c1, data_4, data_5);

            // (mdta{t}(:,1) - Xt * b) ./ tY2
            double alpha = 1.0 / tY2;
            double[] data_6 = new double[row_Xt];
            MathNet.Numerics.Control.LinearAlgebraProvider.ScaleArray(alpha, data_5, data_6);

            // mldivide(Lt', (mdta{t}(:,1) - Xt * b) ./ tY2)
            double[] data_7 = new double[Lt.ColumnCount];
            MathNet.Numerics.Control.LinearAlgebraProvider.MatrixMultiply(data_3, Lt.ColumnCount, Lt.RowCount, data_6, data_6.Length, 1, data_7);

            // mldivide(Lt, mldivide(Lt', (mdta{t}(:,1) - Xt * b) ./ tY2))
            double[] data_8 = new double[Lt.RowCount];
            MathNet.Numerics.Control.LinearAlgebraProvider.MatrixMultiply(data_1, Lt.RowCount, Lt.ColumnCount, data_7, data_7.Length, 1, data_8);

            double[] result = new double[data_8.Length];
            MathNet.Numerics.Control.LinearAlgebraProvider.AddArrays(data_8, data_2, result);

            data_1 = new double[0]; data_2 = new double[0]; data_3 = new double[0]; data_4 = new double[0];
            data_5 = new double[0]; data_6 = new double[0]; data_7 = new double[0]; data_8 = new double[0];
            data_st_random = new double[0];
            return result;
        }
        /// <summary>
        /// 计算 s2b，tY2，按顺序返回；
        /// 计算公式：s2b = 1 / gamrnd(vctt1, 2 / (b0s' * b0sinvmat * b0s + vctt2))；
        /// tY2 = 1 / gamrnd(vctt3, 2 ./ (tY2_b + (mdta{t}(:,1) - b0s)' * (W - Xt * mldivide(Xt' * Xt, Xt')) * (mdta{t}(:,1) - b0s)))；
        /// </summary>
        [Obsolete("以前的版本")]
        public static Tuple<double, double> Compute_s2b_tY2(double vctt1, Matrix<double> b0s, Matrix<double> b0sinvmat, double vctt2, double vctt3, double tY2_b, double[] data_c1, double[] W, Matrix<double> Xt)
        {
            // s2b
            double[] data_b0s = b0s.ToColumnWiseArray();
            double[] data_b0s_trans = b0s.ToRowWiseArray();
            double[] data_b0sinvmat = b0sinvmat.ToColumnWiseArray();
            double[] data_1 = new double[b0s.ColumnCount * b0sinvmat.ColumnCount];
            MathNet.Numerics.Control.LinearAlgebraProvider.MatrixMultiply(data_b0s_trans, b0s.ColumnCount, b0s.RowCount, data_b0sinvmat, b0sinvmat.RowCount, b0sinvmat.ColumnCount, data_1);
            double[] data_2 = new double[1];
            MathNet.Numerics.Control.LinearAlgebraProvider.MatrixMultiply(data_1, b0s.ColumnCount, b0sinvmat.ColumnCount, data_b0s, b0s.RowCount, b0s.ColumnCount, data_2);
            double s2b = data_2[0] + vctt2;
            s2b = 2.0 / s2b;
            s2b = MathNet.Numerics.Distributions.Gamma.WithShapeScale(vctt1, s2b).Sample();
            s2b = 1.0 / s2b;
            // tY2
            // (mdta{t}(:,1) - b0s)
            data_1 = new double[data_c1.Length];
            MathNet.Numerics.Control.LinearAlgebraProvider.SubtractArrays(data_c1, data_b0s, data_1);
            // (W - Xt * mldivide(Xt' * Xt, Xt'))
            double[] data_Xt = Xt.ToColumnWiseArray();
            double[] data_Xt_trans = Xt.ToRowWiseArray();
            data_2 = new double[Xt.ColumnCount * Xt.ColumnCount];
            MathNet.Numerics.Control.LinearAlgebraProvider.MatrixMultiply(data_Xt_trans, Xt.ColumnCount, Xt.RowCount, data_Xt, Xt.RowCount, Xt.ColumnCount, data_2);
            MathNet.Numerics.Control.LinearAlgebraProvider.LUInverse(data_2, Xt.ColumnCount);
            double[] data_3 = new double[data_Xt.Length];
            MathNet.Numerics.Control.LinearAlgebraProvider.MatrixMultiply(data_2, Xt.ColumnCount, Xt.ColumnCount, data_Xt_trans, Xt.ColumnCount, Xt.RowCount, data_3);
            double[] data_4 = new double[Xt.RowCount * Xt.RowCount];
            MathNet.Numerics.Control.LinearAlgebraProvider.MatrixMultiply(data_Xt, Xt.RowCount, Xt.ColumnCount, data_3, Xt.ColumnCount, Xt.RowCount, data_4);
            double[] data_5 = new double[W.Length];
            MathNet.Numerics.Control.LinearAlgebraProvider.SubtractArrays(W, data_4, data_5);
            double[] result1 = new double[Xt.RowCount];
            MathNet.Numerics.Control.LinearAlgebraProvider.MatrixMultiply(data_1, 1, data_1.Length, data_5, Xt.RowCount, Xt.RowCount, result1);
            double[] result2 = new double[1];
            MathNet.Numerics.Control.LinearAlgebraProvider.MatrixMultiply(result1, 1, result1.Length, data_1, data_1.Length, 1, result2);
            double tY2 = tY2_b + result2[0];
            tY2 = 2.0 / tY2;
            tY2 = MathNet.Numerics.Distributions.Gamma.WithShapeScale(vctt3, tY2).Sample();
            tY2 = 1.0 / tY2;

            return Tuple.Create<double, double>(s2b, tY2);
        }
        /// <summary>
        /// 计算 s2b；
        /// 计算公式：s2b = 1 / gamrnd(vctt1, 2 / (b0s' * b0sinvmat * b0s + vctt2))；
        /// </summary>
        public static double Compute_s2b(double vctt1, double[] data_b0s, double[] data_b0sinvmat, int order_b0sinvmat, double vctt2)
        {
            // b0s' * b0sinvmat * b0s
            double[] data_1 = new double[order_b0sinvmat];
            MathNet.Numerics.Control.LinearAlgebraProvider.MatrixMultiply(data_b0s, 1, order_b0sinvmat, data_b0sinvmat, order_b0sinvmat, order_b0sinvmat, data_1);
            double[] data_2 = new double[1];
            MathNet.Numerics.Control.LinearAlgebraProvider.MatrixMultiply(data_1, 1, order_b0sinvmat, data_b0s, order_b0sinvmat, 1, data_2);

            double s2b = data_2[0] + vctt2;
            s2b = 2.0 / s2b;
            s2b = MathNet.Numerics.Distributions.Gamma.WithShapeScale(vctt1, s2b).Sample();
            s2b = 1.0 / s2b;

            data_1 = new double[0]; data_2 = new double[0];
            return s2b;
        }
        /// <summary>
        /// 计算 tY2；
        /// 计算公式：tY2 = 1 / gamrnd(vctt3, 2 ./ (tY2_b + (mdta{t}(:,1) - b0s)' * (W - Xt * mldivide(Xt' * Xt, Xt')) * (mdta{t}(:,1) - b0s)))；
        /// </summary>
        public static double Compute_tY2(double vctt3, double tY2_b, double[] data_c1_b0s, double[] W, double[] data_Xt, int row_Xt, int column_Xt, double[] data_XtXt, int order_XtXt, double[] data_Xt_trans)
        {
            // inv(Xt' * Xt)
            double[] data_1 = (double[])data_XtXt.Clone();
            MathNet.Numerics.Control.LinearAlgebraProvider.LUInverse(data_1, order_XtXt);

            // mldivide(Xt' * Xt, Xt')
            double[] data_2 = new double[data_Xt.Length];
            MathNet.Numerics.Control.LinearAlgebraProvider.MatrixMultiply(data_1, order_XtXt, order_XtXt, data_Xt_trans, column_Xt, row_Xt, data_2);

            // Xt * mldivide(Xt' * Xt, Xt')
            double[] data_3 = new double[row_Xt * row_Xt];
            MathNet.Numerics.Control.LinearAlgebraProvider.MatrixMultiply(data_Xt, row_Xt, column_Xt, data_2, column_Xt, row_Xt, data_3);

            // W - Xt * mldivide(Xt' * Xt, Xt')
            double[] data_4 = new double[W.Length];
            MathNet.Numerics.Control.LinearAlgebraProvider.SubtractArrays(W, data_3, data_4);

            // (mdta{t}(:,1) - b0s)' * (W - Xt * mldivide(Xt' * Xt, Xt'))
            double[] data_5 = new double[row_Xt];
            MathNet.Numerics.Control.LinearAlgebraProvider.MatrixMultiply(data_c1_b0s, 1, row_Xt, data_4, row_Xt, row_Xt, data_5);

            // (mdta{t}(:,1) - b0s)' * (W - Xt * mldivide(Xt' * Xt, Xt')) * (mdta{t}(:,1) - b0s)
            double[] data_6 = new double[1];
            MathNet.Numerics.Control.LinearAlgebraProvider.MatrixMultiply(data_5, 1, row_Xt, data_c1_b0s, row_Xt, 1, data_6);

            double tY2 = tY2_b + data_6[0];
            tY2 = 2.0 / tY2;
            tY2 = MathNet.Numerics.Distributions.Gamma.WithShapeScale(vctt3, tY2).Sample();
            tY2 = 1.0 / tY2;

            data_1 = new double[0]; data_2 = new double[0]; data_3 = new double[0]; data_4 = new double[0]; data_5 = new double[0];
            return tY2;
        }
        /// <summary>
        /// 计算矩阵 Lt；
        /// 计算公式：Lt = chol(Xt' * Xt ./ tY2 + vctt4)；
        /// </summary>
        [Obsolete("以前的版本")]
        public static Matrix<double> ComputeMatrix_Lt_2(Matrix<double> Xt, double tY2, double[] vctt4)
        {
            double[] data_Xt = Xt.ToColumnWiseArray();
            double[] data_Xt_trans = Xt.ToRowWiseArray();
            double[] result1 = new double[Xt.ColumnCount * Xt.ColumnCount];
            MathNet.Numerics.Control.LinearAlgebraProvider.MatrixMultiply(data_Xt_trans, Xt.ColumnCount, Xt.RowCount, data_Xt, Xt.RowCount, Xt.ColumnCount, result1);
            double[] result2 = new double[result1.Length];
            MathNet.Numerics.Control.LinearAlgebraProvider.ScaleArray(1.0 / tY2, result1, result2);
            double[] result = new double[vctt4.Length];
            MathNet.Numerics.Control.LinearAlgebraProvider.AddArrays(result2, vctt4, result);
            MathNet.Numerics.Control.LinearAlgebraProvider.CholeskyFactor(result, Xt.ColumnCount);
            var Lt = Matrix<double>.Build.DenseOfColumnMajor(Xt.ColumnCount, Xt.ColumnCount, result);
            return Lt.Transpose();
        }
        /// <summary>
        /// 计算矩阵 Lt；
        /// 计算公式：Lt = chol(Xt' * Xt ./ tY2 + vctt4)；
        /// </summary>
        public static Matrix<double> ComputeMatrix_Lt_2(double[] data_XtXt, int order_XtXt, double tY2, double[] vctt4)
        {
            // Xt' * Xt ./ tY2
            double alpha = 1.0 / tY2;
            double[] data_1 = new double[data_XtXt.Length];
            MathNet.Numerics.Control.LinearAlgebraProvider.ScaleArray(alpha, data_XtXt, data_1);

            // chol(Xt' * Xt ./ tY2 + vctt4)
            double[] data_2 = new double[data_1.Length];
            MathNet.Numerics.Control.LinearAlgebraProvider.AddArrays(data_1, vctt4, data_2);
            MathNet.Numerics.Control.LinearAlgebraProvider.CholeskyFactor(data_2, order_XtXt);

            var Lt = Matrix<double>.Build.DenseOfColumnMajor(order_XtXt, order_XtXt, data_2);
            data_1 = new double[0]; data_2 = new double[0];
            return Lt.Transpose();
        }
        /// <summary>
        /// 计算数组 b；
        /// 计算公式：b = mldivide(Lt, mldivide(Lt', Xt' * (mdta{t}(:,1) - b0s) ./ tY2)) + mldivide(Lt, randn(2,1))；
        /// </summary>
        [Obsolete("以前的版本")]
        public static double[] Compute_b(Matrix<double> Lt, Matrix<double> Xt, double[] data_c1, Matrix<double> b0s, double tY2)
        {
            double[] data_Lt = Lt.ToColumnWiseArray();
            double[] data_Lt_trans = Lt.ToRowWiseArray();
            double[] data_Xt_trans = Xt.ToRowWiseArray();
            double[] data_b0s = b0s.ToColumnWiseArray();
            // inv(Lt)
            double[] result_1 = (double[])data_Lt.Clone();
            MathNet.Numerics.Control.LinearAlgebraProvider.LUInverse(result_1, Lt.RowCount);
            // inv(Lt')
            double[] result_2 = (double[])data_Lt_trans.Clone();
            MathNet.Numerics.Control.LinearAlgebraProvider.LUInverse(result_2, Lt.RowCount);
            // mldivide(Lt', Xt' * (mdta{t}(:,1) - b0s) ./ tY2)
            double[] data_1 = new double[data_c1.Length];
            MathNet.Numerics.Control.LinearAlgebraProvider.SubtractArrays(data_c1, data_b0s, data_1);
            double[] data_2 = new double[Xt.ColumnCount];
            MathNet.Numerics.Control.LinearAlgebraProvider.MatrixMultiply(data_Xt_trans, Xt.ColumnCount, Xt.RowCount, data_1, data_1.Length, 1, data_2);
            double[] data_3 = new double[Xt.ColumnCount];
            MathNet.Numerics.Control.LinearAlgebraProvider.ScaleArray(1.0 / tY2, data_2, data_3);
            double[] result_3 = new double[Lt.ColumnCount];
            MathNet.Numerics.Control.LinearAlgebraProvider.MatrixMultiply(result_2, Lt.ColumnCount, Lt.RowCount, data_3, data_3.Length, 1, result_3);
            // mldivide(Lt, mldivide(Lt', Xt' * (mdta{t}(:,1) - b0s) ./ tY2))
            double[] result_4 = new double[Lt.RowCount];
            MathNet.Numerics.Control.LinearAlgebraProvider.MatrixMultiply(result_1, Lt.RowCount, Lt.ColumnCount, result_3, result_3.Length, 1, result_4);
            // mldivide(Lt, randn(2,1))
            data_1 = MathNet.Numerics.Generate.Standard(2);
            double[] result_5 = new double[Lt.RowCount];
            MathNet.Numerics.Control.LinearAlgebraProvider.MatrixMultiply(result_1, Lt.RowCount, Lt.ColumnCount, data_1, data_1.Length, 1, result_5);

            double[] result = new double[result_5.Length];
            MathNet.Numerics.Control.LinearAlgebraProvider.AddArrays(result_4, result_5, result);
            return result;
        }
        /// <summary>
        /// 计算数组 b；
        /// 计算公式：b = mldivide(Lt, mldivide(Lt', Xt' * (mdta{t}(:,1) - b0s) ./ tY2)) + mldivide(Lt, randn(2,1))；
        /// </summary>
        public static double[] Compute_b(Matrix<double> Lt, double[] data_Xt_trans, int row_Xt, int column_Xt, double[] data_c1_b0s, double tY2)
        {
            // inv(Lt)
            double[] data_1 = Lt.ToColumnWiseArray();
            MathNet.Numerics.Control.LinearAlgebraProvider.LUInverse(data_1, Lt.RowCount);

            // inv(Lt')
            double[] data_2 = Lt.ToRowWiseArray();
            MathNet.Numerics.Control.LinearAlgebraProvider.LUInverse(data_2, Lt.ColumnCount);

            // Xt' * (mdta{t}(:,1) - b0s)
            double[] data_3 = new double[column_Xt];
            MathNet.Numerics.Control.LinearAlgebraProvider.MatrixMultiply(data_Xt_trans, column_Xt, row_Xt, data_c1_b0s, data_c1_b0s.Length, 1, data_3);

            // Xt' * (mdta{t}(:,1) - b0s) ./ tY2
            double alpha = 1.0 / tY2;
            double[] data_4 = new double[column_Xt];
            MathNet.Numerics.Control.LinearAlgebraProvider.ScaleArray(alpha, data_3, data_4);

            // mldivide(Lt', Xt' * (mdta{t}(:,1) - b0s) ./ tY2)
            double[] data_5 = new double[Lt.ColumnCount];
            MathNet.Numerics.Control.LinearAlgebraProvider.MatrixMultiply(data_2, Lt.ColumnCount, Lt.RowCount, data_4, column_Xt, 1, data_5);

            // mldivide(Lt, mldivide(Lt', Xt' * (mdta{t}(:,1) - b0s) ./ tY2))
            double[] data_6 = new double[Lt.RowCount];
            MathNet.Numerics.Control.LinearAlgebraProvider.MatrixMultiply(data_1, Lt.RowCount, Lt.ColumnCount, data_5, Lt.ColumnCount, 1, data_6);

            // mldivide(Lt, randn(2,1))
            double[] data_7 = MathNet.Numerics.Generate.Standard(2);
            double[] data_8 = new double[Lt.RowCount];
            MathNet.Numerics.Control.LinearAlgebraProvider.MatrixMultiply(data_1, Lt.RowCount, Lt.ColumnCount, data_7, 2, 1, data_8);

            double[] result = new double[data_8.Length];
            MathNet.Numerics.Control.LinearAlgebraProvider.AddArrays(data_6, data_8, result);

            data_1 = new double[0]; data_2 = new double[0]; data_3 = new double[0]; data_4 = new double[0];
            data_5 = new double[0]; data_6 = new double[0]; data_7 = new double[0]; data_8 = new double[0];
            return result;
        }
        /// <summary>
        /// 计算数组 prop；
        /// 计算公式：prop = s{4} * normrnd(0, 1, length(Q), 1)；
        /// prop = prop - mean(prop)；
        /// </summary>
        [Obsolete("以前的版本")]
        public static double[] Compute_prop(Matrix<double> s4, double[] Q)
        {
            double[] data_random = MathNet.Numerics.Generate.Normal(Q.Length, 0.0, 1.0);
            double[] data_s4 = s4.ToColumnWiseArray();
            double[] prop = new double[s4.RowCount];
            MathNet.Numerics.Control.LinearAlgebraProvider.MatrixMultiply(data_s4, s4.RowCount, s4.ColumnCount, data_random, data_random.Length, 1, prop);
            //double mean = prop.Average();
            //prop = prop.Select((double a) => a - mean).ToArray();
            double mean = 0.0;
            int len = prop.Length;
            for (int i = 0; i < len; i++)
                mean = mean + prop[i];
            mean = mean / len;
            for (int i = 0; i < len; i++)
                prop[i] = prop[i] - mean;
            return prop;
        }
        /// <summary>
        /// 计算数组 prop；
        /// 计算公式：prop = s{4} * normrnd(0, 1, length(Q), 1)；
        /// prop = prop - mean(prop)；
        /// </summary>
        public static double[] Compute_prop(double[] data_s4, int row_s4, int column_s4, double[] Q)
        {
            // normrnd(0, 1, length(Q), 1)
            double[] data_1 = MathNet.Numerics.Generate.Normal(Q.Length, 0.0, 1.0);

            // s{4} * normrnd(0, 1, length(Q), 1)
            double[] prop = new double[row_s4];
            MathNet.Numerics.Control.LinearAlgebraProvider.MatrixMultiply(data_s4, row_s4, column_s4, data_1, data_1.Length, 1, prop);
            data_1 = new double[0];

            // mean(prop)
            double mean = prop.Average();

            // prop - mean(prop)
            prop = prop.Select((double a) => a - mean).ToArray();
            return prop;
        }
        /// <summary>
        /// 计算参数 alph；
        /// 计算公式：alph = -((mdta{t}(:,1) - Xtp * b - b0s)' * (mdta{t}(:,1) - Xtp * b - b0s)) + ((mdta{t}(:,1) - Xt * b - b0s)' * (mdta{t}(:,1) - Xt * b - b0s))；
        /// </summary>
        [Obsolete("以前的版本")]
        public static double Compute_alph(double[] data_c1, Matrix<double> Xtp, double[] b, Matrix<double> b0s, Matrix<double> Xt)
        {
            double[] data_Xtp = Xtp.ToColumnWiseArray();
            double[] data_Xt = Xt.ToColumnWiseArray();
            double[] data_b0s = b0s.ToColumnWiseArray();
            // mdta{t}(:,1) - b0s
            double[] result1 = new double[data_c1.Length];
            MathNet.Numerics.Control.LinearAlgebraProvider.SubtractArrays(data_c1, data_b0s, result1);
            // Xtp * b
            double[] result2 = new double[Xtp.RowCount];
            MathNet.Numerics.Control.LinearAlgebraProvider.MatrixMultiply(data_Xtp, Xtp.RowCount, Xtp.ColumnCount, b, b.Length, 1, result2);
            // Xt * b
            double[] result3 = new double[Xt.RowCount];
            MathNet.Numerics.Control.LinearAlgebraProvider.MatrixMultiply(data_Xt, Xt.RowCount, Xt.ColumnCount, b, b.Length, 1, result3);
            // (mdta{t}(:,1) - Xtp * b - b0s)
            double[] result4 = new double[data_c1.Length];
            MathNet.Numerics.Control.LinearAlgebraProvider.SubtractArrays(result1, result2, result4);
            // (mdta{t}(:,1) - Xt * b - b0s)
            double[] result5 = new double[data_c1.Length];
            MathNet.Numerics.Control.LinearAlgebraProvider.SubtractArrays(result1, result3, result5);

            result1 = new double[1];
            MathNet.Numerics.Control.LinearAlgebraProvider.MatrixMultiply(result4, 1, Xtp.RowCount, result4, Xtp.RowCount, 1, result1);
            double alph = -result1[0];
            result2 = new double[1];
            MathNet.Numerics.Control.LinearAlgebraProvider.MatrixMultiply(result5, 1, Xt.RowCount, result5, Xt.RowCount, 1, result2);
            alph += result2[0];
            return alph;
        }
        /// <summary>
        /// 计算参数 alph；
        /// 计算公式：alph = -((mdta{t}(:,1) - Xtp * b - b0s)' * (mdta{t}(:,1) - Xtp * b - b0s)) + ((mdta{t}(:,1) - Xt * b - b0s)' * (mdta{t}(:,1) - Xt * b - b0s))；
        /// </summary>
        public static double Compute_alph(double[] data_c1_b0s, Matrix<double> Xtp, double[] b, double[] data_Xt, int row_Xt, int column_Xt)
        {
            double[] data_Xtp = Xtp.ToColumnWiseArray();
            // Xtp * b
            double[] data_1 = new double[Xtp.RowCount];
            MathNet.Numerics.Control.LinearAlgebraProvider.MatrixMultiply(data_Xtp, Xtp.RowCount, Xtp.ColumnCount, b, b.Length, 1, data_1);

            // Xt * b
            double[] data_2 = new double[row_Xt];
            MathNet.Numerics.Control.LinearAlgebraProvider.MatrixMultiply(data_Xt, row_Xt, column_Xt, b, b.Length, 1, data_2);

            // (mdta{t}(:,1) - Xtp * b - b0s)
            double[] data_3 = new double[data_c1_b0s.Length];
            MathNet.Numerics.Control.LinearAlgebraProvider.SubtractArrays(data_c1_b0s, data_1, data_3);

            // (mdta{t}(:,1) - Xt * b - b0s)
            double[] data_4 = new double[data_c1_b0s.Length];
            MathNet.Numerics.Control.LinearAlgebraProvider.SubtractArrays(data_c1_b0s, data_2, data_4);

            data_1 = new double[1];
            MathNet.Numerics.Control.LinearAlgebraProvider.MatrixMultiply(data_3, 1, Xtp.RowCount, data_3, Xtp.RowCount, 1, data_1);
            double alph = -data_1[0];
            data_2 = new double[1];
            MathNet.Numerics.Control.LinearAlgebraProvider.MatrixMultiply(data_4, 1, row_Xt, data_4, row_Xt, 1, data_2);
            alph += data_2[0];

            data_Xtp = new double[0]; data_3 = new double[0]; data_4 = new double[0];
            return alph;
        }
        /// <summary>
        /// 判断 alph；
        /// 计算公式：all(eig(tY2 .* mldivide(Xtp' * Xtp, eye(2))) > 0)；
        /// </summary>
        public static bool Judge_tY2_Xtp(double tY2, Matrix<double> Xtp)
        {
            double[] data_Xtp = Xtp.ToColumnWiseArray();
            double[] data_Xtp_trans = Xtp.ToRowWiseArray();
            // inv(Xtp' * Xtp)
            double[] data_1 = new double[4];
            MathNet.Numerics.Control.LinearAlgebraProvider.MatrixMultiply(data_Xtp_trans, 2, Xtp.RowCount, data_Xtp, Xtp.RowCount, 2, data_1);
            MathNet.Numerics.Control.LinearAlgebraProvider.LUInverse(data_1, 2);

            // mldivide(Xtp' * Xtp, eye(2))
            //double[] data_2 = new double[Xtp.ColumnCount * 2];
            //MathNet.Numerics.Control.LinearAlgebraProvider.MatrixMultiply(data_1, Xtp.ColumnCount, Xtp.ColumnCount, new double[] { 1.0, 0.0, 0.0, 1.0 }, 2, 2, data_2);

            double[] data_3 = new double[4];
            MathNet.Numerics.Control.LinearAlgebraProvider.ScaleArray(tY2, data_1, data_3);
            double[] matrixEv = new double[4];
            System.Numerics.Complex[] vectorEv = new System.Numerics.Complex[2];
            double[] matrixD = new double[4];
            MathNet.Numerics.Control.LinearAlgebraProvider.EigenDecomp(true, 2, data_3, matrixEv, vectorEv, matrixD);

            data_Xtp = new double[0]; data_Xtp_trans = new double[0];
            return (vectorEv[0].Real > 0.0) && (vectorEv[1].Real > 0.0);
        }
        /// <summary>
        /// 计算数组 pmean pvar，按顺序返回；
        /// 计算公式：pmean = b(1) + b(2) .* cmaq{2}(:,1) + krigmat * b0s；
        /// pvar = s2b .* (1 - mvars) + tY2；
        /// </summary>
        [Obsolete("以前的版本")]
        public static Tuple<double[], double[]> Compute_pmean_pvar(double[] b, Matrix<double> matrixModelConc, Matrix<double> krigmat, Matrix<double> b0s, double s2b, Matrix<double> mvars, double tY2)
        {
            double[] data_cmaq2 = matrixModelConc.Column(0).ToArray();
            double[] data_krigmat = krigmat.ToColumnWiseArray();
            double[] data_b0s = b0s.ToColumnWiseArray();
            double[] data_mvars = mvars.ToColumnWiseArray();
            // pmean
            double[] data_1 = new double[data_cmaq2.Length];
            MathNet.Numerics.Control.LinearAlgebraProvider.ScaleArray(b[1], data_cmaq2, data_1);
            double[] data_2 = new double[krigmat.RowCount * b0s.ColumnCount];
            MathNet.Numerics.Control.LinearAlgebraProvider.MatrixMultiply(data_krigmat, krigmat.RowCount, krigmat.ColumnCount, data_b0s, b0s.RowCount, b0s.ColumnCount, data_2);
            double[] pmean = new double[data_1.Length];
            MathNet.Numerics.Control.LinearAlgebraProvider.AddArrays(data_1, data_2, pmean);
            //pmean = pmean.Select((double a) => b[0] + a).ToArray();
            int len = pmean.Length;
            for (int i = 0; i < len; i++)
                pmean[i] = pmean[i] + b[0];
            // pvar
            //data_1 = data_mvars.Select((double a) => 1.0 - a).ToArray();
            len = data_mvars.Length;
            data_1 = new double[len];
            for (int i = 0; i < len; i++)
                data_1[i] = 1.0 - data_mvars[i];
            double[] pvar = new double[len];
            MathNet.Numerics.Control.LinearAlgebraProvider.ScaleArray(s2b, data_1, pvar);
            //pvar = pvar.Select((double a) => Math.Sqrt(a + tY2)).ToArray();
            for (int i = 0; i < len; i++)
                pvar[i] = Math.Sqrt(pvar[i] + tY2);

            return Tuple.Create<double[], double[]>(pmean, pvar);
        }
        /// <summary>
        /// 计算数组 pmean；
        /// 计算公式：pmean = b(1) + b(2) .* cmaq{2}(:,1) + krigmat * b0s；
        /// </summary>
        public static double[] Compute_pmean(double[] b, double[] data_matrixModelConc_c1, double[] data_krigmat, int row_krigmat, int column_krigmat, double[] data_b0s)
        {
            // cmaq{2}(:,1)
            //double[] data_1 = matrixModelConc.Column(0).ToArray();
            //double[] data_1 = (double[])data_matrixModelConc_c1.Clone();

            // b(2) .* cmaq{2}(:,1)
            double[] data_1 = new double[data_matrixModelConc_c1.Length];
            MathNet.Numerics.Control.LinearAlgebraProvider.ScaleArray(b[1], data_matrixModelConc_c1, data_1);

            // krigmat * b0s
            double[] data_2 = new double[row_krigmat];
            MathNet.Numerics.Control.LinearAlgebraProvider.MatrixMultiply(data_krigmat, row_krigmat, column_krigmat, data_b0s, data_b0s.Length, 1, data_2);

            double[] result = new double[row_krigmat];
            MathNet.Numerics.Control.LinearAlgebraProvider.AddArrays(data_1, data_2, result);

            data_1 = new double[0]; data_2 = new double[0];
            result = result.Select((double a) => b[0] + a).ToArray();
            return result;
        }
        /// <summary>
        /// 计算数组 pvar；
        /// 计算公式：pvar = s2b .* (1 - mvars) + tY2；
        /// </summary>
        public static double[] Compute_pvar(double s2b, double[] data_mvars, double tY2)
        {
            int len = data_mvars.Length;
            double[] data_1 = data_mvars.Select((double a) => 1.0 - a).ToArray();

            double[] result = new double[len];
            MathNet.Numerics.Control.LinearAlgebraProvider.ScaleArray(s2b, data_1, result);
            data_1 = new double[0];

            for (int i = 0; i < len; i++)
                result[i] = Math.Sqrt(result[i] + tY2);

            return result;
        }
        /// <summary>
        /// 计算数组 kY 和 kY 的平方，前面一步已经对 pvar 进行了 sqrt 操作；
        /// 计算公式：kY = max(0, normrnd(pmean, sqrt(pvar)))；
        /// </summary>
        public static Tuple<double[], double[]> Compute_kY(double[] pmean, double[] pvar)
        {
            int len = pmean.Length;
            double[] kY = new double[len];
            double[] kY2 = new double[len];
            for (int i = 0; i < len; i++)
            {
                kY[i] = Math.Max(0.0, MathNet.Numerics.Distributions.Normal.WithMeanStdDev(pmean[i], pvar[i]).Sample());
                kY2[i] = Math.Pow(kY[i], 2.0);
            }
            return Tuple.Create<double[], double[]>(kY, kY2);
        }
        /// <summary>
        /// 计算数组 Ybar Y2bar，按顺序返回；
        /// 计算公式：Ybar = (Ybar .* (kp - 1) + kY) / kp；
        /// Y2bar = (Y2bar .* (kp - 1) + (kY .^ 2)) / kp；
        /// </summary>
        [Obsolete("以前的版本")]
        public static Tuple<double[], double[]> Compute_Ybar_Y2bar(double[] Ybar, int kp, double[] kY, double[] Y2bar, double[] kY2)
        {
            // Ybar
            double[] data_1 = new double[Ybar.Length];
            MathNet.Numerics.Control.LinearAlgebraProvider.ScaleArray(kp - 1, Ybar, data_1);
            double[] data_2 = new double[data_1.Length];
            MathNet.Numerics.Control.LinearAlgebraProvider.AddArrays(data_1, kY, data_2);
            double[] result1 = new double[data_2.Length];
            MathNet.Numerics.Control.LinearAlgebraProvider.ScaleArray(1.0 / kp, data_2, result1);
            // Y2bar
            data_1 = new double[Y2bar.Length];
            MathNet.Numerics.Control.LinearAlgebraProvider.ScaleArray(kp - 1, Y2bar, data_1);
            data_2 = new double[data_1.Length];
            MathNet.Numerics.Control.LinearAlgebraProvider.AddArrays(data_1, kY2, data_2);
            double[] result2 = new double[data_2.Length];
            MathNet.Numerics.Control.LinearAlgebraProvider.ScaleArray(1.0 / kp, data_2, result2);

            return Tuple.Create<double[], double[]>(result1, result2);
        }
        /// <summary>
        /// 计算数组 Ybar；
        /// 计算公式：Ybar = (Ybar .* (kp - 1) + kY) / kp；
        /// </summary>
        public static double[] Compute_Ybar(double[] Ybar, int kp, double[] kY)
        {
            double alpha = kp - 1;
            double[] data_1 = new double[Ybar.Length];
            MathNet.Numerics.Control.LinearAlgebraProvider.ScaleArray(alpha, Ybar, data_1);

            double[] data_2 = new double[data_1.Length];
            MathNet.Numerics.Control.LinearAlgebraProvider.AddArrays(data_1, kY, data_2);

            alpha = 1.0 / kp;
            double[] result = new double[data_2.Length];
            MathNet.Numerics.Control.LinearAlgebraProvider.ScaleArray(alpha, data_2, result);

            data_1 = new double[0]; data_2 = new double[0];
            return result;
        }
        /// <summary>
        /// 计算数组 Y2bar；
        /// 计算公式：Y2bar = (Y2bar .* (kp - 1) + (kY .^ 2)) / kp；
        /// </summary>
        public static double[] Compute_Y2bar(double[] Y2bar, int kp, double[] kY2)
        {
            double alpha = kp - 1;
            double[] data_1 = new double[Y2bar.Length];
            MathNet.Numerics.Control.LinearAlgebraProvider.ScaleArray(alpha, Y2bar, data_1);

            double[] data_2 = new double[data_1.Length];
            MathNet.Numerics.Control.LinearAlgebraProvider.AddArrays(data_1, kY2, data_2);

            alpha = 1.0 / kp;
            double[] result = new double[data_2.Length];
            MathNet.Numerics.Control.LinearAlgebraProvider.ScaleArray(alpha, data_2, result);

            data_1 = new double[0]; data_2 = new double[0];
            return result;
        }
        /// <summary>
        /// 计算结果输出数据；
        /// 计算公式：sqrt((kp .* Y2bar - kp .* (Ybar .^ 2)) / (kp - 1))；
        /// </summary>
        public static double[] Compute_Result(int kp, double[] Y2bar, double[] Ybar)
        {
            double[] data_1 = (double[])Ybar.Clone();
            double[] data_2 = (double[])Y2bar.Clone();
            //data_1 = data_1.Select((double a) => Math.Pow(a, 2.0)).ToArray();
            int len = data_1.Length;
            for (int i = 0; i < len; i++)
                data_1[i] = Math.Pow(data_1[i], 2.0);
            double[] data_3 = new double[data_1.Length];
            MathNet.Numerics.Control.LinearAlgebraProvider.SubtractArrays(data_2, data_1, data_3);

            double alpha = (double)kp / (kp - 1);
            double[] data_4 = new double[data_3.Length];
            MathNet.Numerics.Control.LinearAlgebraProvider.ScaleArray(alpha, data_3, data_4);
            //data_4 = data_4.Select((double a) => Math.Sqrt(a)).ToArray();
            len = data_4.Length;
            for (int i = 0; i < len; i++)
                data_4[i] = Math.Sqrt(data_4[i]);

            data_1 = new double[0]; data_2 = new double[0]; data_3 = new double[0];
            return data_4;
        }
    }
}
