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 vcsMDA
{
    public partial class Form1 : Form
    {
        public Form1()
        {
            InitializeComponent();
        }
        //--------------------------------------------------------------------

        private void toolStripButton1_Click(object sender, EventArgs e)
        {
            //重判別分析（multiple discriminant analysis）
            int i, j, k, n;
            int mm, ndata, nclass, ndf, iflg;
            double tr, aa;
            string scom;

            System.IO.StreamReader sr;
            System.IO.StreamWriter sw;
            string fnameR = "";
            string fnameW = "";
            string dat;
            string[] sbuf;
            char delim = ',';

            //入出力ファイル指定
            openFileDialog1.InitialDirectory = System.IO.Directory.GetCurrentDirectory();
            if (openFileDialog1.ShowDialog() == DialogResult.OK) { fnameR = openFileDialog1.FileName; }
            if (saveFileDialog1.ShowDialog() == DialogResult.OK) { fnameW = saveFileDialog1.FileName; }

            //データ読み込み
            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); mm = int.Parse(sbuf[0]); ndata = int.Parse(sbuf[1]);

            double[,] xd = new double[ndata, mm];
            int[] iclass = new int[ndata];
            int[] members = new int[ndata];
            double[,] classmean = new double[ndata, mm];
            double[] xm = new double[mm];
            double[] temp = new double[mm];
            double[,] a = new double[mm, mm];
            double[] ev = new double[mm];
            double[,] vec = new double[mm, mm];

            k = 0;
            while (!sr.EndOfStream)
            {
                dat = sr.ReadLine(); sbuf = dat.Split(delim);
                iclass[k] = int.Parse(sbuf[0]);
                for (i = 0; i <= mm - 1; i++)
                {
                    xd[k, i] = double.Parse(sbuf[i + 1]);
                }
                k = k + 1;
            }
            sr.Close();
            ndata = k - 1;

            //クラス毎データ組数（クラスは0から小さい順に並んでいる前提）
            nclass = 0;
            for (k = 0; k <= ndata; k++)
            {
                if (iclass[k] >= nclass) nclass = iclass[k];
            }
            //平均値の計算
            for (i = 0; i <= mm - 1; i++)
            {
                xm[i] = 0.0;
                for (k = 0; k <= ndata; k++)
                {
                    xm[i] = xm[i] + xd[k, i];
                }
                xm[i] = xm[i] / (double)(ndata + 1);
            }
            //平均値を引く
            for (i = 0; i <= mm - 1; i++)
            {
                for (k = 0; k <= ndata; k++)
                {
                    xd[k, i] = xd[k, i] - xm[i];
                }
            }
            //クラス平均
            for (k = 0; k <= ndata; k++)
            {
                j = iclass[k];
                members[j] = members[j] + 1;
                for (i = 0; i <= mm - 1; i++)
                {
                    classmean[j, i] = classmean[j, i] + xd[k, i];
                }
            }
            for (i = 0; i <= nclass; i++)
            {
                for (j = 0; j <= mm - 1; j++)
                {
                    classmean[i, j] = classmean[i, j] / (double)members[i];
                }
            }
            //固有方程式
            for (i = 0; i <= mm - 1; i++)
            {
                for (j = i; j <= mm - 1; j++)
                {
                    aa = 0.0;
                    for (k = 0; k <= ndata; k++)
                    {
                        aa = aa + (xd[k, i] - classmean[iclass[k], i]) * (xd[k, j] - classmean[iclass[k], j]);
                    }
                    a[i, j] = aa; a[j, i] = a[i, j];
                }
            }
            //**********************************************
            iflg = 0; MJACOBI(mm, a, ev, vec, iflg);
            //**********************************************
            ndf = ndata - nclass;
            for (i = 0; i <= mm - 1; i++)
            {
                aa = Math.Sqrt(ev[i] / (double)ndf);
                for (j = 0; j <= mm - 1; j++)
                {
                    vec[j, i] = vec[j, i] / aa;
                    a[i, j] = 0.0;
                }
            }
            for (k = 0; k <= nclass; k++)
            {
                for (i = 0; i <= mm - 1; i++)
                {
                    temp[i] = 0.0;
                    for (j = 0; j <= mm - 1; j++)
                    {
                        temp[i] = temp[i] + vec[j, i] * classmean[k, j];
                    }
                }
                for (i = 0; i <= mm - 1; i++)
                {
                    for (j = i; j <= mm - 1; j++)
                    {
                        a[i, j] = a[i, j] + (double)members[k] * temp[i] * temp[j];
                        a[j, i] = a[i, j];
                    }
                }
            }
            //**********************************************
            iflg = 1; MJACOBI(mm, a, ev, vec, iflg);
            //**********************************************
            tr = 0.0;
            for (i = 0; i <= mm - 1; i++)
            {
                tr = tr + ev[i];
            }

            double[,] z = new double[ndata + 1, mm];
            double[,] d = new double[ndata + 1, nclass + 1];
            //クラス平均座標表示
            for (n = 0; n <= nclass; n++)
            {
                for (i = 0; i <= mm - 1; i++)
                {
                    z[n, i] = 0.0;
                    for (j = 0; j <= mm - 1; j++)
                    {
                        z[n, i] = z[n, i] + vec[j, i] * classmean[n, j];
                    }
                }
            }
            for (k = 0; k <= ndata; k++)
            {
                //データ毎座標表示
                for (i = 0; i <= mm - 1; i++)
                {
                    temp[i] = xd[k, i];
                }
                for (i = 0; i <= mm - 1; i++)
                {
                    xd[k, i] = 0.0;
                    for (j = 0; j <= mm - 1; j++)
                    {
                        xd[k, i] = xd[k, i] + vec[j, i] * temp[j];
                    }
                }
                //クラス平均からの距離計算
                for (n = 0; n <= nclass; n++)
                {
                    d[k, n] = 0.0;
                    for (j = 0; j <= mm - 1; j++)
                    {
                        d[k, n] = d[k, n] + (xd[k, j] - z[n, j]) * (xd[k, j] - z[n, j]);
                    }
                    d[k, n] = Math.Sqrt(d[k, n]);
                }
            }

            //データ出力
            sw = new System.IO.StreamWriter(fnameW, false, System.Text.Encoding.Default);
            dat = scom; sw.WriteLine(dat);
            dat = "No,"; for (i = 0; i <= mm - 1; i++) { dat = dat + "," + (i + 1).ToString("0"); }; sw.WriteLine(dat);
            dat = "固有値,"; for (i = 0; i <= mm - 1; i++) { dat = dat + "," + ev[i].ToString("0.000"); }; sw.WriteLine(dat);
            dat = "寄与率,"; for (i = 0; i <= mm - 1; i++) { dat = dat + "," + (ev[i] / tr).ToString("0.000"); }; sw.WriteLine(dat);
            dat = "クラス平均正準座標"; sw.WriteLine(dat);
            for (n = 0; n <= nclass; n++)
            {
                dat = "size=" + members[n].ToString("0") + ",class" + n.ToString("0");
                for (i = 0; i <= mm - 1; i++)
                {
                    dat = dat + "," + z[n, i].ToString("0.000");
                }
                sw.WriteLine(dat);
            }
            dat = "データ毎正準座標,";
            for (i = 0; i <= mm - 1; i++) { dat = dat + "," + (i + 1).ToString("0"); }
            for (n = 0; n <= nclass; n++) { dat = dat + ",dis-c" + n.ToString("0"); }
            sw.WriteLine(dat);
            for (k = 0; k <= ndata; k++)
            {
                dat = "data" + (k + 1).ToString("00") + ",class" + iclass[k].ToString("0");
                for (i = 0; i <= mm - 1; i++)
                {
                    dat = dat + "," + xd[k, i].ToString("0.000");
                }
                for (n = 0; n <= nclass; n++)
                {
                    dat = dat + "," + d[k, n].ToString("0.000");
                }
                sw.WriteLine(dat);
            }
            sw.Close();

            //新規データ評価実施判断
            string newcase = "";
            switch (MessageBox.Show("新規データ評価実施", "ファイル読込", MessageBoxButtons.YesNo, MessageBoxIcon.Question))
            {
                case DialogResult.Yes: newcase = "yes"; break;
                case DialogResult.No: newcase = "no"; break;
            }
            if (newcase == "yes")
            {
                int kdata;
                //新規評価データ読み込み
                openFileDialog1.InitialDirectory = System.IO.Directory.GetCurrentDirectory();
                if (openFileDialog1.ShowDialog() == DialogResult.OK) { fnameR = openFileDialog1.FileName; }
                sr = new System.IO.StreamReader(fnameR, System.Text.Encoding.Default);
                dat = sr.ReadLine();
                dat = sr.ReadLine(); sbuf = dat.Split(delim); mm = int.Parse(sbuf[0]); kdata = int.Parse(sbuf[1]);
                double[,] xnew = new double[kdata, mm];
                k = 0;
                while (!sr.EndOfStream)
                {
                    dat = sr.ReadLine(); sbuf = dat.Split(delim);
                    for (i = 0; i <= mm - 1; i++)
                    {
                        xnew[k, i] = double.Parse(sbuf[i]);
                    }
                    k = k + 1;
                }
                sr.Close();
                kdata = k - 1;

                for (k = 0; k <= kdata; k++)
                {
                    //座標計算
                    for (i = 0; i <= mm - 1; i++)
                    {
                        z[k, i] = 0.0;
                        for (j = 0; j <= mm - 1; j++)
                        {
                            z[k, i] = z[k, i] + vec[j, i] * (xnew[k, j] - xm[j]);
                        }
                    }
                    //クラス平均からの距離計算
                    for (n = 0; n <= nclass; n++)
                    {
                        d[k, n] = 0.0;
                        for (i = 0; i <= mm - 1; i++)
                        {
                            aa = 0.0;
                            for (j = 0; j <= mm - 1; j++)
                            {
                                aa = aa + vec[j, i] * (xnew[k, j] - xm[j] - classmean[n, j]);
                            }
                            d[k, n] = d[k, n] + aa * aa;
                        }
                        d[k, n] = Math.Sqrt(d[k, n]);
                    }
                }

                //データ追加出力
                sw = new System.IO.StreamWriter(fnameW, true, System.Text.Encoding.Default);
                if (newcase == "yes")
                {
                    //評価データ諸元書き込み
                    dat = "新規評価データ正準座標"; sw.WriteLine(dat);
                    for (k = 0; k <= kdata; k++)
                    {
                        dat = "new," + k.ToString("00");
                        for (i = 0; i <= mm - 1; i++)
                        {
                            dat = dat + "," + z[k, i].ToString("0.000");
                        }
                        for (n = 0; n <= nclass; n++)
                        {
                            dat = dat + "," + d[k, n].ToString("0.000");
                        }
                        sw.WriteLine(dat);
                    }
                }
                sw.Close();
            }

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

        private void MJACOBI(int n, double[,] a, double[] ev, double[,] vec, int iflg)
        {
            //Jacobi法による固有値・固有ベクトル計算（奥村先生）
            //固有値ev[k],固有ベクトルvec[i,k]
            int i, j, k, im = 0;
            double s, c, aoff, eps, t, aa = 0.0, aj, ak, vj, vk, vmax, bb;

            if (iflg == 0)
            { //固有ベクトル初期化
                for (i = 0; i <= n - 1; i++)
                {
                    for (j = 0; j <= n - 1; j++)
                    {
                        vec[i, j] = 0.0;
                    }
                    vec[i, i] = 1.0;
                }
            }

            aoff = 0.0;
            for (j = 0; j <= n - 2; j++)
            {
                for (k = j + 1; k <= n - 1; k++)
                {
                    aoff = aoff + a[j, k] * a[j, k];
                }
            }
            eps = aoff * 2.0;
            for (j = 0; j <= n - 1; j++)
            {
                eps = eps + a[j, j] * a[j, j];
            }
            eps = eps * 0.00000000005;

            while (aoff > eps)
            {
                for (j = 0; j <= n - 2; j++)
                {
                    for (k = j + 1; k <= n - 1; k++)
                    {
                        if (a[j, k] == 0.0)
                        {
                            t = 0.0;
                        }
                        else
                        {
                            aa = (a[k, k] - a[j, j]) / (2.0 * a[j, k]);
                        }
                        if (aa >= 0)
                        {
                            t = 1.0 / (aa + Math.Sqrt(1.0 + aa * aa));
                        }
                        else
                        {
                            t = 1.0 / (aa - Math.Sqrt(1.0 + aa * aa));
                        }
                        c = 1.0 / Math.Sqrt(1.0 + t * t);
                        s = t * c;
                        a[j, j] = a[j, j] - t * a[j, k];
                        a[k, k] = a[k, k] + t * a[j, k];
                        a[j, k] = 0.0;
                        for (i = 0; i <= j - 1; i++)
                        {
                            aj = a[i, j];
                            ak = a[i, k];
                            a[i, j] = aj * c - ak * s;
                            a[i, k] = aj * s + ak * c;
                        }
                        for (i = j + 1; i <= k - 1; i++)
                        {
                            aj = a[j, i];
                            ak = a[i, k];
                            a[j, i] = aj * c - ak * s;
                            a[i, k] = aj * s + ak * c;
                        }
                        for (i = k + 1; i <= n - 1; i++)
                        {
                            aj = a[j, i];
                            ak = a[k, i];
                            a[j, i] = aj * c - ak * s;
                            a[k, i] = aj * s + ak * c;
                        }
                        for (i = 0; i <= n - 1; i++)
                        {
                            vj = vec[i, j];
                            vk = vec[i, k];
                            vec[i, j] = vj * c - vk * s;
                            vec[i, k] = vj * s + vk * c;
                        }
                    }
                }
                aoff = 0.0;
                for (j = 0; j <= n - 2; j++)
                {
                    for (k = j + 1; k <= n - 1; k++)
                    {
                        aoff = aoff + a[j, k] * a[j, k];
                    }
                }
            }

            for (k = 0; k <= n - 1; k++)
            {
                ev[k] = a[k, k];
            }
            for (k = 0; k <= n - 2; k++)
            {
                vmax = -1.0E+30;
                for (i = k; i <= n - 1; i++)
                {
                    if (ev[i] >= vmax)
                    {
                        im = i;
                        vmax = ev[i];
                    }
                }
                aa = ev[k]; ev[k] = ev[im]; ev[im] = aa;
                for (i = 0; i <= n - 1; i++)
                {
                    bb = vec[i, k];
                    vec[i, k] = vec[i, im];
                    vec[i, im] = bb;
                }
            }

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