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

namespace SMAT_CE.DownscaleR.AnnualPM
{
    /// <summary>
    /// Annual PM downscaler analysis
    /// </summary>
    public class AnnualPMAnalysis : IDisposable
    {
        #region 字段
        // 外界传进来参与计算的参数
        private Dictionary<string, PmModelDataOutput> _dicModelData;
        private Dictionary<string, PMMonitorOutput> _dicSpeciesMonitorData;
        private Dictionary<string, PMMonitorOutput> _dicPmMonitorData;
        private CommonSettings _settings;

        #region 中间参数

        #endregion

        // 输出参数
        private StringBuilder _errormsg;
        private Dictionary<string, AnnualPMDSOutput> _dicDSOutputValue;
        private string[] _modelMissingSpecies;

        // 其他参数
        private bool _disposed = false;// Disposable状态
        #endregion

        #region 构造函数
        /// <summary>
        /// 初始化
        /// </summary>
        /// <param name="settings">设置参数</param>
        private AnnualPMAnalysis(CommonSettings settings)
        {
            _settings = (CommonSettings)settings.Clone();
            _errormsg = new StringBuilder();
            _dicDSOutputValue = new Dictionary<string, AnnualPMDSOutput>();
            _modelMissingSpecies = new string[] { };
        }
        /// <summary>
        /// 初始化
        /// </summary>
        /// <param name="settings">设置参数</param>
        /// <param name="dicModelData">模型数据，包括PM和Species</param>
        /// <param name="dicSpeciesMonitorData">Species监测数据</param>
        /// <param name="dicPmMonitorData">PM监测数据</param>
        public AnnualPMAnalysis(
            CommonSettings settings,
            Dictionary<string, PmModelDataOutput> dicModelData,
            Dictionary<string, PMMonitorOutput> dicSpeciesMonitorData,
            Dictionary<string, PMMonitorOutput> dicPmMonitorData)
            : this(settings)
        {
            _dicModelData = dicModelData.Select(k => k.Value).ToDictionary(l => l.id);
            _dicSpeciesMonitorData = dicSpeciesMonitorData.Select(k => k.Value).ToDictionary(l => l.id);
            _dicPmMonitorData = dicPmMonitorData.Select(k => k.Value).ToDictionary(l => l.id);
        }
        #endregion

        #region 属性
        /// <summary>
        /// Downscaler设置参数
        /// </summary>
        public CommonSettings Settings
        {
            get { return _settings; }
        }
        /// <summary>
        /// 错误消息
        /// </summary>
        public string ErrorMsg
        {
            get
            {
                if (_errormsg != null)
                    return _errormsg.ToString();
                return null;
            }
        }
        /// <summary>
        /// 输出结果，季度DS预测值
        /// </summary>
        public Dictionary<string, AnnualPMDSOutput> DSOutputValue
        {
            get { return _dicDSOutputValue; }
        }
        /// <summary>
        /// 模型文件没有的成分
        /// </summary>
        public string[] ModelMissingSpecies
        {
            get { return _modelMissingSpecies; }
        }
        #endregion

        #region 公共方法
        /// <summary>
        /// 开始分析
        /// </summary>
        public bool StartAnalysis()
        {
            try
            {
                if ((_dicModelData == null || _dicModelData.Count <= 0) ||
                    (_dicSpeciesMonitorData == null || _dicSpeciesMonitorData.Count <= 0) ||
                    (_dicPmMonitorData == null || _dicPmMonitorData.Count <= 0))
                {
                    _errormsg.AppendLine("smat-ce input datas are null or empty.");
                    return false;
                }
                if (_settings == null)
                {
                    _errormsg.AppendLine("DS setting is null.");
                    return false;
                }

                #region 获取模型文件没有的成分，主要是salt
                if (!string.IsNullOrEmpty(_settings.ModelFilePath))
                {
                    string species = "so4,no3r,oc,ec,crustal,salt,no3,nh4,don";
                    //string try_miss = "salt,don";
                    using (FileStream fs = new FileStream(_settings.ModelFilePath, FileMode.Open, FileAccess.Read, FileShare.ReadWrite))
                    {
                        using (StreamReader csv = new StreamReader(fs, System.Text.Encoding.UTF8))
                        {
                            string strLine = csv.ReadLine();// day or quarter
                            strLine = csv.ReadLine();
                            string[] strLineArray = strLine.Split(new char[] { ',' }, StringSplitOptions.RemoveEmptyEntries);
                            for (int i = 0; i < strLineArray.Length; i++)
                            {
                                string col = strLineArray[i].Trim().ToLower().Replace("\"", "");
                                if (col == "so4" || col == "no3r" || col == "oc" || col == "ec" || col == "crustal" || col == "salt" || col == "no3" || col == "nh4" || col == "don")
                                {
                                    if (species.Contains(col))
                                    {
                                        if (col == "no3r" || col == "no3") { species = species.Replace("no3r", ""); species = species.Replace("no3", ""); }
                                        else { species = species.Replace(col, ""); }
                                    }
                                }
                            }
                            _modelMissingSpecies = species.Split(new char[] { ',' }, StringSplitOptions.RemoveEmptyEntries);
                        }
                    }
                }
                #endregion

                #region 中间变量
                Matrix<double> matrixLatLonPM_Model, matrixLatLonPM_Monitor, matrixPM_Model, matrixPM_Monitor,
                matrixLatLonSO4_Model, matrixLatLonSO4_Monitor, matrixSO4_Model, matrixSO4_Monitor,
                matrixLatLonNO3R_Model, matrixLatLonNO3R_Monitor, matrixNO3R_Model, matrixNO3R_Monitor,
                matrixLatLonOCB_Model, matrixLatLonOCB_Monitor, matrixOCB_Model, matrixOCB_Monitor,
                matrixLatLonEC_Model, matrixLatLonEC_Monitor, matrixEC_Model, matrixEC_Monitor,
                matrixLatLonCrustal_Model, matrixLatLonCrustal_Monitor, matrixCrustal_Model, matrixCrustal_Monitor,
                matrixLatLonSalt_Model, matrixLatLonSalt_Monitor, matrixSalt_Model, matrixSalt_Monitor,
                matrixLatLonNO3_Model, matrixLatLonNO3_Monitor, matrixNO3_Model, matrixNO3_Monitor,
                matrixLatLonNH4_Model, matrixLatLonNH4_Monitor, matrixNH4_Model, matrixNH4_Monitor;
                #endregion

                #region 多线程，一个case要保证4个线程，因为计算过程本身就用了多线程
                if (_settings.MultiThreading)
                {

                }
                #endregion
                #region 单线程
                else
                {
                    string[] quarters = new string[4] { "1", "2", "3", "4" };
                    Tuple<double[], double[]> result;
                    double[] arrPrediction, arrSepred;
                    string q;
                    bool ok = false;
                    for (int i = 0; i < 4; i++)
                    {
                        matrixLatLonPM_Model = null; matrixLatLonPM_Monitor = null; matrixPM_Model = null; matrixPM_Monitor = null;
                        matrixLatLonSO4_Model = null; matrixLatLonSO4_Monitor = null; matrixSO4_Model = null; matrixSO4_Monitor = null;
                        matrixLatLonNO3R_Model = null; matrixLatLonNO3R_Monitor = null; matrixNO3R_Model = null; matrixNO3R_Monitor = null;
                        matrixLatLonOCB_Model = null; matrixLatLonOCB_Monitor = null; matrixOCB_Model = null; matrixOCB_Monitor = null;
                        matrixLatLonEC_Model = null; matrixLatLonEC_Monitor = null; matrixEC_Model = null; matrixEC_Monitor = null;
                        matrixLatLonCrustal_Model = null; matrixLatLonCrustal_Monitor = null; matrixCrustal_Model = null; matrixCrustal_Monitor = null;
                        matrixLatLonSalt_Model = null; matrixLatLonSalt_Monitor = null; matrixSalt_Model = null; matrixSalt_Monitor = null;
                        matrixLatLonNO3_Model = null; matrixLatLonNO3_Monitor = null; matrixNO3_Model = null; matrixNO3_Monitor = null;
                        matrixLatLonNH4_Model = null; matrixLatLonNH4_Monitor = null; matrixNH4_Model = null; matrixNH4_Monitor = null;
                        #region DS计算，没有的成分用VNA代替
                        q = quarters[i];
                        ok = MatrixingInputData(q, ref _dicDSOutputValue,
                                    ref matrixLatLonPM_Model, ref matrixLatLonPM_Monitor, ref matrixPM_Model, ref matrixPM_Monitor,
                                    ref matrixLatLonSO4_Model, ref matrixLatLonSO4_Monitor, ref matrixSO4_Model, ref matrixSO4_Monitor,
                                    ref matrixLatLonNO3R_Model, ref matrixLatLonNO3R_Monitor, ref matrixNO3R_Model, ref matrixNO3R_Monitor,
                                    ref matrixLatLonOCB_Model, ref matrixLatLonOCB_Monitor, ref matrixOCB_Model, ref matrixOCB_Monitor,
                                    ref matrixLatLonEC_Model, ref matrixLatLonEC_Monitor, ref matrixEC_Model, ref matrixEC_Monitor,
                                    ref matrixLatLonCrustal_Model, ref matrixLatLonCrustal_Monitor, ref matrixCrustal_Model, ref matrixCrustal_Monitor,
                                    ref matrixLatLonSalt_Model, ref matrixLatLonSalt_Monitor, ref matrixSalt_Model, ref matrixSalt_Monitor,
                                    ref matrixLatLonNO3_Model, ref matrixLatLonNO3_Monitor, ref matrixNO3_Model, ref matrixNO3_Monitor,
                                    ref matrixLatLonNH4_Model, ref matrixLatLonNH4_Monitor, ref matrixNH4_Model, ref matrixNH4_Monitor);
                        if (ok)
                        {
                            #region 计算PM25
                            if (matrixLatLonPM_Model != null && matrixLatLonPM_Monitor != null && matrixPM_Model != null && matrixPM_Monitor != null)
                            {
                                //SaveMatlabInputData(matrixLatLonPM_Model, matrixLatLonPM_Monitor, matrixPM_Model, matrixPM_Monitor, "pm25", _settings.Cmaqres.ToString(), q);
                                result = SingleRun(matrixLatLonPM_Model, matrixLatLonPM_Monitor, matrixPM_Model, matrixPM_Monitor);
                                if (result != null)
                                {
                                    arrPrediction = result.Item1;
                                    arrSepred = result.Item2;
                                    //SaveDSOutputData(matrixLatLonPM_Model, arrPrediction, arrSepred, "pm25", _settings.Cmaqres.ToString(), q);
                                    foreach (KeyValuePair<string, AnnualPMDSOutput> kv in _dicDSOutputValue)
                                    {
                                        for (int j = 0; j < matrixLatLonPM_Model.RowCount; j++)
                                        {
                                            if (kv.Value.Lat == matrixLatLonPM_Model.At(j, 0) && kv.Value.Lon == matrixLatLonPM_Model.At(j, 1))
                                            {
                                                if (kv.Value.DicQuarterlyDSValue.ContainsKey(q))
                                                {
                                                    kv.Value.DicQuarterlyDSValue[q].PM25_DS_Prediction = arrPrediction[j];
                                                    kv.Value.DicQuarterlyDSValue[q].PM25_DS_Sepred = arrSepred[j];
                                                }
                                                break;
                                            }
                                        }
                                    }
                                }
                                GC.Collect();
                            }
                            #endregion
                            #region 计算SO4
                            if (!_modelMissingSpecies.Contains("so4"))
                            {
                                if (matrixLatLonSO4_Model != null && matrixLatLonSO4_Monitor != null && matrixSO4_Model != null && matrixSO4_Monitor != null)
                                {
                                    //SaveMatlabInputData(matrixLatLonSO4_Model, matrixLatLonSO4_Monitor, matrixSO4_Model, matrixSO4_Monitor, "so4", _settings.Cmaqres.ToString(), q);
                                    result = SingleRun(matrixLatLonSO4_Model, matrixLatLonSO4_Monitor, matrixSO4_Model, matrixSO4_Monitor);
                                    if (result != null)
                                    {
                                        arrPrediction = result.Item1;
                                        arrSepred = result.Item2;
                                        //SaveDSOutputData(matrixLatLonSO4_Model, arrPrediction, arrSepred, "so4", _settings.Cmaqres.ToString(), q);
                                        foreach (KeyValuePair<string, AnnualPMDSOutput> kv in _dicDSOutputValue)
                                        {
                                            for (int j = 0; j < matrixLatLonSO4_Model.RowCount; j++)
                                            {
                                                if (kv.Value.Lat == matrixLatLonSO4_Model.At(j, 0) && kv.Value.Lon == matrixLatLonSO4_Model.At(j, 1))
                                                {
                                                    if (kv.Value.DicQuarterlyDSValue.ContainsKey(q))
                                                    {
                                                        kv.Value.DicQuarterlyDSValue[q].SO4_DS_Prediction = arrPrediction[j];
                                                        kv.Value.DicQuarterlyDSValue[q].SO4_DS_Sepred = arrSepred[j];
                                                    }
                                                    break;
                                                }
                                            }
                                        }
                                    }
                                    GC.Collect();
                                }
                            }
                            #endregion
                            #region 计算NO3R
                            if (!_modelMissingSpecies.Contains("no3r"))
                            {
                                if (matrixLatLonNO3R_Model != null && matrixLatLonNO3R_Monitor != null && matrixNO3R_Model != null && matrixNO3R_Monitor != null)
                                {
                                    //SaveMatlabInputData(matrixLatLonNO3R_Model, matrixLatLonNO3R_Monitor, matrixNO3R_Model, matrixNO3R_Monitor, "no3r", _settings.Cmaqres.ToString(), q);
                                    result = SingleRun(matrixLatLonNO3R_Model, matrixLatLonNO3R_Monitor, matrixNO3R_Model, matrixNO3R_Monitor);
                                    if (result != null)
                                    {
                                        arrPrediction = result.Item1;
                                        arrSepred = result.Item2;
                                        //SaveDSOutputData(matrixLatLonNO3R_Model, arrPrediction, arrSepred, "no3r", _settings.Cmaqres.ToString(), q);
                                        foreach (KeyValuePair<string, AnnualPMDSOutput> kv in _dicDSOutputValue)
                                        {
                                            for (int j = 0; j < matrixLatLonNO3R_Model.RowCount; j++)
                                            {
                                                if (kv.Value.Lat == matrixLatLonNO3R_Model.At(j, 0) && kv.Value.Lon == matrixLatLonNO3R_Model.At(j, 1))
                                                {
                                                    if (kv.Value.DicQuarterlyDSValue.ContainsKey(q))
                                                    {
                                                        kv.Value.DicQuarterlyDSValue[q].NO3R_DS_Prediction = arrPrediction[j];
                                                        kv.Value.DicQuarterlyDSValue[q].NO3R_DS_Sepred = arrSepred[j];
                                                    }
                                                    break;
                                                }
                                            }
                                        }
                                    }
                                    GC.Collect();
                                }
                            }
                            #endregion
                            #region 计算OCB
                            if (!_modelMissingSpecies.Contains("oc"))
                            {
                                if (matrixLatLonOCB_Model != null && matrixLatLonOCB_Monitor != null && matrixOCB_Model != null && matrixOCB_Monitor != null)
                                {
                                    //SaveMatlabInputData(matrixLatLonOCB_Model, matrixLatLonOCB_Monitor, matrixOCB_Model, matrixOCB_Monitor, "ocb", _settings.Cmaqres.ToString(), q);
                                    result = SingleRun(matrixLatLonOCB_Model, matrixLatLonOCB_Monitor, matrixOCB_Model, matrixOCB_Monitor);
                                    if (result != null)
                                    {
                                        arrPrediction = result.Item1;
                                        arrSepred = result.Item2;
                                        //SaveDSOutputData(matrixLatLonOCB_Model, arrPrediction, arrSepred, "ocb", _settings.Cmaqres.ToString(), q);
                                        foreach (KeyValuePair<string, AnnualPMDSOutput> kv in _dicDSOutputValue)
                                        {
                                            for (int j = 0; j < matrixLatLonOCB_Model.RowCount; j++)
                                            {
                                                if (kv.Value.Lat == matrixLatLonOCB_Model.At(j, 0) && kv.Value.Lon == matrixLatLonOCB_Model.At(j, 1))
                                                {
                                                    if (kv.Value.DicQuarterlyDSValue.ContainsKey(q))
                                                    {
                                                        kv.Value.DicQuarterlyDSValue[q].OCB_DS_Prediction = arrPrediction[j];
                                                        kv.Value.DicQuarterlyDSValue[q].OCB_DS_Sepred = arrSepred[j];
                                                    }
                                                    break;
                                                }
                                            }
                                        }
                                    }
                                    GC.Collect();
                                }
                            }
                            #endregion
                            #region 计算EC
                            if (!_modelMissingSpecies.Contains("ec"))
                            {
                                if (matrixLatLonEC_Model != null && matrixLatLonEC_Monitor != null && matrixEC_Model != null && matrixEC_Monitor != null)
                                {
                                    //SaveMatlabInputData(matrixLatLonEC_Model, matrixLatLonEC_Monitor, matrixEC_Model, matrixEC_Monitor, "ec", _settings.Cmaqres.ToString(), q);
                                    result = SingleRun(matrixLatLonEC_Model, matrixLatLonEC_Monitor, matrixEC_Model, matrixEC_Monitor);
                                    if (result != null)
                                    {
                                        arrPrediction = result.Item1;
                                        arrSepred = result.Item2;
                                        //SaveDSOutputData(matrixLatLonEC_Model, arrPrediction, arrSepred, "ec", _settings.Cmaqres.ToString(), q);
                                        foreach (KeyValuePair<string, AnnualPMDSOutput> kv in _dicDSOutputValue)
                                        {
                                            for (int j = 0; j < matrixLatLonEC_Model.RowCount; j++)
                                            {
                                                if (kv.Value.Lat == matrixLatLonEC_Model.At(j, 0) && kv.Value.Lon == matrixLatLonEC_Model.At(j, 1))
                                                {
                                                    if (kv.Value.DicQuarterlyDSValue.ContainsKey(q))
                                                    {
                                                        kv.Value.DicQuarterlyDSValue[q].EC_DS_Prediction = arrPrediction[j];
                                                        kv.Value.DicQuarterlyDSValue[q].EC_DS_Sepred = arrSepred[j];
                                                    }
                                                    break;
                                                }
                                            }
                                        }
                                    }
                                    GC.Collect();
                                }
                            }
                            #endregion
                            #region 计算Crustal
                            if (!_modelMissingSpecies.Contains("crustal"))
                            {
                                if (matrixLatLonCrustal_Model != null && matrixLatLonCrustal_Monitor != null && matrixCrustal_Model != null && matrixCrustal_Monitor != null)
                                {
                                    //SaveMatlabInputData(matrixLatLonCrustal_Model, matrixLatLonCrustal_Monitor, matrixCrustal_Model, matrixCrustal_Monitor, "crustal", _settings.Cmaqres.ToString(), q);
                                    result = SingleRun(matrixLatLonCrustal_Model, matrixLatLonCrustal_Monitor, matrixCrustal_Model, matrixCrustal_Monitor);
                                    if (result != null)
                                    {
                                        arrPrediction = result.Item1;
                                        arrSepred = result.Item2;
                                        //SaveDSOutputData(matrixLatLonCrustal_Model, arrPrediction, arrSepred, "crustal", _settings.Cmaqres.ToString(), q);
                                        foreach (KeyValuePair<string, AnnualPMDSOutput> kv in _dicDSOutputValue)
                                        {
                                            for (int j = 0; j < matrixLatLonCrustal_Model.RowCount; j++)
                                            {
                                                if (kv.Value.Lat == matrixLatLonCrustal_Model.At(j, 0) && kv.Value.Lon == matrixLatLonCrustal_Model.At(j, 1))
                                                {
                                                    if (kv.Value.DicQuarterlyDSValue.ContainsKey(q))
                                                    {
                                                        kv.Value.DicQuarterlyDSValue[q].Crustal_DS_Prediction = arrPrediction[j];
                                                        kv.Value.DicQuarterlyDSValue[q].Crustal_DS_Sepred = arrSepred[j];
                                                    }
                                                    break;
                                                }
                                            }
                                        }
                                    }
                                    GC.Collect();
                                }
                            }
                            #endregion
                            #region 计算Salt
                            if (!_modelMissingSpecies.Contains("salt"))
                            {
                                if (matrixLatLonSalt_Model != null && matrixLatLonSalt_Monitor != null && matrixSalt_Model != null && matrixSalt_Monitor != null)
                                {
                                    //SaveMatlabInputData(matrixLatLonSalt_Model, matrixLatLonSalt_Monitor, matrixSalt_Model, matrixSalt_Monitor, "salt", _settings.Cmaqres.ToString(), q);
                                    result = SingleRun(matrixLatLonSalt_Model, matrixLatLonSalt_Monitor, matrixSalt_Model, matrixSalt_Monitor);
                                    if (result != null)
                                    {
                                        arrPrediction = result.Item1;
                                        arrSepred = result.Item2;
                                        //SaveDSOutputData(matrixLatLonSalt_Model, arrPrediction, arrSepred, "salt", _settings.Cmaqres.ToString(), q);
                                        foreach (KeyValuePair<string, AnnualPMDSOutput> kv in _dicDSOutputValue)
                                        {
                                            for (int j = 0; j < matrixLatLonSalt_Model.RowCount; j++)
                                            {
                                                if (kv.Value.Lat == matrixLatLonSalt_Model.At(j, 0) && kv.Value.Lon == matrixLatLonSalt_Model.At(j, 1))
                                                {
                                                    if (kv.Value.DicQuarterlyDSValue.ContainsKey(q))
                                                    {
                                                        kv.Value.DicQuarterlyDSValue[q].Salt_DS_Prediction = arrPrediction[j];
                                                        kv.Value.DicQuarterlyDSValue[q].Salt_DS_Sepred = arrSepred[j];
                                                    }
                                                    break;
                                                }
                                            }
                                        }
                                    }
                                    GC.Collect();
                                }
                            }
                            #endregion
                            #region 计算NO3
                            if (!_modelMissingSpecies.Contains("no3"))
                            {
                                if (matrixLatLonNO3_Model != null && matrixLatLonNO3_Monitor != null && matrixNO3_Model != null && matrixNO3_Monitor != null)
                                {
                                    //SaveMatlabInputData(matrixLatLonNO3_Model, matrixLatLonNO3_Monitor, matrixNO3_Model, matrixNO3_Monitor, "no3", _settings.Cmaqres.ToString(), q);
                                    result = SingleRun(matrixLatLonNO3_Model, matrixLatLonNO3_Monitor, matrixNO3_Model, matrixNO3_Monitor);
                                    if (result != null)
                                    {
                                        arrPrediction = result.Item1;
                                        arrSepred = result.Item2;
                                        //SaveDSOutputData(matrixLatLonNO3_Model, arrPrediction, arrSepred, "no3", _settings.Cmaqres.ToString(), q);
                                        foreach (KeyValuePair<string, AnnualPMDSOutput> kv in _dicDSOutputValue)
                                        {
                                            for (int j = 0; j < matrixLatLonNO3_Model.RowCount; j++)
                                            {
                                                if (kv.Value.Lat == matrixLatLonNO3_Model.At(j, 0) && kv.Value.Lon == matrixLatLonNO3_Model.At(j, 1))
                                                {
                                                    if (kv.Value.DicQuarterlyDSValue.ContainsKey(q))
                                                    {
                                                        kv.Value.DicQuarterlyDSValue[q].NO3_DS_Prediction = arrPrediction[j];
                                                        kv.Value.DicQuarterlyDSValue[q].NO3_DS_Sepred = arrSepred[j];
                                                    }
                                                    break;
                                                }
                                            }
                                        }
                                    }
                                    GC.Collect();
                                }
                            }
                            #endregion
                            #region 计算NH4
                            if (!_modelMissingSpecies.Contains("nh4"))
                            {
                                if (matrixLatLonNH4_Model != null && matrixLatLonNH4_Monitor != null && matrixNH4_Model != null && matrixNH4_Monitor != null)
                                {
                                    //SaveMatlabInputData(matrixLatLonNH4_Model, matrixLatLonNH4_Monitor, matrixNH4_Model, matrixNH4_Monitor, "nh4", _settings.Cmaqres.ToString(), q);
                                    result = SingleRun(matrixLatLonNH4_Model, matrixLatLonNH4_Monitor, matrixNH4_Model, matrixNH4_Monitor);
                                    if (result != null)
                                    {
                                        arrPrediction = result.Item1;
                                        arrSepred = result.Item2;
                                        //SaveDSOutputData(matrixLatLonNH4_Model, arrPrediction, arrSepred, "nh4", _settings.Cmaqres.ToString(), q);
                                        foreach (KeyValuePair<string, AnnualPMDSOutput> kv in _dicDSOutputValue)
                                        {
                                            for (int j = 0; j < matrixLatLonNH4_Model.RowCount; j++)
                                            {
                                                if (kv.Value.Lat == matrixLatLonNH4_Model.At(j, 0) && kv.Value.Lon == matrixLatLonNH4_Model.At(j, 1))
                                                {
                                                    if (kv.Value.DicQuarterlyDSValue.ContainsKey(q))
                                                    {
                                                        kv.Value.DicQuarterlyDSValue[q].NH4_DS_Prediction = arrPrediction[j];
                                                        kv.Value.DicQuarterlyDSValue[q].NH4_DS_Sepred = arrSepred[j];
                                                    }
                                                    break;
                                                }
                                            }
                                        }
                                    }
                                    GC.Collect();
                                }
                            }
                            #endregion
                            #region 计算DON，实际上没有算

                            #endregion
                        }
                        #endregion
                        GC.Collect();
                    }
                    matrixLatLonPM_Model = null; matrixLatLonPM_Monitor = null; matrixPM_Model = null; matrixPM_Monitor = null;
                    matrixLatLonSO4_Model = null; matrixLatLonSO4_Monitor = null; matrixSO4_Model = null; matrixSO4_Monitor = null;
                    matrixLatLonNO3R_Model = null; matrixLatLonNO3R_Monitor = null; matrixNO3R_Model = null; matrixNO3R_Monitor = null;
                    matrixLatLonOCB_Model = null; matrixLatLonOCB_Monitor = null; matrixOCB_Model = null; matrixOCB_Monitor = null;
                    matrixLatLonEC_Model = null; matrixLatLonEC_Monitor = null; matrixEC_Model = null; matrixEC_Monitor = null;
                    matrixLatLonCrustal_Model = null; matrixLatLonCrustal_Monitor = null; matrixCrustal_Model = null; matrixCrustal_Monitor = null;
                    matrixLatLonSalt_Model = null; matrixLatLonSalt_Monitor = null; matrixSalt_Model = null; matrixSalt_Monitor = null;
                    matrixLatLonNO3_Model = null; matrixLatLonNO3_Monitor = null; matrixNO3_Model = null; matrixNO3_Monitor = null;
                    matrixLatLonNH4_Model = null; matrixLatLonNH4_Monitor = null; matrixNH4_Model = null; matrixNH4_Monitor = null;
                    result = null; arrPrediction = new double[0]; arrSepred = new double[0];
                    GC.Collect();
                }
                #endregion

                return true;
            }
            catch (Exception ex)
            {
                _errormsg.AppendFormat("start annual pm analysis error: {0}.", ex.Message);
            }
            return false;
        }
        #region 垃圾回收
        /// <summary>
        /// 释放资源，全部清空
        /// </summary>
        public void Dispose()
        {
            Dispose(true);
            GC.SuppressFinalize(this);
        }
        /// <summary>
        /// 析构函数
        /// </summary>
        ~AnnualPMAnalysis()
        {
            Dispose(false);
        }
        protected virtual void Dispose(bool disposing)
        {
            if (!_disposed)
            {
                if (disposing)
                {
                    // Dispose the managed resources
                    _dicModelData.Clear(); _dicModelData = null;
                    _dicSpeciesMonitorData.Clear(); _dicSpeciesMonitorData = null;
                    _dicPmMonitorData.Clear(); _dicPmMonitorData = null;
                    _settings = null;
                    _errormsg.Clear(); _errormsg = null;
                    _dicDSOutputValue.Clear(); _dicDSOutputValue = null;
                    _modelMissingSpecies = new string[0];
                }
                // Dispose the unmanaged resources

            }
            _disposed = true;
        }
        #endregion
        #endregion

        #region 私有方法
        /// <summary>
        /// 根据季度处理输入参数，矩阵化
        /// </summary>
        private bool MatrixingInputData(
            string quarter,
            ref Dictionary<string, AnnualPMDSOutput> dicDSOutputValue,
            ref Matrix<double> matrixLatLonPM_Model, ref Matrix<double> matrixLatLonPM_Monitor, ref Matrix<double> matrixPM_Model, ref Matrix<double> matrixPM_Monitor,
            ref Matrix<double> matrixLatLonSO4_Model, ref Matrix<double> matrixLatLonSO4_Monitor, ref Matrix<double> matrixSO4_Model, ref Matrix<double> matrixSO4_Monitor,
            ref Matrix<double> matrixLatLonNO3R_Model, ref Matrix<double> matrixLatLonNO3R_Monitor, ref Matrix<double> matrixNO3R_Model, ref Matrix<double> matrixNO3R_Monitor,
            ref Matrix<double> matrixLatLonOCB_Model, ref Matrix<double> matrixLatLonOCB_Monitor, ref Matrix<double> matrixOCB_Model, ref Matrix<double> matrixOCB_Monitor,
            ref Matrix<double> matrixLatLonEC_Model, ref Matrix<double> matrixLatLonEC_Monitor, ref Matrix<double> matrixEC_Model, ref Matrix<double> matrixEC_Monitor,
            ref Matrix<double> matrixLatLonCrustal_Model, ref Matrix<double> matrixLatLonCrustal_Monitor, ref Matrix<double> matrixCrustal_Model, ref Matrix<double> matrixCrustal_Monitor,
            ref Matrix<double> matrixLatLonSalt_Model, ref Matrix<double> matrixLatLonSalt_Monitor, ref Matrix<double> matrixSalt_Model, ref Matrix<double> matrixSalt_Monitor,
            ref Matrix<double> matrixLatLonNO3_Model, ref Matrix<double> matrixLatLonNO3_Monitor, ref Matrix<double> matrixNO3_Model, ref Matrix<double> matrixNO3_Monitor,
            ref Matrix<double> matrixLatLonNH4_Model, ref Matrix<double> matrixLatLonNH4_Monitor, ref Matrix<double> matrixNH4_Model, ref Matrix<double> matrixNH4_Monitor)
        {
            try
            {
                #region 一些中间变量
                List<double> lst_LatLon_PM = new List<double>();
                List<double> lst_PM = new List<double>();
                List<double> lst_LatLon_SO4 = new List<double>();
                List<double> lst_SO4 = new List<double>();
                List<double> lst_LatLon_NO3R = new List<double>();
                List<double> lst_NO3R = new List<double>();
                List<double> lst_LatLon_OCB = new List<double>();
                List<double> lst_OCB = new List<double>();
                List<double> lst_LatLon_EC = new List<double>();
                List<double> lst_EC = new List<double>();
                List<double> lst_LatLon_Crustal = new List<double>();
                List<double> lst_Crustal = new List<double>();
                List<double> lst_LatLon_Salt = new List<double>();
                List<double> lst_Salt = new List<double>();
                List<double> lst_LatLon_NO3 = new List<double>();
                List<double> lst_NO3 = new List<double>();
                List<double> lst_LatLon_NH4 = new List<double>();
                List<double> lst_NH4 = new List<double>();
                #endregion 一些中间变量

                #region 处理Model数据
                Dictionary<string, ModelDataSpecies> dicModelQuarter = null;
                ModelDataSpecies modeldata = null;
                foreach (KeyValuePair<string, PmModelDataOutput> kv in _dicModelData)
                {
                    if (!dicDSOutputValue.ContainsKey(kv.Key))
                    {
                        dicDSOutputValue.Add(kv.Key, new AnnualPMDSOutput()
                        {
                            Id = kv.Key,
                            Lat = kv.Value.lat,
                            Lon = kv.Value.longitude,
                            DicQuarterlyDSValue = new Dictionary<string, AnnualPMDSQuarterlyValue>()
                        });
                    }
                    dicModelQuarter = kv.Value.dicBaselineModel;
                    if (dicModelQuarter == null || dicModelQuarter.Count <= 0 || !dicModelQuarter.ContainsKey(quarter))
                        continue;
                    modeldata = dicModelQuarter[quarter];
                    if (!dicDSOutputValue[kv.Key].DicQuarterlyDSValue.ContainsKey(quarter))
                    {
                        dicDSOutputValue[kv.Key].DicQuarterlyDSValue.Add(quarter, new AnnualPMDSQuarterlyValue()
                        {
                            SO4_Base = modeldata.so4,
                            OCB_Base = modeldata.oc,
                            EC_Base = modeldata.ec,
                            Crustal_Base = modeldata.crustal,
                            Salt_Base = modeldata.salt,
                            NO3_Base = modeldata.no3,
                            NH4_Base = modeldata.nh4,
                            PM25_Base = modeldata.pm25
                        });
                    }
                    #region 过滤各类数据
                    // PM25
                    if (modeldata.pm25 >= 0)
                    {
                        lst_LatLon_PM.Add(kv.Value.lat);
                        lst_LatLon_PM.Add(kv.Value.longitude);
                        lst_PM.Add(modeldata.pm25 == 0 ? 0.1 : modeldata.pm25);
                    }
                    // SO4
                    if (modeldata.so4 >= 0)
                    {
                        lst_LatLon_SO4.Add(kv.Value.lat);
                        lst_LatLon_SO4.Add(kv.Value.longitude);
                        lst_SO4.Add(modeldata.so4 == 0 ? 0.1 : modeldata.so4);
                    }
                    // NO3和NO3R
                    if (modeldata.no3 >= 0)
                    {
                        lst_LatLon_NO3.Add(kv.Value.lat);
                        lst_LatLon_NO3.Add(kv.Value.longitude);
                        double v = modeldata.no3 == 0 ? 0.1 : modeldata.no3;
                        lst_NO3.Add(v);
                        lst_LatLon_NO3R.Add(kv.Value.lat);
                        lst_LatLon_NO3R.Add(kv.Value.longitude);
                        lst_NO3R.Add(v);
                    }
                    // OCB
                    if (modeldata.oc >= 0)
                    {
                        lst_LatLon_OCB.Add(kv.Value.lat);
                        lst_LatLon_OCB.Add(kv.Value.longitude);
                        lst_OCB.Add(modeldata.oc == 0 ? 0.1 : modeldata.oc);
                    }
                    // EC
                    if (modeldata.ec >= 0)
                    {
                        lst_LatLon_EC.Add(kv.Value.lat);
                        lst_LatLon_EC.Add(kv.Value.longitude);
                        lst_EC.Add(modeldata.ec == 0 ? 0.1 : modeldata.ec);
                    }
                    // Crustal
                    if (modeldata.crustal >= 0)
                    {
                        lst_LatLon_Crustal.Add(kv.Value.lat);
                        lst_LatLon_Crustal.Add(kv.Value.longitude);
                        lst_Crustal.Add(modeldata.crustal == 0 ? 0.1 : modeldata.crustal);
                    }
                    // Salt
                    if (modeldata.salt >= 0)
                    {
                        lst_LatLon_Salt.Add(kv.Value.lat);
                        lst_LatLon_Salt.Add(kv.Value.longitude);
                        lst_Salt.Add(modeldata.salt == 0 ? 0.1 : modeldata.salt);
                    }
                    // NH4
                    if (modeldata.nh4 >= 0)
                    {
                        lst_LatLon_NH4.Add(kv.Value.lat);
                        lst_LatLon_NH4.Add(kv.Value.longitude);
                        lst_NH4.Add(modeldata.nh4 == 0 ? 0.1 : modeldata.nh4);
                    }
                    #endregion 过滤各类数据
                }// foreach
                #region 矩阵化
                // PM25
                if (lst_PM.Count > 0)
                {
                    matrixLatLonPM_Model = Matrix<double>.Build.DenseOfColumnMajor(2, lst_LatLon_PM.Count / 2, lst_LatLon_PM.ToArray());
                    matrixLatLonPM_Model = matrixLatLonPM_Model.Transpose();
                    matrixPM_Model = Matrix<double>.Build.DenseOfColumnMajor(lst_PM.Count, 1, lst_PM.ToArray());
                }
                // SO4
                if (lst_SO4.Count > 0)
                {
                    matrixLatLonSO4_Model = Matrix<double>.Build.DenseOfColumnMajor(2, lst_LatLon_SO4.Count / 2, lst_LatLon_SO4.ToArray());
                    matrixLatLonSO4_Model = matrixLatLonSO4_Model.Transpose();
                    matrixSO4_Model = Matrix<double>.Build.DenseOfColumnMajor(lst_SO4.Count, 1, lst_SO4.ToArray());
                }
                // NO3R
                if (lst_NO3R.Count > 0)
                {
                    matrixLatLonNO3R_Model = Matrix<double>.Build.DenseOfColumnMajor(2, lst_LatLon_NO3R.Count / 2, lst_LatLon_NO3R.ToArray());
                    matrixLatLonNO3R_Model = matrixLatLonNO3R_Model.Transpose();
                    matrixNO3R_Model = Matrix<double>.Build.DenseOfColumnMajor(lst_NO3R.Count, 1, lst_NO3R.ToArray());
                }
                // OCB
                if (lst_OCB.Count > 0)
                {
                    matrixLatLonOCB_Model = Matrix<double>.Build.DenseOfColumnMajor(2, lst_LatLon_OCB.Count / 2, lst_LatLon_OCB.ToArray());
                    matrixLatLonOCB_Model = matrixLatLonOCB_Model.Transpose();
                    matrixOCB_Model = Matrix<double>.Build.DenseOfColumnMajor(lst_OCB.Count, 1, lst_OCB.ToArray());
                }
                // EC
                if (lst_EC.Count > 0)
                {
                    matrixLatLonEC_Model = Matrix<double>.Build.DenseOfColumnMajor(2, lst_LatLon_EC.Count / 2, lst_LatLon_EC.ToArray());
                    matrixLatLonEC_Model = matrixLatLonEC_Model.Transpose();
                    matrixEC_Model = Matrix<double>.Build.DenseOfColumnMajor(lst_EC.Count, 1, lst_EC.ToArray());
                }
                // Crustal
                if (lst_Crustal.Count > 0)
                {
                    matrixLatLonCrustal_Model = Matrix<double>.Build.DenseOfColumnMajor(2, lst_LatLon_Crustal.Count / 2, lst_LatLon_Crustal.ToArray());
                    matrixLatLonCrustal_Model = matrixLatLonCrustal_Model.Transpose();
                    matrixCrustal_Model = Matrix<double>.Build.DenseOfColumnMajor(lst_Crustal.Count, 1, lst_Crustal.ToArray());
                }
                // Salt
                if (lst_Salt.Count > 0)
                {
                    matrixLatLonSalt_Model = Matrix<double>.Build.DenseOfColumnMajor(2, lst_LatLon_Salt.Count / 2, lst_LatLon_Salt.ToArray());
                    matrixLatLonSalt_Model = matrixLatLonSalt_Model.Transpose();
                    matrixSalt_Model = Matrix<double>.Build.DenseOfColumnMajor(lst_Salt.Count, 1, lst_Salt.ToArray());
                }
                // NO3
                if (lst_NO3.Count > 0)
                {
                    matrixLatLonNO3_Model = Matrix<double>.Build.DenseOfColumnMajor(2, lst_LatLon_NO3.Count / 2, lst_LatLon_NO3.ToArray());
                    matrixLatLonNO3_Model = matrixLatLonNO3_Model.Transpose();
                    matrixNO3_Model = Matrix<double>.Build.DenseOfColumnMajor(lst_NO3.Count, 1, lst_NO3.ToArray());
                }
                // NH4
                if (lst_NH4.Count > 0)
                {
                    matrixLatLonNH4_Model = Matrix<double>.Build.DenseOfColumnMajor(2, lst_LatLon_NH4.Count / 2, lst_LatLon_NH4.ToArray());
                    matrixLatLonNH4_Model = matrixLatLonNH4_Model.Transpose();
                    matrixNH4_Model = Matrix<double>.Build.DenseOfColumnMajor(lst_NH4.Count, 1, lst_NH4.ToArray());
                }
                #endregion 矩阵化
                #endregion 处理Model数据

                #region 处理Monitor PM数据
                lst_LatLon_PM.Clear();
                lst_PM.Clear();
                lst_LatLon_SO4.Clear();
                lst_SO4.Clear();
                lst_LatLon_NO3R.Clear();
                lst_NO3R.Clear();
                lst_LatLon_OCB.Clear();
                lst_OCB.Clear();
                lst_LatLon_EC.Clear();
                lst_EC.Clear();
                lst_LatLon_Crustal.Clear();
                lst_Crustal.Clear();
                lst_LatLon_Salt.Clear();
                lst_Salt.Clear();
                lst_LatLon_NO3.Clear();
                lst_NO3.Clear();
                lst_LatLon_NH4.Clear();
                lst_NH4.Clear();

                Dictionary<string, float> dicMonitorPMQuarter = null;
                foreach (KeyValuePair<string, PMMonitorOutput> kv in _dicPmMonitorData)
                {
                    dicMonitorPMQuarter = kv.Value.dicUnofficialPM;
                    if (dicMonitorPMQuarter == null || dicMonitorPMQuarter.Count <= 0 || !dicMonitorPMQuarter.ContainsKey(quarter))
                        continue;
                    // PM25
                    float v = dicMonitorPMQuarter[quarter];
                    if (v >= 0)
                    {
                        lst_LatLon_PM.Add(kv.Value.latitude);
                        lst_LatLon_PM.Add(kv.Value.longitude);
                        lst_PM.Add(v == 0 ? 0.1 : v);
                    }
                }
                // PM25
                if (lst_PM.Count > 0)
                {
                    matrixLatLonPM_Monitor = Matrix<double>.Build.DenseOfColumnMajor(2, lst_LatLon_PM.Count / 2, lst_LatLon_PM.ToArray());
                    matrixLatLonPM_Monitor = matrixLatLonPM_Monitor.Transpose();
                    matrixPM_Monitor = Matrix<double>.Build.DenseOfColumnMajor(lst_PM.Count, 1, lst_PM.ToArray());
                }
                #endregion 处理Monitor PM数据

                #region 处理Monitor成分数据
                lst_LatLon_PM.Clear();
                lst_PM.Clear();
                lst_LatLon_SO4.Clear();
                lst_SO4.Clear();
                lst_LatLon_NO3R.Clear();
                lst_NO3R.Clear();
                lst_LatLon_OCB.Clear();
                lst_OCB.Clear();
                lst_LatLon_EC.Clear();
                lst_EC.Clear();
                lst_LatLon_Crustal.Clear();
                lst_Crustal.Clear();
                lst_LatLon_Salt.Clear();
                lst_Salt.Clear();
                lst_LatLon_NO3.Clear();
                lst_NO3.Clear();
                lst_LatLon_NH4.Clear();
                lst_NH4.Clear();

                Dictionary<string, SpeciesMonitors> dicMonitorSpeciesQuarter = null;
                SpeciesMonitors monitordata = null;
                foreach (KeyValuePair<string, PMMonitorOutput> kv in _dicSpeciesMonitorData)
                {
                    dicMonitorSpeciesQuarter = kv.Value.dicSpecies;
                    if (dicMonitorSpeciesQuarter == null || dicMonitorSpeciesQuarter.Count <= 0 || !dicMonitorSpeciesQuarter.ContainsKey(quarter))
                        continue;
                    monitordata = dicMonitorSpeciesQuarter[quarter];
                    #region 过滤各类数据
                    // SO4
                    if (float.IsNaN(monitordata.bSO4Mass)) monitordata.bSO4Mass = 0.001f;
                    if (monitordata.bSO4Mass >= 0)
                    {
                        lst_LatLon_SO4.Add(kv.Value.latitude);
                        lst_LatLon_SO4.Add(kv.Value.longitude);
                        lst_SO4.Add(monitordata.bSO4Mass == 0 ? 0.1 : monitordata.bSO4Mass);
                    }
                    // NO3
                    if (float.IsNaN(monitordata.bNO3Mass)) monitordata.bNO3Mass = 0.001f;
                    if (monitordata.bNO3Mass >= 0)
                    {
                        lst_LatLon_NO3.Add(kv.Value.latitude);
                        lst_LatLon_NO3.Add(kv.Value.longitude);
                        lst_NO3.Add(monitordata.bNO3Mass == 0 ? 0.1 : monitordata.bNO3Mass);
                    }
                    // NO3R
                    if (float.IsNaN(monitordata.bNO3rMass)) monitordata.bNO3rMass = 0.001f;
                    if (monitordata.bNO3rMass >= 0)
                    {
                        lst_LatLon_NO3R.Add(kv.Value.latitude);
                        lst_LatLon_NO3R.Add(kv.Value.longitude);
                        lst_NO3R.Add(monitordata.bNO3rMass == 0 ? 0.1 : monitordata.bNO3rMass);
                    }
                    // OCB
                    if (float.IsNaN(monitordata.bOCBMass)) monitordata.bOCBMass = 0.001f;
                    if (monitordata.bOCBMass >= 0)
                    {
                        lst_LatLon_OCB.Add(kv.Value.latitude);
                        lst_LatLon_OCB.Add(kv.Value.longitude);
                        lst_OCB.Add(monitordata.bOCBMass == 0 ? 0.1 : monitordata.bOCBMass);
                    }
                    // EC
                    if (float.IsNaN(monitordata.bECMass)) monitordata.bECMass = 0.001f;
                    if (monitordata.bECMass >= 0)
                    {
                        lst_LatLon_EC.Add(kv.Value.latitude);
                        lst_LatLon_EC.Add(kv.Value.longitude);
                        lst_EC.Add(monitordata.bECMass == 0 ? 0.1 : monitordata.bECMass);
                    }
                    // Crustal
                    if (float.IsNaN(monitordata.bCrustalMass)) monitordata.bCrustalMass = 0.001f;
                    if (monitordata.bCrustalMass >= 0)
                    {
                        lst_LatLon_Crustal.Add(kv.Value.latitude);
                        lst_LatLon_Crustal.Add(kv.Value.longitude);
                        lst_Crustal.Add(monitordata.bCrustalMass == 0 ? 0.1 : monitordata.bCrustalMass);
                    }
                    // Salt
                    if (float.IsNaN(monitordata.bSaltMass)) monitordata.bSaltMass = 0.001f;
                    if (monitordata.bSaltMass >= 0)
                    {
                        lst_LatLon_Salt.Add(kv.Value.latitude);
                        lst_LatLon_Salt.Add(kv.Value.longitude);
                        lst_Salt.Add(monitordata.bSaltMass == 0 ? 0.1 : monitordata.bSaltMass);
                    }
                    // NH4
                    if (float.IsNaN(monitordata.bNH4Mass)) monitordata.bNH4Mass = 0.001f;
                    if (monitordata.bNH4Mass >= 0)
                    {
                        lst_LatLon_NH4.Add(kv.Value.latitude);
                        lst_LatLon_NH4.Add(kv.Value.longitude);
                        lst_NH4.Add(monitordata.bNH4Mass == 0 ? 0.1 : monitordata.bNH4Mass);
                    }
                    #endregion 过滤各类数据
                }
                #region 矩阵化
                // SO4
                if (lst_SO4.Count > 0)
                {
                    matrixLatLonSO4_Monitor = Matrix<double>.Build.DenseOfColumnMajor(2, lst_LatLon_SO4.Count / 2, lst_LatLon_SO4.ToArray());
                    matrixLatLonSO4_Monitor = matrixLatLonSO4_Monitor.Transpose();
                    matrixSO4_Monitor = Matrix<double>.Build.DenseOfColumnMajor(lst_SO4.Count, 1, lst_SO4.ToArray());
                }
                // NO3R
                if (lst_NO3R.Count > 0)
                {
                    matrixLatLonNO3R_Monitor = Matrix<double>.Build.DenseOfColumnMajor(2, lst_LatLon_NO3R.Count / 2, lst_LatLon_NO3R.ToArray());
                    matrixLatLonNO3R_Monitor = matrixLatLonNO3R_Monitor.Transpose();
                    matrixNO3R_Monitor = Matrix<double>.Build.DenseOfColumnMajor(lst_NO3R.Count, 1, lst_NO3R.ToArray());
                }
                // OCB
                if (lst_OCB.Count > 0)
                {
                    matrixLatLonOCB_Monitor = Matrix<double>.Build.DenseOfColumnMajor(2, lst_LatLon_OCB.Count / 2, lst_LatLon_OCB.ToArray());
                    matrixLatLonOCB_Monitor = matrixLatLonOCB_Monitor.Transpose();
                    matrixOCB_Monitor = Matrix<double>.Build.DenseOfColumnMajor(lst_OCB.Count, 1, lst_OCB.ToArray());
                }
                // EC
                if (lst_EC.Count > 0)
                {
                    matrixLatLonEC_Monitor = Matrix<double>.Build.DenseOfColumnMajor(2, lst_LatLon_EC.Count / 2, lst_LatLon_EC.ToArray());
                    matrixLatLonEC_Monitor = matrixLatLonEC_Monitor.Transpose();
                    matrixEC_Monitor = Matrix<double>.Build.DenseOfColumnMajor(lst_EC.Count, 1, lst_EC.ToArray());
                }
                // Crustal
                if (lst_Crustal.Count > 0)
                {
                    matrixLatLonCrustal_Monitor = Matrix<double>.Build.DenseOfColumnMajor(2, lst_LatLon_Crustal.Count / 2, lst_LatLon_Crustal.ToArray());
                    matrixLatLonCrustal_Monitor = matrixLatLonCrustal_Monitor.Transpose();
                    matrixCrustal_Monitor = Matrix<double>.Build.DenseOfColumnMajor(lst_Crustal.Count, 1, lst_Crustal.ToArray());
                }
                // Salt
                if (lst_Salt.Count > 0)
                {
                    matrixLatLonSalt_Monitor = Matrix<double>.Build.DenseOfColumnMajor(2, lst_LatLon_Salt.Count / 2, lst_LatLon_Salt.ToArray());
                    matrixLatLonSalt_Monitor = matrixLatLonSalt_Monitor.Transpose();
                    matrixSalt_Monitor = Matrix<double>.Build.DenseOfColumnMajor(lst_Salt.Count, 1, lst_Salt.ToArray());
                }
                // NO3
                if (lst_NO3.Count > 0)
                {
                    matrixLatLonNO3_Monitor = Matrix<double>.Build.DenseOfColumnMajor(2, lst_LatLon_NO3.Count / 2, lst_LatLon_NO3.ToArray());
                    matrixLatLonNO3_Monitor = matrixLatLonNO3_Monitor.Transpose();
                    matrixNO3_Monitor = Matrix<double>.Build.DenseOfColumnMajor(lst_NO3.Count, 1, lst_NO3.ToArray());
                }
                // NH4
                if (lst_NH4.Count > 0)
                {
                    matrixLatLonNH4_Monitor = Matrix<double>.Build.DenseOfColumnMajor(2, lst_LatLon_NH4.Count / 2, lst_LatLon_NH4.ToArray());
                    matrixLatLonNH4_Monitor = matrixLatLonNH4_Monitor.Transpose();
                    matrixNH4_Monitor = Matrix<double>.Build.DenseOfColumnMajor(lst_NH4.Count, 1, lst_NH4.ToArray());
                }
                #endregion 矩阵化
                #endregion 处理Monitor成分数据

                lst_LatLon_PM.Clear();
                lst_PM.Clear();
                lst_LatLon_SO4.Clear();
                lst_SO4.Clear();
                lst_LatLon_NO3R.Clear();
                lst_NO3R.Clear();
                lst_LatLon_OCB.Clear();
                lst_OCB.Clear();
                lst_LatLon_EC.Clear();
                lst_EC.Clear();
                lst_LatLon_Crustal.Clear();
                lst_Crustal.Clear();
                lst_LatLon_Salt.Clear();
                lst_Salt.Clear();
                lst_LatLon_NO3.Clear();
                lst_NO3.Clear();
                lst_LatLon_NH4.Clear();
                lst_NH4.Clear();
                return true;
            }
            catch (Exception ex)
            {
                _errormsg.AppendFormat("matrixing smat-ce input datas error: {0}.{1}", ex.Message, Environment.NewLine);
            }
            return false;
        }
        /// <summary>
        /// 单线程运行
        /// </summary>
        /// <param name="matrixLatLon_Model">Model经纬度矩阵</param>
        /// <param name="matrixLatLon_Monitor">Monitor经纬度矩阵</param>
        /// <param name="matrix_Model">Model控制值矩阵</param>
        /// <param name="matrix_Monitor">Monitor控制值矩阵</param>
        /// <returns>返回DS运算后的Prediction和Sepred，按顺序</returns>
        private Tuple<double[], double[]> SingleRun(Matrix<double> matrixLatLon_Model, Matrix<double> matrixLatLon_Monitor, Matrix<double> matrix_Model, Matrix<double> matrix_Monitor)
        {
            try
            {
                #region 计算步骤
                #region 1. 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 2. spatialize model matrix
                double vct2 = MathNet.Numerics.Constants.Pi / 180;
                data_d = matrixLatLon_Model.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(matrixLatLon_Model.RowCount, matrixLatLon_Model.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();
                var qlat = Matrix<double>.Build.DenseOfColumnMajor(20, 1, data_c1);
                var qlon = Matrix<double>.Build.DenseOfColumnMajor(20, 1, data_c2);

                SpecialFunction.SpatializeMatrix(ref s1);

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

                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);

                data_d = new double[0]; data_d_clone = new double[0]; data_c1 = new double[0]; data_c2 = new double[0];
                GC.Collect();
                #endregion

                #region 3. at 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;
                for (int i = 0; i < len; i++)
                {
                    if (data_d[i] > 0.1)
                    {
                        lstConcRow.Add(s2.Row(i));
                        lstLatlonRow.Add(s3.Row(i));
                    }
                }
                s2 = Matrix<double>.Build.DenseOfRowVectors(lstConcRow);
                s3 = Matrix<double>.Build.SparseOfRowVectors(lstLatlonRow);

                data_d = new double[0]; lstConcRow.Clear(); lstLatlonRow.Clear();
                GC.Collect();
                #endregion

                #region 4. get the matrix "mdta"
                data_d = matrix_Monitor.ToColumnWiseArray();
                var tuple = MatrixCompute.MatrixRowIsmemberMatrixRow(matrixLatLon_Monitor, 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();
                        temp = new double[0]; wgt = new double[0];
                    }
                    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);
                    }
                    index_arr = new int[0]; item2 = null; match = new double[0];
                }

                data_d = new double[0]; data_d_clone = new double[0]; data_c1 = new double[0]; data_c2 = new double[0]; data_c3 = new double[0];
                GC.Collect();
                #endregion

                #region 5. 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 6. 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);
                lstLatlonRow.Clear();
                var ds = MatrixCompute.Pdist2(s2_row, s2_row);
                double[] tmp = new double[mdta.RowCount];
                double max;
                data_d = matrix_Model.Column(0).ToArray();
                double[] data_d2;
                for (int i = 0; i < data_c2.Length; i++)
                {
                    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;
                }
                data_d2 = new double[0];
                LinearRegression sl = MatrixCompute.Regstats(tmp, data_c1);
                if (sl.ResidualVariance == 0.0)
                    sl.ResidualVariance = 0.01;
                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];
                }
                tmp = new double[0];
                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;

                    matrixEv = new double[0]; vectorEv = new System.Numerics.Complex[0]; matrixD = new double[0];
                    tmp2 = new double[0]; tmp3 = new double[0];
                }// for 9
                double phi_b0s = 3.0 / ((Array.IndexOf(probs, probs.Max()) * 0.1 + 0.1) * max);
                Debug.WriteLine("phi_b0s = {0}", phi_b0s);

                probs = new double[0]; data_d = new double[0]; tmp1 = new double[0];
                b0sinvmat_d = new double[0];
                GC.Collect();
                #endregion

                #region 7. 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 8. 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);
                lstLatlonRow.Clear();
                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);
                s3_row2 = null;
                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);
                s3_row = null;
                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, matrix_Model);
                ind = new bool[0];
                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();

                data_d = new double[0]; W_data = new double[0];
                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 9. 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 = matrix_Model.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();
                    }
                }// for
                Lt = null; Xtp = null;
                prop = new double[0]; pmean = new double[0]; pvar = new double[0]; kY = new double[0]; kY2 = new double[0];
                data_Xt = new double[0]; data_Xt_trans = new double[0]; data_XtXt = new double[0]; data_b0s = new double[0]; data_c1_b0s = new double[0];
                data_s4 = new double[0]; data_matrixModelConc_c1 = new double[0];
                lst2.Clear();
                #endregion
                #endregion
                double[] result = SpecialFunction.Compute_Result(kp, Y2bar, Ybar);
                return Tuple.Create<double[], double[]>(Ybar, result);
            }
            catch (Exception ex)
            {
                _errormsg.AppendFormat("single run ds error: {0}.{1}", ex.Message, Environment.NewLine);
            }
            return null;
        }
        [Obsolete("保存Matlab输入格式数据，只做数据对比用")]
        private void SaveMatlabInputData(Matrix<double> matrixLatLon_Model, Matrix<double> matrixLatLon_Monitor, Matrix<double> matrix_Model, Matrix<double> matrix_Monitor, string name, string cmaqres, string q)
        {
            try
            {
                int count;
                FileStream fs;
                StreamWriter streamWriter;
                StringBuilder sb = new StringBuilder();
                string filepath = @"D:\SMAT_InputData\Annual_PM\";
                if (!Directory.Exists(filepath))
                    Directory.CreateDirectory(filepath);

                // 监测数据
                string filename = "monitor_" + name + "_" + cmaqres + "km_" + q + ".csv";
                filename = filepath + filename;
                if (File.Exists(filename))
                    File.Delete(filename);
                fs = new FileStream(filename, FileMode.Create, FileAccess.Write);
                streamWriter = new StreamWriter(fs, Encoding.Default);
                sb.Append("Site,POC,Date,Lat,Lon,Conc");
                streamWriter.WriteLine(sb.ToString());
                count = matrixLatLon_Monitor.RowCount;
                for (int i = 0; i < count; i++)
                {
                    sb.Clear();
                    sb.AppendFormat("{0},{1},{2},{3},{4},{5}", new object[] { i, i, "2011-1-1", matrixLatLon_Monitor.At(i, 0), matrixLatLon_Monitor.At(i, 1), matrix_Monitor.At(i, 0) });
                    streamWriter.WriteLine(sb.ToString());
                }
                streamWriter.Close();
                fs.Close();
                sb.Clear();

                // 模型数据
                filename = "model_" + name + "_" + cmaqres + "km_" + q + ".csv";
                filename = filepath + filename;
                if (File.Exists(filename))
                    File.Delete(filename);
                fs = new FileStream(filename, FileMode.Create, FileAccess.Write);
                streamWriter = new StreamWriter(fs, Encoding.Default);
                sb.Append("Date,Lat,Lon,Conc,id");
                streamWriter.WriteLine(sb.ToString());
                count = matrixLatLon_Model.RowCount;
                for (int i = 0; i < count; i++)
                {
                    sb.Clear();
                    sb.AppendFormat("{0},{1},{2},{3},{4}", new object[] { "2011-1-1", matrixLatLon_Model.At(i, 0), matrixLatLon_Model.At(i, 1), matrix_Model.At(i, 0), i });
                    streamWriter.WriteLine(sb.ToString());
                }
                streamWriter.Close();
                fs.Close();
                sb.Clear();
                streamWriter.Dispose();
                fs.Dispose();
                GC.Collect();
            }
            catch (Exception ex)
            {
                throw ex;
            }
        }
        [Obsolete("保存DS输出数据，只做数据对比用")]
        private void SaveDSOutputData(Matrix<double> matrixLatLon_Model, double[] arrPrediction, double[] arrSepred, string name, string cmaqres, string q)
        {
            try
            {
                string filename, filepath;
                filename = name + "_DS_" + cmaqres + "km_" + q + ".csv";
                filepath = AppDomain.CurrentDomain.BaseDirectory + @"\Result\";
                if (!Directory.Exists(filepath))
                    Directory.CreateDirectory(filepath);
                filename = filepath + filename;
                if (File.Exists(filename))
                    File.Delete(filename);
                FileStream fs = new FileStream(filename, 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());
                int count = matrixLatLon_Model.RowCount;
                for (int i = 0; i < count; i++)
                {
                    sb.Clear();
                    sb.AppendFormat("{0},{1},{2},{3}", new object[] { matrixLatLon_Model.At(i, 0), matrixLatLon_Model.At(i, 1), arrPrediction[i], arrSepred[i] });
                    streamWriter.WriteLine(sb.ToString());
                }
                streamWriter.Close();
                fs.Close();
                sb.Clear();
                streamWriter.Dispose();
                fs.Dispose();
                GC.Collect();
            }
            catch (Exception ex)
            {
                throw ex;
            }
        }
        #endregion
    }
}
