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

        private void Form1_Load(object sender, EventArgs e)
        {
            listView1.Clear();
            listView1.View = View.Details;
            listView1.Columns.Add("No[i]", 100, HorizontalAlignment.Right);
            listView1.Columns.Add("1", 70, HorizontalAlignment.Right);
            listView1.Columns.Add("2", 70, HorizontalAlignment.Right);
            listView1.Columns.Add("3", 70, HorizontalAlignment.Right);
            listView1.Columns.Add("4", 70, HorizontalAlignment.Right);
        }

        //-----------------------------------------------------------------------
        private void toolStripButton1_Click(object sender, EventArgs e)
        {
	        //主成分分析（principal component analysis）
    	    System.IO.StreamReader sr;
        	System.IO.StreamWriter sw;
        	string dat="";
        	string[] sbuf;
        	char delim=',';
        	string fnameR = "";
        	string fnameW = "";
        
			int i,j,k;
        
			int    mm;    //変数の数−１
        	string scom;  //コメント
        	int    ndata; //データ組数−１
        	double tr;    //固有値総和（分散共分散行列対角項総和）

            //入出力ファイル指定
            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]) - 1;
        
			double[,] xd  =new double[ndata+1, mm]; //データ入力用配列
        	double[,] a   =new double[mm, mm];      //分散共分散行列格納配列
        	double[]  xm  =new double[mm];          //平均値格納配列
        	double[]  sd  =new double[mm];          //標準偏差格納配列
        	double[]  ev  =new double[mm];          //固有値
        	double[,] vec =new double[mm, mm];      //固有ベクトル
        	double[,] zd  =new double[ndata+1, mm]; //主成分得点（標準化変数×固有ベクトル）
        	
			k = 0;
            while(!sr.EndOfStream){
                dat = sr.ReadLine(); sbuf = dat.Split(delim);
            	for( i = 0;i<=mm - 1;i++){
                	xd[k, i] = double.Parse(sbuf[i]);
            	}
            	k = k + 1;
        	}
        	sr.Close();
        	ndata = k - 1;

            //平均値の計算
            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++){
                sd[i] = 0.0;
                for( k = 0;k<=ndata;k++){
                    sd[i] = sd[i] + (xd[k, i] - xm[i]) * (xd[k, i] - xm[i]);
                }
                sd[i] = Math.Sqrt(sd[i] / (double)(ndata));
            }
            //変数の標準化
            for( i = 0;i<=mm - 1;i++){
                for( k = 0;k<=ndata;k++){
                    xd[k, i] = (xd[k, i] - xm[i]) / sd[i];
                }
            }
            //分散共分散行列(相関行列)の計算
            tr = 0.0;
            for(i = 0;i<=mm - 1;i++){
                for(j = 0;j<=mm - 1;j++){
                    a[i, j] = 0.0;
                    for(k = 0;k<=ndata;k++){
                        a[i, j] = a[i, j] + xd[k, i] * xd[k, j];
                    }
                    a[i, j] = a[i, j] / (double)ndata;
                }
                tr = tr + a[i, i]; //対角項の総和（=固有値の総和）
            }
            //固有値・固有ベクトルの計算
            JACOBI(mm, a, ev, vec);
            //主成分得点の計算
            for(k = 0;k<=ndata;k++){
                for(i = 0;i<=mm - 1;i++){
                    zd[k, i] = 0.0;
                    for(j = 0;j<=mm - 1;j++){
                        zd[k, i] = zd[k, i] + xd[k, j] * vec[j, i];
                    }
                }
            }

            //リストビューへの書き出し
            ListViewItem lvitem;
            lvitem = new ListViewItem("固有値");
            lvitem.SubItems.Add(ev[0].ToString("0.000e+00"));
            lvitem.SubItems.Add(ev[1].ToString("0.000e+00"));
            if(3 <= mm)lvitem.SubItems.Add(ev[2].ToString("0.000e+00"));
            if(4 <= mm)lvitem.SubItems.Add(ev[3].ToString("0.000e+00"));
            listView1.Items.Add(lvitem);
            lvitem = new ListViewItem("寄与率");
            lvitem.SubItems.Add((ev[0] / tr).ToString("0.000"));
            lvitem.SubItems.Add((ev[1] / tr).ToString("0.000"));
            if(3 <= mm)lvitem.SubItems.Add((ev[2] / tr).ToString("0.000"));
            if(4 <= mm)lvitem.SubItems.Add((ev[3] / tr).ToString("0.000"));
            listView1.Items.Add(lvitem);
            for(j = 0;j<=mm - 1;j++){
                lvitem = new ListViewItem("固有ベクトル[i," + (j + 1).ToString() + "]");
                lvitem.SubItems.Add(vec[j, 0].ToString("0.000e+00"));
                lvitem.SubItems.Add(vec[j, 1].ToString("0.000e+00"));
                if(3 <= mm)lvitem.SubItems.Add(vec[j, 2].ToString("0.000e+00"));
                if(4 <= mm)lvitem.SubItems.Add(vec[j, 3].ToString("0.000e+00"));
                listView1.Items.Add(lvitem);
            }

            //ファイル出力
            sw = new System.IO.StreamWriter(fnameW, false, System.Text.Encoding.Default);
            sw.WriteLine(scom);
            dat = "No."    ; for (i = 0; i <= mm - 1; i++) { dat = dat + "," + (i + 1).ToString("0"); }            ;sw.WriteLine(dat); //左から固有値の大きい順
            dat = "e-value"; for (i = 0; i <= mm - 1; i++) { dat = dat + "," + ev[i].ToString("0.000000"); }       ;sw.WriteLine(dat);  //固有値
            dat = "Prop"   ; for (i = 0; i <= mm - 1; i++) { dat = dat + "," + (ev[i] / tr).ToString("0.000000"); };sw.WriteLine(dat); //寄与率
            for(j=0;j<=mm-1;j++){
                dat = "e-vector" + (j + 1).ToString("0"); for (i = 0; i <= mm - 1; i++) { dat = dat + "," + vec[j, i].ToString("0.000000"); };sw.WriteLine(dat); //固有ベクトル
            }
            for(k=0;k<=ndata;k++){
                dat = "data" + (k + 1).ToString("0") ;
                for (i = 0; i <= mm - 1; i++) { dat = dat + "," + zd[k, i].ToString("0.000000"); }; sw.WriteLine(dat); //各データの主成分
            }
            sw.Close();

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

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

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

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

    }
}