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
        '重判別分析（multiple discriminant analysis）
        Dim i As Integer
        Dim j As Integer
        Dim k As Integer
        Dim n As Integer
        Dim fname1 As String = ""
        Dim fname2 As String = ""

        Dim mm As Integer
        Dim ndata As Integer
        Dim nclass As Integer
        Dim ndf As Integer
        Dim xd(,) As Double
        Dim xm() As Double
        Dim temp() As Double
        Dim classmean(,) As Double
        Dim iclass() As Integer
        Dim members() As Integer

        Dim a(,) As Double
        Dim ev() As Double
        Dim vec(,) As Double
        Dim iflg As Integer
        Dim tr As Double
        Dim aa As Double
        Dim z(,) As Double
        Dim d(,) As Double

        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 str As String

        'データ読み込み
        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) : str = sbuf(0)
        dat = sr.ReadLine() : sbuf = dat.Split(delim) : mm = CInt(sbuf(0)) : ndata = CInt(sbuf(1))
        ReDim xd(ndata - 1, mm - 1)
        ReDim iclass(ndata - 1)
        ReDim members(ndata - 1)
        ReDim classmean(ndata - 1, mm - 1)
        ReDim xm(mm - 1)
        ReDim temp(mm - 1)
        ReDim a(mm - 1, mm - 1)
        ReDim ev(mm - 1)
        ReDim vec(mm - 1, mm - 1)
        k = 0
        Do While sr.Peek() >= 0
            dat = sr.ReadLine() : sbuf = dat.Split(delim)
            iclass(k) = CInt(sbuf(0))
            For i = 0 To mm - 1
                xd(k, i) = CDbl(sbuf(i + 1))
            Next i
            k = k + 1
        Loop
        sr.Close()
        ndata = k - 1

        'クラス毎データ組数（クラスは0から小さい順に並んでいる前提）
        nclass = 0
        For k = 0 To ndata
            If iclass(k) >= nclass Then nclass = iclass(k)
        Next k
        '平均値の計算
        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
            For k = 0 To ndata
                xd(k, i) = xd(k, i) - xm(i)
        Next k, i
        'クラス平均
        For k = 0 To ndata
            j = iclass(k)
            members(j) = members(j) + 1
            For i = 0 To mm - 1
                classmean(j, i) = classmean(j, i) + xd(k, i)
        Next i, k
        For i = 0 To nclass
            For j = 0 To mm - 1
                classmean(i, j) = classmean(i, j) / CDbl(members(i))
        Next j, i
        '固有方程式
        For i = 0 To mm - 1
            For j = i To mm - 1
                aa = 0.0
                For k = 0 To ndata
                    aa = aa + (xd(k, i) - classmean(iclass(k), i)) * (xd(k, j) - classmean(iclass(k), j))
                Next k
                a(i, j) = aa : a(j, i) = a(i, j)
        Next j, i
        '**********************************************
        iflg = 0 : Call MJACOBI(mm, a, ev, vec, iflg)
        '**********************************************
        ndf = ndata - nclass
        For i = 0 To mm - 1
            aa = Math.Sqrt(ev(i) / CDbl(ndf))
            For j = 0 To mm - 1
                vec(j, i) = vec(j, i) / aa
                a(i, j) = 0.0
        Next j, i
        For k = 0 To nclass
            For i = 0 To mm - 1
                temp(i) = 0.0
                For j = 0 To mm - 1
                    temp(i) = temp(i) + vec(j, i) * classmean(k, j)
            Next j, i
            For i = 0 To mm - 1
                For j = i To mm - 1
                    a(i, j) = a(i, j) + CDbl(members(k)) * temp(i) * temp(j)
                    a(j, i) = a(i, j)
            Next j, i
        Next k
        '**********************************************
        iflg = 1 : Call MJACOBI(mm, a, ev, vec, iflg)
        '**********************************************
        tr = 0.0
        For i = 0 To mm - 1
            tr = tr + ev(i)
        Next i

        ReDim z(ndata, mm - 1)
        ReDim d(ndata, nclass)
        'クラス平均座標表示
        For n = 0 To nclass
            For i = 0 To mm - 1
                z(n, i) = 0.0
                For j = 0 To mm - 1
                    z(n, i) = z(n, i) + vec(j, i) * classmean(n, j)
            Next j, i
        Next n
        For k = 0 To ndata
            'データ毎座標表示
            For i = 0 To mm - 1
                temp(i) = xd(k, i)
            Next i
            For i = 0 To mm - 1
                xd(k, i) = 0.0
                For j = 0 To mm - 1
                    xd(k, i) = xd(k, i) + vec(j, i) * temp(j)
            Next j, i
            'クラス平均からの距離計算
            For n = 0 To nclass
                d(k, n) = 0.0
                For j = 0 To mm - 1
                    d(k, n) = d(k, n) + (xd(k, j) - z(n, j)) * (xd(k, j) - z(n, j))
                Next j
                d(k, n) = Math.Sqrt(d(k, n))
            Next n
        Next k

        'データ出力
        If SaveFileDialog1.ShowDialog() = Windows.Forms.DialogResult.OK Then fname2 = SaveFileDialog1.FileName
        sw = New System.IO.StreamWriter(fname2, False, System.Text.Encoding.Default)
        dat = str : sw.WriteLine(dat)
        dat = "No," : For i = 0 To mm - 1 : dat = dat & "," & CStr(i + 1) : Next i : sw.WriteLine(dat)
        dat = "固有値," : For i = 0 To mm - 1 : dat = dat & "," & ev(i).ToString("0.000") : Next i : sw.WriteLine(dat)
        dat = "寄与率," : For i = 0 To mm - 1 : dat = dat & "," & (ev(i) / tr).ToString("0.000") : Next i : sw.WriteLine(dat)
        dat = "クラス平均正準座標" : sw.WriteLine(dat)
        For n = 0 To nclass
            dat = "size=" & members(n) & ",class" & CStr(n)
            For i = 0 To mm - 1
                dat = dat & "," & z(n, i).ToString("0.000")
            Next i
            sw.WriteLine(dat)
        Next n
        dat = "データ毎正準座標,"
        For i = 0 To mm - 1
            dat = dat & "," & CStr(i + 1)
        Next i
        For n = 0 To nclass
            dat = dat & ",dis-c" & CStr(n)
        Next n
        sw.WriteLine(dat)
        For k = 0 To ndata
            dat = "data" & (k + 1).ToString("00") & ",class" & CStr(iclass(k))
            For i = 0 To mm - 1
                dat = dat & "," & xd(k, i).ToString("0.000")
            Next i
            For n = 0 To nclass
                dat = dat & "," & d(k, n).ToString("0.000")
            Next n
            sw.WriteLine(dat)
        Next k
        sw.Close()

        '新規データ評価実施判断
        Dim newcase As String = ""
        Select Case MessageBox.Show("新規データ評価実施", "ファイル読込", MessageBoxButtons.YesNo, MessageBoxIcon.Question)
            Case Windows.Forms.DialogResult.Yes : newcase = "yes"
            Case Windows.Forms.DialogResult.No : newcase = "no"
        End Select
        If newcase = "yes" Then
            Dim kdata As Integer
            Dim xnew(,) 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()
            dat = sr.ReadLine() : sbuf = dat.Split(delim) : mm = CInt(sbuf(0)) : kdata = CInt(sbuf(1))
            ReDim xnew(kdata - 1, mm - 1)
            k = 0
            Do While sr.Peek() >= 0
                dat = sr.ReadLine() : sbuf = dat.Split(delim)
                For i = 0 To mm - 1
                    xnew(k, i) = CDbl(sbuf(i))
                Next i
                k = k + 1
            Loop
            sr.Close()
            kdata = k - 1
            For k = 0 To kdata
                '座標計算
                For i = 0 To mm - 1
                    z(k, i) = 0.0
                    For j = 0 To mm - 1
                        z(k, i) = z(k, i) + vec(j, i) * (xnew(k, j) - xm(j))
                Next j, i
                'クラス平均からの距離計算
                For n = 0 To nclass
                    d(k, n) = 0.0
                    For i = 0 To mm - 1
                        aa = 0.0
                        For j = 0 To mm - 1
                            aa = aa + vec(j, i) * (xnew(k, j) - xm(j) - classmean(n, j))
                        Next j
                        d(k, n) = d(k, n) + aa * aa
                    Next i
                    d(k, n) = Math.Sqrt(d(k, n))
                Next n
            Next k

            'データ追加出力
            sw = New System.IO.StreamWriter(fname2, True, System.Text.Encoding.Default)
            If newcase = "yes" Then
                '評価データ諸元書き込み
                dat = "新規評価データ正準座標" : sw.WriteLine(dat)
                For k = 0 To kdata
                    dat = "new," & k.ToString("00")
                    For i = 0 To mm - 1
                        dat = dat & "," & z(k, i).ToString("0.000")
                    Next i
                    For n = 0 To nclass
                        dat = dat & "," & d(k, n).ToString("0.000")
                    Next n
                    sw.WriteLine(dat)
                Next k
            End If
            sw.Close()
        End If

    End Sub


    Private Sub MJACOBI(ByVal n As Integer, ByRef a(,) As Double, ByRef ev() As Double, ByRef vec(,) As Double, ByVal iflg As Integer)
        '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

        If iflg = 0 Then '固有ベクトル初期化
            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
        End If

        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


End Class