using System;
using System.Collections.Generic;
using System.ComponentModel;
using System.Data;
using System.Drawing;
using System.Linq;
using System.Text;
using System.Windows.Forms;

namespace vcsSIMPLEX
{
    public partial class Form1 : Form
    {
        struct DRAWPARA
        {          //作図用データ構造体
            public int kpt;       //拡大倍率（画面表示：1，ファイル保存：2）
            public string sxjiku; //ｘ軸名
            public string syjiku; //ｙ軸名
            public double xmin;   //ｘ軸最小値
            public double xmax;   //ｘ軸最大値
            public double dx;     //ｘ軸増分
            public double ymin;   //ｙ軸最小値
            public double ymax;   //ｙ軸最大値
            public double dy;     //ｙ軸増分
            public int kxi0;      //ｘ軸最小値座標（px）
            public int kxf0;      //ｘ軸最大値座標（px）
            public int kyi0;      //ｙ軸最小値座標（px）
            public int kyf0;      //ｙ軸最大値座標（px）
            public int HSIZE;     //画像横寸法
            public int DHL;       //画像左余白（縦軸描画スペース：px）
            public int DHR;       //画像右余白（px）
            public int VSIZE;     //画像縦サイズ（px）
            public int DVL;       //画像下余白（横軸描画スペース：px）
            public int DVU;       //画像上余白（px）
        }

        public Form1()
        {
            InitializeComponent();
        }
        //--------------------------------------------------------------------------------

        private void toolStripButton1_Click(object sender, EventArgs e)
        {
            main_part();
        }
        //--------------------------------------------------------------------------------
        private void main_part()
        {
            int i;
            string fname = "";
            string fnameR = ""; //入力ファイル名
            string fnameW = ""; //出力ファイル名
            string fnameP = ""; //png画像出力ファイル名
            System.IO.StreamReader sr;
            System.IO.StreamWriter sw;
            string dat;
            string[] sbuf;
            char delim = ',';
            int NDATA;
            string scom;
            DRAWPARA pltdata;

            int iflag;
            int iteration;
            int MPARA;

            INI_PLTDATA(out pltdata);

            MPARA = MPARA_SET(); //回帰係数の数
            double[] parameter = new double[MPARA];

            openFileDialog1.InitialDirectory = System.IO.Directory.GetCurrentDirectory();
            if (openFileDialog1.ShowDialog() == DialogResult.OK) { fnameR = openFileDialog1.FileName; }
            if (saveFileDialog1.ShowDialog() == DialogResult.OK) { fnameW = saveFileDialog1.FileName; }
            fname = System.IO.Path.GetFileNameWithoutExtension(fnameW);
            fnameP = System.IO.Path.GetDirectoryName(fnameW) + "\\" + fname + ".png";

            //データ入力
            sr = new System.IO.StreamReader(fnameR, System.Text.Encoding.Default);
            dat = sr.ReadLine(); sbuf = dat.Split(delim); scom = sbuf[0];
            dat = sr.ReadLine(); sbuf = dat.Split(delim); pltdata.sxjiku = sbuf[0]; pltdata.syjiku = sbuf[1];
            dat = sr.ReadLine(); sbuf = dat.Split(delim); pltdata.HSIZE = int.Parse(sbuf[0]); pltdata.DHL = int.Parse(sbuf[1]); pltdata.DHR = int.Parse(sbuf[2]);
            dat = sr.ReadLine(); sbuf = dat.Split(delim); pltdata.VSIZE = int.Parse(sbuf[0]); pltdata.DVL = int.Parse(sbuf[1]); pltdata.DVU = int.Parse(sbuf[2]);
            dat = sr.ReadLine(); sbuf = dat.Split(delim); pltdata.xmin = double.Parse(sbuf[0]); pltdata.xmax = double.Parse(sbuf[1]); pltdata.dx = double.Parse(sbuf[2]);
            dat = sr.ReadLine(); sbuf = dat.Split(delim); pltdata.ymin = double.Parse(sbuf[0]); pltdata.ymax = double.Parse(sbuf[1]); pltdata.dy = double.Parse(sbuf[2]);
            dat = sr.ReadLine(); sbuf = dat.Split(delim); NDATA = int.Parse(sbuf[0]) - 1;
            double[] datax = new double[NDATA + 1];
            double[] datay = new double[NDATA + 1];
            for (i = 0; i <= NDATA; i++)
            {
                dat = sr.ReadLine(); sbuf = dat.Split(delim);
                datax[i] = double.Parse(sbuf[0]);
                datay[i] = double.Parse(sbuf[1]);
            }
            sr.Close();

            //回帰係数計算
            SIMPLEX(MPARA, parameter, NDATA, datax, datay, out iteration, out iflag);

            //結果出力
            sw = new System.IO.StreamWriter(fnameW, false, System.Text.Encoding.Default);
            sw.WriteLine(scom);
            sw.WriteLine("iflag=," + iflag.ToString("0"));
            sw.WriteLine("iteration=," + iteration.ToString("0"));
            sw.WriteLine("MPARA=," + MPARA.ToString("0"));
            for (i = 0; i <= MPARA - 1; i++)
            {
                sw.WriteLine("parameter[" + i.ToString("0") + "]=," + parameter[i].ToString("E6"));
            }
            sw.Close();

            //画像出力
            PLTCTL(pltdata, NDATA, datax, datay, scom, parameter);
            pictureBox2.Image.Save(fnameP, System.Drawing.Imaging.ImageFormat.Png);

            //計算結果表示
            dat = scom + "\n";
            dat = dat + "iflag=" + iflag.ToString("0") + "\n";
            dat = dat + "iteration=" + iteration.ToString("0") + "\n";
            dat = dat + "MPARA=" + MPARA.ToString("0") + "\n";
            for (i = 0; i <= MPARA - 1; i++)
            {
                dat = dat + "parameter[" + i.ToString("0") + "]=" + parameter[i].ToString("E6") + "\n";
            }
            MessageBox.Show(dat, "計算結果表示");

        }
        //--------------------------------------------------------------------------------

        private int MPARA_SET()
        {
            //****************************************************
            //パラメータ数：MPARA（ここでは変数：mpara）のセット
            //****************************************************
            int mpara;

            mpara = 5;
            //mpara = 4;

            return mpara;
        }
        //--------------------------------------------------------------------------------
        private double Vfunc(double x, double[] parameter)
        {
            //**********************************************
            //関数値Vfuncの定義・計算
            //parameter[i]：回帰係数(i=0〜MPARA-1)
            //**********************************************

            //パソコンによるデータ解析入門p236(MPARA = 5)
            double gamma2, dx;
            gamma2 = parameter[1] * parameter[1];
            dx = x - parameter[0];
            return parameter[2] * gamma2 / (dx * dx + gamma2) + parameter[3] + parameter[4] * x;

            //減衰振動(MPARA=4)
            //return parameter[0] * Math.Exp(-parameter[1] * x) * Math.Sin(parameter[2] * x + parameter[3]);
        }

        //--------------------------------------------------------------------------------
        private void INIVAL(int MPARA, double[,] vertex, double[] stepsize, double[] tolerance)
        {
            //*******************************************
            //解析ﾊﾟﾗﾒ-ﾀ初期値等入力(0〜MPARA-1)
            //vertex：初期値
            //stepsize：初期の検索ステップ
            //tolerance：収束判定基準
            //*******************************************
            int j;
            for (j = 0; j <= MPARA - 1; j++)
            {
                vertex[0, j] = 1.0; stepsize[j] = 0.1 * vertex[0, j]; tolerance[j] = 0.000001;
            }
        }

        //--------------------------------------------------------------------------------
        private void SIMPLEX(int MPARA, double[] parameter, int NDATA, double[] datax, double[] datay,
                             out int iteration, out int iflag)
        {
            double[,] vertex = new double[MPARA + 1, MPARA];
            double[] criterion = new double[MPARA + 1];
            double[] stepsize = new double[MPARA];
            double[] tolerance = new double[MPARA];
            double s, v1, v2, cr, v, vmax, vmin, vsum, p, x, y;
            int i, j, k, kworst, kbest, ks;

            //*******************************************
            //解析ﾊﾟﾗﾒ-ﾀ初期値等入力
            //*******************************************
            INIVAL(MPARA, vertex, stepsize, tolerance);

            iflag = 0;
            ks = 0;
            s = Math.Sqrt(0.5);
            v1 = s * (Math.Sqrt((double)(MPARA + 1)) - 1) / (double)MPARA;
            for (j = 0; j <= MPARA - 1; j++)
            {
                v2 = vertex[0, j] + v1 * stepsize[j];
                for (k = 1; k <= MPARA; k++)
                {
                    vertex[k, j] = v2;
                }
                vertex[j + 1, j] = v2 + s * stepsize[j];
            }
            for (k = 0; k <= MPARA; k++)
            {
                for (j = 0; j <= MPARA - 1; j++)
                {
                    parameter[j] = vertex[k, j];
                }
                cr = 0;                 //Vcrit-Start
                for (i = 0; i <= NDATA; i++)
                {
                    x = datax[i];
                    y = Vfunc(x, parameter);
                    cr = cr + (y - datay[i]) * (y - datay[i]);
                }                     //Vcrit-End
                criterion[k] = cr;
            }
            iteration = 0;
            while (iteration <= 100000)
            {
                iteration = iteration + 1;
                kworst = 0; kbest = 0;
                for (k = 1; k <= MPARA; k++)
                {
                    if (criterion[k] > criterion[kworst]) kworst = k;
                    if (criterion[k] < criterion[kbest]) kbest = k;
                }
                iflag = 0;
                for (j = 0; j <= MPARA - 1; j++)
                {
                    vmax = vertex[0, j]; vmin = vmax;
                    for (k = 1; k <= MPARA; k++)
                    {
                        v = vertex[k, j];
                        if (v >= vmax) vmax = v;
                        if (v < vmin) vmin = v;
                    }
                    if (vmax - vmin > tolerance[j]) iflag = 1;
                }
                if (iflag == 0) break;

                for (j = 0; j <= MPARA - 1; j++)
                {
                    vsum = 0;
                    for (k = 0; k <= MPARA; k++)
                    {
                        if (k != kworst) vsum = vsum + vertex[k, j];
                    }
                    parameter[j] = 2 * vsum / MPARA - vertex[kworst, j];
                }
                cr = 0;                 //Vcrit-Start
                for (i = 0; i <= NDATA; i++)
                {
                    x = datax[i];
                    y = Vfunc(x, parameter);
                    cr = cr + (y - datay[i]) * (y - datay[i]);
                }                     //Vcrit-End
                if (cr < criterion[kbest]) ks = 0;
                if (criterion[kbest] <= cr && cr <= criterion[kworst]) ks = 1;
                if (criterion[kworst] < cr) ks = 2;
                switch (ks)
                {
                    case 0:
                        for (j = 0; j <= MPARA - 1; j++)
                        {
                            p = 1.5 * parameter[j] - 0.5 * vertex[kworst, j];
                            vertex[kworst, j] = parameter[j];
                            parameter[j] = p;
                        }
                        criterion[kworst] = cr;
                        cr = 0;                 //Vcrit-Start
                        for (i = 0; i <= NDATA; i++)
                        {
                            x = datax[i];
                            y = Vfunc(x, parameter);
                            cr = cr + (y - datay[i]) * (y - datay[i]);
                        }                     //Vcrit-End
                        if (cr < criterion[kworst])
                        {
                            for (j = 0; j <= MPARA - 1; j++)
                            {  //REPLACE-Start
                                vertex[kworst, j] = parameter[j];
                            }
                            criterion[kworst] = cr; //REPLACE-End
                        }
                        break;
                    case 1:
                        for (j = 0; j <= MPARA - 1; j++)
                        {  //REPLACE-Start
                            vertex[kworst, j] = parameter[j];
                        }
                        criterion[kworst] = cr; //REPLACE-End
                        break;
                    case 2:
                        for (j = 0; j <= MPARA - 1; j++)
                        {
                            parameter[j] = 0.75 * vertex[kworst, j] + 0.25 * parameter[j];
                        }
                        cr = 0;                 //Vcrit-Start
                        for (i = 0; i <= NDATA; i++)
                        {
                            x = datax[i];
                            y = Vfunc(x, parameter);
                            cr = cr + (y - datay[i]) * (y - datay[i]);
                        }                     //Vcrit-End
                        if (cr < criterion[kworst])
                        {
                            for (j = 0; j <= MPARA - 1; j++)
                            {  //REPLACE-Start
                                vertex[kworst, j] = parameter[j];
                            }
                            criterion[kworst] = cr; //REPLACE-End
                        }
                        else
                        {
                            for (k = 0; k <= MPARA; k++)
                            {
                                if (k != kbest)
                                {
                                    for (j = 0; j <= MPARA - 1; j++)
                                    {
                                        parameter[j] = 0.5 * (vertex[k, j] + vertex[kbest, j]);
                                        vertex[k, j] = parameter[j];
                                    }
                                    cr = 0;                 //Vcrit-Start
                                    for (i = 0; i <= NDATA; i++)
                                    {
                                        x = datax[i];
                                        y = Vfunc(x, parameter);
                                        cr = cr + (y - datay[i]) * (y - datay[i]);
                                    }                     //Vcrit-End
                                    criterion[k] = cr;
                                }
                            }
                        }
                        break;
                }
            }

        }
        //---------------------------------------------------------------------------
        private void PLTCTL(DRAWPARA pltdata, int nd, double[] datax, double[] datay, string scom, double[] parameter)
        {
            //nd：データ数（入力はデータ数，プログラム上数値はデータ数−１）
            //datax[]：プロットデータｘ座標
            //datay[]：プロットデータｙ座標
            Bitmap bmp;
            Graphics g;
            PictureBox thePicBox = pictureBox1;
            int i;

            //プロット領域定義
            pltdata.kxi0 = pltdata.DHL;
            pltdata.kxf0 = pltdata.HSIZE - pltdata.DHR;
            pltdata.kyi0 = pltdata.DVU;
            pltdata.kyf0 = pltdata.VSIZE - pltdata.DVL;

            for (i = 0; i <= 1; i++)
            {
                switch (i)
                {
                    case 0: //画面描画
                        pltdata.kpt = 1;
                        thePicBox = pictureBox1;
                        thePicBox.Visible = true;
                        break;
                    case 1: //画像書き出し
                        pltdata.kpt = 2;
                        thePicBox = pictureBox2;
                        thePicBox.Visible = false;
                        break;
                }
                thePicBox.Size = new Size(pltdata.kpt * pltdata.HSIZE, pltdata.kpt * pltdata.VSIZE);
                bmp = new Bitmap(thePicBox.Width, thePicBox.Height);
                thePicBox.Image = bmp;
                g = Graphics.FromImage(thePicBox.Image);
                g.FillRectangle(Brushes.White, 0, 0, thePicBox.Width, thePicBox.Height);
                PLOT(g, pltdata, nd, datax, datay, scom, parameter);
                g.Dispose();
            }
            pictureBox1.Left = 0;
            pictureBox1.Top = toolStrip1.Height;
            ClientSize = new Size(pictureBox1.Width, pictureBox1.Height + toolStrip1.Height);
        }
        //---------------------------------------------------------------------------------
        private void PLOT(Graphics g, DRAWPARA pltdata, int nd, double[] datax, double[] datay, string scom, double[] parameter)
        {
            int i, nn;
            int kxi, kxf, kyi, kyf;
            double xx, yy;
            int kxx, kyy;
            int kxx1, kyy1, kxx2, kyy2;
            string str;
            System.Drawing.Font f = new Font("ＭＳゴシック", pltdata.kpt * 10);
            System.Drawing.SizeF TextSize1;
            System.Drawing.SizeF TextSize2;

            kxi = pltdata.kpt * pltdata.kxi0;
            kxf = pltdata.kpt * pltdata.kxf0;
            kyi = pltdata.kpt * pltdata.kyi0;
            kyf = pltdata.kpt * pltdata.kyf0;
            //座標軸
            System.Drawing.Pen LPen = new Pen(System.Drawing.Color.Black);
            LPen.DashStyle = System.Drawing.Drawing2D.DashStyle.Dot;
            DRN_XJIKU(g, f, LPen, pltdata, kxi, kxf, kyi, kyf);
            DRN_YJIKU(g, f, LPen, pltdata, kxi, kxf, kyi, kyf);
            //枠線描画
            g.DrawRectangle(new Pen(Color.Black, 2), kxi, kyi, kxf - kxi, kyf - kyi);
            //ｙ軸名描画
            str = pltdata.syjiku;
            TextSize2 = g.MeasureString(str, f);
            TextSize1 = g.MeasureString("-1.0", f);
            kxx = (int)(kxi - pltdata.kpt * 10 - TextSize1.Width - TextSize2.Height);
            kyy = (int)((kyi + kyf) / 2 + TextSize2.Width / 2);
            INC_STR(g, f, str, kxx, kyy, -90);
            //ｘ軸名描画
            str = pltdata.sxjiku;
            TextSize2 = g.MeasureString(str, f);
            kxx = (int)((kxi + kxf) / 2 - TextSize2.Width / 2);
            kyy = (int)(kyf + pltdata.kpt * 5 + TextSize1.Height);
            INC_STR(g, f, str, kxx, kyy, 0);
            //データプロット
            int dc;
            dc = pltdata.kpt * 5;
            for (i = 0; i <= nd; i++)
            {
                xx = datax[i];
                yy = datay[i];
                kxx = kxi + (int)((xx - pltdata.xmin) * (kxf - kxi) / (pltdata.xmax - pltdata.xmin)) - (int)(dc / 2);
                kyy = kyf - (int)((yy - pltdata.ymin) * (kyf - kyi) / (pltdata.ymax - pltdata.ymin)) - (int)(dc / 2);
                g.DrawEllipse(new Pen(Color.Red, 1), kxx, kyy, dc, dc);
            }
            //回帰推定線
            nn = 500;
            xx = pltdata.xmin;
            yy = Vfunc(xx, parameter);
            kxx1 = kxi + (int)((xx - pltdata.xmin) * (kxf - kxi) / (pltdata.xmax - pltdata.xmin));
            kyy1 = kyf - (int)((yy - pltdata.ymin) * (kyf - kyi) / (pltdata.ymax - pltdata.ymin));
            for (i = 1; i <= nn; i++)
            {
                xx = xx + (pltdata.xmax - pltdata.xmin) / (double)nn;
                yy = Vfunc(xx, parameter);
                kxx2 = kxi + (int)((xx - pltdata.xmin) * (kxf - kxi) / (pltdata.xmax - pltdata.xmin));
                kyy2 = kyf - (int)((yy - pltdata.ymin) * (kyf - kyi) / (pltdata.ymax - pltdata.ymin));
                g.DrawLine(new Pen(Color.Blue, 2), kxx1, kyy1, kxx2, kyy2);
                kxx1 = kxx2;
                kyy1 = kyy2;
            }
            //右上コメント
            str = scom;
            TextSize1 = g.MeasureString(str, f);
            kxx = kxf - 2 * pltdata.kpt - (int)TextSize1.Width;
            kyy = kyi - 2 * pltdata.kpt - (int)TextSize1.Height;
            g.DrawString(str, f, Brushes.Black, kxx, kyy);
            f.Dispose();
        }
        //---------------------------------------------------------------------------------
        private void INC_STR(Graphics g, System.Drawing.Font f, string str, int kxx, int kyy, float ang)
        {
            //軸名描画
            g.ScaleTransform(1, 1);//横・縦の表示比率を設定
            g.TranslateTransform(kxx, kyy);//表示位置の設定（表示位置を原点とする座標移動）
            g.RotateTransform(ang);//表示角度を指定
            g.DrawString(str, f, Brushes.Black, 0, 0);//描画実行
            g.ResetTransform();//単位行列にリセット
        }
        //---------------------------------------------------------------------------------
        private void DRN_XJIKU(Graphics g, System.Drawing.Font f, System.Drawing.Pen LPen,
                               DRAWPARA pltdata, int kxi, int kxf, int kyi, int kyf)
        {
            int i, ix, kxx;
            double xx, wv;
            string str;
            System.Drawing.SizeF TextSize1;
            //普通ｘ軸描画
            ix = (int)((pltdata.xmax - pltdata.xmin) / pltdata.dx);
            for (i = 0; i <= ix; i++)
            {
                xx = pltdata.xmin + (double)(i) * pltdata.dx;
                kxx = kxi + (int)((xx - pltdata.xmin) * (kxf - kxi) / (pltdata.xmax - pltdata.xmin));
                wv = pltdata.xmin + (double)(i) * pltdata.dx;
                str = wv.ToString("0");
                if ((int)(pltdata.dx * 1000.0) % 10 != 0) { str = wv.ToString("0.000"); }
                if ((int)(pltdata.dx * 1000.0) % 10 == 0 && (int)(pltdata.dx * 100.0) % 10 != 0) { str = wv.ToString("0.00"); }
                if ((int)(pltdata.dx * 1000.0) % 10 == 0 && (int)(pltdata.dx * 100.0) % 10 == 0 && (int)(pltdata.dx * 10.0) % 10 != 0) { str = wv.ToString("0.0"); }
                TextSize1 = g.MeasureString(str, f);
                g.DrawLine(LPen, kxx, kyi, kxx, kyf);
                g.DrawString(str, f, Brushes.Black, kxx - TextSize1.Width / 2, kyf + pltdata.kpt * 3);
            }
        }
        //---------------------------------------------------------------------------------
        private void DRN_YJIKU(Graphics g, System.Drawing.Font f, System.Drawing.Pen LPen,
                               DRAWPARA pltdata, int kxi, int kxf, int kyi, int kyf)
        {
            int i, iy, kyy;
            double yy, wv;
            string str;
            System.Drawing.SizeF TextSize1;
            //普通ｙ軸描画
            iy = (int)((pltdata.ymax - pltdata.ymin) / pltdata.dy);
            for (i = 0; i <= iy; i++)
            {
                yy = pltdata.ymin + (double)(i) * pltdata.dy;
                kyy = kyf - (int)((yy - pltdata.ymin) * (kyf - kyi) / (pltdata.ymax - pltdata.ymin));
                wv = pltdata.ymin + (int)(i) * pltdata.dy;
                str = wv.ToString("0");
                if ((int)(pltdata.dy * 1000.0) % 10 != 0) { str = wv.ToString("0.000"); }
                if ((int)(pltdata.dy * 1000.0) % 10 == 0 && (int)(pltdata.dy * 100.0) % 10 != 0) { str = wv.ToString("0.00"); }
                if ((int)(pltdata.dy * 1000.0) % 10 == 0 && (int)(pltdata.dy * 100.0) % 10 == 0 && (int)(pltdata.dy * 10.0) % 10 != 0) { str = wv.ToString("0.0"); }
                TextSize1 = g.MeasureString(str, f);
                g.DrawLine(LPen, kxi, kyy, kxf, kyy);
                g.DrawString(str, f, Brushes.Black, kxi - TextSize1.Width - pltdata.kpt * 2, kyy - TextSize1.Height / 2);
            }
        }
        //---------------------------------------------------------------------------------

        private void INI_PLTDATA(out DRAWPARA pltdata)
        {
            //pltdata初期化
            pltdata.kpt = 0;
            pltdata.HSIZE = 0;
            pltdata.DHL = 0;
            pltdata.DHR = 0;
            pltdata.VSIZE = 0;
            pltdata.DVL = 0;
            pltdata.DVU = 0;
            pltdata.sxjiku = "";
            pltdata.xmin = 0.0;
            pltdata.xmax = 0.0;
            pltdata.dx = 0.0;
            pltdata.syjiku = "";
            pltdata.ymin = 0.0;
            pltdata.ymax = 0.0;
            pltdata.dy = 0.0;
            pltdata.kxi0 = 0;
            pltdata.kxf0 = 0;
            pltdata.kyi0 = 0;
            pltdata.kyf0 = 0;
        }
        //---------------------------------------------------------------------------------

    }
}