Option Explicit On
Option Strict On
Public Class Form1

    Private Sub ToolStripButton1_Click(ByVal sender As System.Object, ByVal e As System.EventArgs) Handles ToolStripButton1.Click
        '主成分分析（principal component analysis）
        Dim sr As System.IO.StreamReader
        Dim sw As System.IO.StreamWriter
        Dim dat As String
        Dim sbuf() As String
        Dim delim() As Char = {","c}
        Dim fname1 As String = ""
        Dim fname2 As String = ""
        Dim i As Integer
        Dim j As Integer
        Dim k As Integer
        Dim mm As Integer     '変数の数−１
        Dim scom As String    'コメント
        Dim ndata As Integer  'データ組数−１
        Dim xd(,) As Double   'データ入力用配列
        Dim xm() As Double    '平均値格納配列
        Dim sd() As Double    '標準偏差格納配列
        Dim a(,) As Double    '分散共分散行列格納配列
        Dim ev() As Double    '固有値
        Dim vec(,) As Double  '固有ベクトル
        Dim zd(,) As Double   '主成分得点（標準化変数×固有ベクトル）
        Dim tr As Double      '固有値総和（分散共分散行列対角項総和）

        If OpenFileDialog1.ShowDialog() = Windows.Forms.DialogResult.OK Then fname1 = OpenFileDialog1.FileName
        sr = New System.IO.StreamReader(fname1, System.Text.Encoding.Default)
        dat = sr.ReadLine() : sbuf = dat.Split(delim) : scom = sbuf(0)
        dat = sr.ReadLine() : sbuf = dat.Split(delim) : mm = CInt(sbuf(0)) : ndata = CInt(sbuf(1)) - 1
        ReDim xd(ndata, mm - 1)
        k = 0
        Do While sr.Peek() >= 0
            dat = sr.ReadLine() : sbuf = dat.Split(delim)
            For i = 0 To mm - 1
                xd(k, i) = CDbl(sbuf(i))
            Next i
            k = k + 1
        Loop
        sr.Close()
        ndata = k - 1

        ReDim a(mm - 1, mm - 1)
        ReDim xm(mm - 1)
        ReDim sd(mm - 1)
        ReDim ev(mm - 1)
        ReDim vec(mm - 1, mm - 1)
        ReDim zd(ndata, mm - 1)

        '平均値の計算
        For i = 0 To mm - 1
            xm(i) = 0.0
            For k = 0 To ndata
                xm(i) = xm(i) + xd(k, i)
            Next k
            xm(i) = xm(i) / CDbl(ndata + 1)
        Next i
        '標準偏差の計算
        For i = 0 To mm - 1
            sd(i) = 0.0
            For k = 0 To ndata
                sd(i) = sd(i) + (xd(k, i) - xm(i)) * (xd(k, i) - xm(i))
            Next k
            sd(i) = Math.Sqrt(sd(i) / CDbl(ndata))
        Next i
        '変数の標準化
        For i = 0 To mm - 1
            For k = 0 To ndata
                xd(k, i) = (xd(k, i) - xm(i)) / sd(i)
            Next k
        Next i
        '分散共分散行列(相関行列)の計算
        tr = 0.0
        For i = 0 To mm - 1
            For j = 0 To mm - 1
                a(i, j) = 0.0
                For k = 0 To ndata
                    a(i, j) = a(i, j) + xd(k, i) * xd(k, j)
                Next k
                a(i, j) = a(i, j) / CDbl(ndata)
            Next j
            tr = tr + a(i, i) '対角項の総和（=固有値の総和）
        Next i
        '固有値・固有ベクトルの計算
        Call JACOBI(mm, a, ev, vec)
        '主成分得点の計算
        For k = 0 To ndata
            For i = 0 To mm - 1
                zd(k, i) = 0.0
                For j = 0 To mm - 1
                    zd(k, i) = zd(k, i) + xd(k, j) * vec(j, i)
                Next j
            Next i
        Next k

        'リストビューへの書き出し
        Dim lvitem As ListViewItem
        lvitem = New ListViewItem("固有値")
        lvitem.SubItems.Add(ev(0).ToString("0.000e+00"))
        lvitem.SubItems.Add(ev(1).ToString("0.000e+00"))
        If 3 <= mm Then lvitem.SubItems.Add(ev(2).ToString("0.000e+00"))
        If 4 <= mm Then 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 Then lvitem.SubItems.Add((ev(2) / tr).ToString("0.000"))
        If 4 <= mm Then lvitem.SubItems.Add((ev(3) / tr).ToString("0.000"))
        ListView1.Items.Add(lvitem)
        For j = 0 To mm - 1
            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 Then lvitem.SubItems.Add(vec(j, 2).ToString("0.000e+00"))
            If 4 <= mm Then lvitem.SubItems.Add(vec(j, 3).ToString("0.000e+00"))
            ListView1.Items.Add(lvitem)
        Next j

        If SaveFileDialog1.ShowDialog() = Windows.Forms.DialogResult.OK Then fname2 = SaveFileDialog1.FileName
        sw = New System.IO.StreamWriter(fname2, False, System.Text.Encoding.Default)
        sw.WriteLine(scom)
        dat = "No." : For i = 0 To mm - 1 : dat = dat & "," & (i + 1).ToString("0") : Next i : sw.WriteLine(dat) '左から固有値の大きい順
        dat = "e-value" : For i = 0 To mm - 1 : dat = dat & "," & ev(i).ToString("0.000000") : Next i : sw.WriteLine(dat)  '固有値
        dat = "Prop" : For i = 0 To mm - 1 : dat = dat & "," & (ev(i) / tr).ToString("0.000000") : Next i : sw.WriteLine(dat) '寄与率
        For j = 0 To mm - 1
            dat = "e-vector" & (j + 1).ToString("0")
            For i = 0 To mm - 1
                dat = dat & "," & vec(j, i).ToString("0.000000")
            Next i
            sw.WriteLine(dat) '固有ベクトル
        Next j
        For k = 0 To ndata
            dat = "data" & (k + 1).ToString("0") : For i = 0 To mm - 1 : dat = dat & "," & zd(k, i).ToString("0.000000") : Next i : sw.WriteLine(dat) '各データの主成分
        Next k
        sw.Close()

    End Sub


    Private Sub JACOBI(ByVal n As Integer, ByRef a(,) As Double, ByRef ev() As Double, ByRef vec(,) As Double)
        'Jacobi法による固有値・固有ベクトル計算（奥村先生）
        '固有値ev(k),固有ベクトルvec(k,i)
        Dim i As Integer
        Dim j As Integer
        Dim k As Integer
        Dim s As Double
        Dim c As Double
        Dim aoff As Double
        Dim eps As Double
        Dim t As Double
        Dim aa As Double
        Dim aj As Double
        Dim ak As Double
        Dim vj As Double
        Dim vk As Double
        Dim vmax As Double
        Dim im As Integer
        Dim bb As Double

        For i = 0 To n - 1
            For j = 0 To n - 1
                vec(i, j) = 0.0
            Next j
            vec(i, i) = 1.0
        Next i

        aoff = 0.0
        For j = 0 To n - 2
            For k = j + 1 To n - 1
                aoff = aoff + a(j, k) * a(j, k)
            Next k
        Next j
        eps = aoff * 2.0
        For j = 0 To n - 1
            eps = eps + a(j, j) * a(j, j)
        Next j
        eps = eps * 0.00000000005

        Do While aoff > eps
            For j = 0 To n - 2
                For k = j + 1 To n - 1
                    If a(j, k) = 0.0 Then
                        t = 0.0
                    Else
                        aa = (a(k, k) - a(j, j)) / (2.0 * a(j, k))
                    End If
                    If aa >= 0 Then
                        t = 1.0 / (aa + Math.Sqrt(1.0 + aa * aa))
                    Else
                        t = 1.0 / (aa - Math.Sqrt(1.0 + aa * aa))
                    End If
                    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 To j - 1
                        aj = a(i, j)
                        ak = a(i, k)
                        a(i, j) = aj * c - ak * s
                        a(i, k) = aj * s + ak * c
                    Next i

                    For i = j + 1 To k - 1
                        aj = a(j, i)
                        ak = a(i, k)
                        a(j, i) = aj * c - ak * s
                        a(i, k) = aj * s + ak * c
                    Next i

                    For i = k + 1 To n - 1
                        aj = a(j, i)
                        ak = a(k, i)
                        a(j, i) = aj * c - ak * s
                        a(k, i) = aj * s + ak * c
                    Next i

                    For i = 0 To n - 1
                        vj = vec(i, j)
                        vk = vec(i, k)
                        vec(i, j) = vj * c - vk * s
                        vec(i, k) = vj * s + vk * c
                    Next i

                Next k
            Next j

            aoff = 0.0
            For j = 0 To n - 2
                For k = j + 1 To n - 1
                    aoff = aoff + a(j, k) * a(j, k)
                Next k
            Next j
        Loop

        For k = 0 To n - 1
            ev(k) = a(k, k)
        Next k

        For k = 0 To n - 2
            vmax = -1.0E+30
            For i = k To n - 1
                If ev(i) >= vmax Then
                    im = i : vmax = ev(i)
                End If
            Next i
            aa = ev(k) : ev(k) = ev(im) : ev(im) = aa
            For i = 0 To n - 1
                bb = vec(i, k)
                vec(i, k) = vec(i, im)
                vec(i, im) = bb
            Next i
        Next k

    End Sub


    Private Sub Form1_Load(ByVal sender As System.Object, ByVal e As System.EventArgs) Handles MyBase.Load
        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)
    End Sub
End Class