' ============================================================================
' ===   ベイズ最小二乗法（BLS）プログラム：点滴２ーコンパートメントモデル  ===
' ===   Ver.1.04     by MyAngel-labo.   https://www.myangel-labo.com        ===
' ===   このソフトウェアを用いて生じるいかなる問題点についても             ===
' ===   製作者は一切の責任を負いませんので、あらかじめご了承ください。     ===
' ============================================================================
' このプログラムは「VCM-TDM」の機能の参考に作成しました。変数名やアルゴリズム
' は適宜改訂しています。また今後の開発のために現時点で不必要な部分もあります。

Option Explicit
' グローバル変数の配列宣言
' フラグ
Dim flagCheckData As Integer, flagCheckCalc As Integer, flagCheckRecom As Integer
Dim flagR As Integer, PopParamSet As Integer

' メッセージ文
Dim temp As String, Temp2 As VbMsgBoxResult

' 患者データ
Dim GEN As Integer, AGE As Single, BWT As Single, CLcr As Single

' 薬物動態パラメータ
Dim P(20) As Single, Pref(20) As Single, AL As Single, BE As Single
Dim THETA(20) As Single, OMEGA(20) As Single, SIGMA(20) As Single
Dim PM(20) As Single, OM(20) As Single, SG(20) As Single
Dim Ke As Single, Vc As Single, Vss As Single, K12 As Single, K21 As Single

' モデルの定義
Dim nModel As Integer, nComp As Integer, nRoute As Integer, nReg As Integer
Dim nParam As Integer, nRef As Integer

' 最小二乗法の種類（通常、ベイズ最小二乗法）
Dim nMethod As Integer:  ' nMethod - (1: OLS), (2: BLS)

' 分布容積の定義（Vc と Vss のどちらを推定パラメータとするか）
Dim nVss As Integer:     ' nVss - (1: Vc, 2: Vss)

' 投与・採血履歴関連
Dim maxrec As Integer
Dim recNo(200) As Single, recDate(200) As String, recTime(200) As String
Dim recT(200) As Single, recEVID(200) As Integer
Dim recRoute(200) As Integer, recDose(200) As Single, recTinf(200) As Single
Dim recDN(200) As Integer, recII(200) As Single, recCp(200) As Single
Dim v0(200) As Single, vv(200, 10) As Single:  ' v0 ... dosing number

' 薬物濃度データ関連
Dim nDose As Single, recN8 As Integer, recN9 As Integer
Dim tx(20) As Single, cy(20) As Single, tx8(10) As Single, cy8(10) As Single
Dim Cp As Single, SS As Single, tn As Integer
        
Sub start_Dosage()
   Dim i As Integer, j As Integer
   
   flagCheckData = 0: ' データが正しく入力されているかのフラフをリセット
   Call ClearCells: ' 計算開始にあたり関連セルの表示をクリア
   Call set_Model: ' モデルの設定
   Call set_paramName: ' パラメータ名の設定
   Call set_Individual: ' 患者背景データの読み込み
   
   ' データが正しく入力されていない場合、フラグを 1 にして先に進まない
   If flagCheckData = 1 Then Exit Sub
   
   ' 既に推奨用量の設定が行われた跡がある場合、それを削除
   For i = 1 To 200
      If Cells(67 + i, 13).Value = "*" Then
         For j = 3 To 13: Cells(67 + i, j) = "": Next j
      End If
   Next i
   
   ' 血中濃度計算値のセルをクリア
   For i = 1 To 200: Cells(67 + i, 13) = "": Next i
   ' セルを移動
   Cells(63, 3).Select
End Sub

Sub ClearCells()
   Dim i As Integer, j As Integer
   
   ' 関連セルの表示をリセット
   Cells(6, 17) = ""
   For i = 1 To 2: For j = 1 To 5: Cells(6 + j, 16 + i).Value = "": Next j: Next i
   For i = 1 To 3: Cells(14 + i, 17).Value = "": Next i
   For i = 1 To 4: For j = 1 To 8: Cells(20 + j, 15 + i).Value = "": Next j: Next i
   For i = 1 To 5: For j = 1 To 5: Cells(30 + j, 15 + i).Value = "": Next j: Next i
   i = j
End Sub

Sub set_Model()
   
   ' モデルの定義
   nComp = 2: nRoute = 4: nReg = 2: ' 2-Comp. 点滴モデル（不定用量・不定間隔）
   nModel = nComp * 100 + nRoute * 10 + nReg
   
   ' 母集団パラメータの設定
   PopParamSet = Cells(52, 16).Value
   If PopParamSet <> 1 And PopParamSet <> 2 Then
      temp = MsgBox("母集団パラメータの準備ができていません", vbOKCancel + vbExclamation)
      nComp = 2
      Exit Sub
   End If
   If PopParamSet = 1 Then nVss = 1: ' Vc  を推定パラメータとして用いる
   If PopParamSet = 2 Then nVss = 2: ' Vss を推定パラメータとして用いる
   
   nParam = 4: nRef = 3: ' 推定パラメータ4つ、二次パラメータ3つ。
   nMethod = Cells(5, 5) ' 1: OLS, 2: BLS
End Sub

Sub set_paramName()
    Dim i As Integer, Pname(20) As String, Rname(20) As String
    
    ' 推定パラメータ名、二次パラメータ名を設定
    Pname(1) = "CL (L/hr)": Pname(2) = "K12 (/hr)": Pname(3) = "K21 (/hr)"
    Rname(1) = "*Ke (/hr)": Rname(2) = "*t1/2b (hr)"
    If nVss = 1 Then Pname(4) = "Vc (L)": Rname(3) = "*Vss (L)"
    If nVss = 2 Then Pname(4) = "Vss (L)": Rname(3) = "*Vc (L)"
    For i = 1 To 4: Cells(52 + i, 9) = Pname(i): Next i
    For i = 1 To 3: Cells(52 + i, 12) = Rname(i): Next i
End Sub

Sub set_Individual()
   Dim i As Integer
   
   ' 患者背景データの読み込み
   ' ID = Cells(7, 5).Value
   GEN = Cells(8, 5).Value: AGE = Cells(8, 8).Value: BWT = Cells(8, 11).Value
   
   ' クレアチニンクリアランスの計算方法を選択 -> 1: 実測値, 2: Cockcroft, 3: eGFR
   CLcr = 0
   For i = 1 To 3:
      If Cells(9 + i, 11) = 1 Then CLcr = Cells(9 + i, 10).Value
   Next i
   
   ' エラーチェック
   If GEN <> 1 And GEN <> 2 Then
      Temp2 = MsgBox("性別が正しく設定されていません", vbOKOnly + vbExclamation)
      flagCheckData = 1: Exit Sub
   End If
   If AGE <= 0 Then
      Temp2 = MsgBox("年齢が正しく設定されていません", vbOKOnly + vbExclamation)
      flagCheckData = 1: Exit Sub
   End If
   If BWT <= 0 Then
      Temp2 = MsgBox("体重が正しく設定されていません", vbOKOnly + vbExclamation)
      flagCheckData = 1: Exit Sub
   End If
   If CLcr <= 0 Then
      Temp2 = MsgBox("CLcr が正しく設定されていません", vbOKOnly + vbExclamation)
      flagCheckData = 1: Exit Sub
   End If
End Sub

Sub set_Dosage()
   Dim i As Integer, j As Integer, k As Integer, maxT As Single
   Dim dosing As Integer, count9 As Integer, count8 As Integer
   
   flagCheckData = 0
   maxrec = 0
   
   ' ベイズ推定のための母集団パラメータ設定
   If nMethod = 2 Then Call set_PopParam(BWT, CLcr)
   
   Dim iBase As Integer: iBase = 67
   ' 投与・採血記録の読み取り、最大行数 maxrec
   For i = 1 To 200
      k = Cells(iBase + i, 6).Value
      If ((k = 1) Or (k = 2) Or (k = 3) Or (k = 8) Or (k = 9)) Then
         maxrec = i
      Else
         If k <> 0 Then
            Temp2 = MsgBox("EVID 行が正しく設定されていません", vbOKOnly + vbExclamation)
            flagCheckData = 1: Exit Sub
         End If
      End If
   Next i
   If maxrec <= 0 Then
      Temp2 = MsgBox("投与の記録がありません", vbOKOnly + vbExclamation)
      flagCheckData = 1: Exit Sub
   End If

   ' 投与・採血記録を変数 recXX[] に読み込み
   For i = 1 To maxrec
      recDate(i) = Cells(iBase + i, 3).Value:                ' 日付：Column C
      If recDate(i) = "" Then recDate(i) = recDate(i - 1)
      recTime(i) = Cells(iBase + i, 4).Value:                ' 時刻：Column D
      
      ' 日付、時刻の入力ができているか確認
      If recDate(1) = "" Then
         Temp2 = MsgBox("一行目の日付が設定されていません", vbOKOnly + vbExclamation)
        flagCheckData = 1: Exit Sub
      End If
      ' 日付は空白の場合もあり得るので１行目の記録のみ空白でないかを確認する
      If recTime(i) = "" Then
         Temp2 = MsgBox("時刻が設定されていません", vbOKOnly + vbExclamation)
         flagCheckData = 1: Exit Sub
      End If
      
      Call calc_Time(i): ' 日付と時刻データをもとに、初回投与からの経過時間を計算
      If flagCheckData = 1 Then Exit Sub
      
      Cells(iBase + i, 5).Value = recT(i):                   ' 経過時間：Column E
      recEVID(i) = Cells(iBase + i, 6).Value:                ' EVID：Column F
      If recEVID(i) = 1 Then
         recRoute(i) = Cells(iBase + i, 7).Value:            ' 投与経路：Column G
         recDose(i) = Cells(iBase + i, 8).Value:             ' 投与量：Column H
         recTinf(i) = Cells(iBase + i, 9).Value:             ' 点滴時間：Column I
         recDN(i) = Cells(iBase + i, 10).Value:              ' 繰り返し回数：Column J
         If recDN(i) = 0 Then
            recDN(i) = 1
         Else
            recII(i) = Cells(iBase + i, 11).Value:           ' 投与間隔；Column K
         End If
         recCp(i) = 0: ' 念のため濃度測定値を 0 にしておく
      End If
      
      ' 採血記録の場合
      If ((recEVID(i) = 8) Or (recEVID(i) = 9)) Then
         recRoute(i) = 0: recDose(i) = 0: recTinf(i) = 0     ' 投与に関するデータをクリア
         recDN(i) = 0: recII(i) = 0
         recCp(i) = Cells(iBase + i, 12).Value: ' 血中濃度の読み込み
         If recCp(i) <= 0 Then
            Temp2 = MsgBox("血中濃度データが0または負の値です", vbOKOnly + vbExclamation)
            flagCheckData = 1: Exit Sub
         End If
      End If
   Next i

   ' Step 1 読み込んだ情報の確認
   nReg = 2: ' このプログラムでは不定用量・不定間隔の場合のみ使用可
   For i = 1 To maxrec
      flagCheckData = 0: ' エラーフラグを一旦リセット
      If recEVID(i) = 1 Then
         ' 投与条件の入力確認
         If recRoute(i) = 0 Then temp = "投与経路が設定されていません": flagCheckData = 1
         If recRoute(i) <> 4 Then temp = "点滴以外の経路が選択されています": flagCheckData = 1
         If ((recDN(i) > 1) And (recII(i) <= 0)) Then
            temp = "投与間隔が設定されていません": flagCheckData = 1
         End If
         If (recDose(i) = 0) Then temp = "投与量が設定されていません": flagCheckData = 1
         If ((recRoute(i) = 4) And (recTinf(i) = 0)) Then
            temp = "点滴時間が設定されていません": flagCheckData = 1
         End If
      End If
      If ((recEVID(i) = 8) Or (recEVID(i) = 9)) Then temp = ""
      If ((recEVID(i) = 2) Or (recEVID(i) = 3)) Then
         temp = "このプログラムでは EVID = 1のみ使えます": flagCheckData = 1
      End If
      
      ' メッセージの表示
      If flagCheckData = 1 Then
         temp = temp + " Line = " + Str$(iBase + i) + "."
         Temp2 = MsgBox(temp, vbOKOnly + vbExclamation)
         Exit Sub
      End If
   Next i
                
   ' Step 2 - 読み込んだ情報を指定配列に代入
   dosing = 0: count9 = 0: count8 = 0: ' 投与回数、血中濃度データ数等を初期化
   For i = 1 To maxrec
   ' 1(time), 2(dose), 3(Tinf), 4(Rate), 5(II), 6(NA), 7(route)
      If recEVID(i) = 1 Then
         For j = 1 To recDN(i)
            dosing = dosing + 1
            v0(dosing) = i:                               ' v0() 投与回
            vv(dosing, 1) = recT(i) + recII(i) * (j - 1): ' vv(*, 1) 投与時間
            vv(dosing, 2) = recDose(i):                   ' vv(*, 2) 投与量
            vv(dosing, 3) = recTinf(i):                   ' vv(*, 3) 点滴時間
            If recTinf(i) <> 0 Then
               vv(dosing, 4) = recDose(i) / recTinf(i):   ' vv(*, 4) 点滴速度
            Else
               vv(dosing, 4) = 0
            End If
            vv(dosing, 5) = recII(i):                     ' vv(*, 5) 投与間隔
            If recRoute(i) = 4 Then vv(dosing, 7) = 4:    ' vv(*, 7) 投与経路（4のみ）
         Next j
      Else
         If recEVID(i) = 8 Then
            count8 = count8 + 1: ' EVID = 8 解析には用いない血中濃度測定値
            tx8(count8) = recT(i): cy8(count8) = recCp(i)
            recN8 = count8
         End If
         If recEVID(i) = 9 Then
            count9 = count9 + 1 ' EVID = 9 解析に用いる血中濃度測定値
            tx(count9) = recT(i): cy(count9) = recCp(i)
            recN9 = count9
         End If
      End If
   Next i
   nDose = dosing: ' 全投与回数

   ' Step 3 - 配列変数 v0 と vv を投与時間順に並べ替え
   Dim work As String
   For i = 1 To nDose
      For j = i + 1 To nDose
         If vv(i, 1) > vv(j, 1) Then
            work = vv(i, 1): vv(i, 1) = vv(j, 1): vv(j, 1) = work
            For k = 2 To 9
               work = vv(i, k): vv(i, k) = vv(j, k): vv(j, k) = work
            Next k
            work = v0(i): v0(i) = v0(j): v0(j) = work
         End If
      Next j
   Next i
        
   maxT = vv(nDose, 1) - vv(1, 1): ' 投与時刻の最大値
   If tx(recN9) >= maxT Then maxT = tx(recN9)
   If tx8(recN8) >= maxT Then maxT = tx8(recN8)
   Cells(46, 11).Value = maxT + 24: ' グラフ横軸最大値を最終データより24時間大きくしておく
  
   ' 投与間隔を投与ごとに設定
   For i = 2 To nDose + 1: vv(i, 5) = vv(i, 1) - vv(i - 1, 1): Next i
   
   ' 推奨用量計算時には次のコメントは表示しない
   If flagR <> 1 Then
      Temp2 = MsgBox("投与・採血記録が正しく設定できました", vbOKOnly + vbInformation)
   End If
   
   ' == 投与・採血履歴データの確認（シート「Check」に表示される ==
   GoTo Jump:
   
   ' 確認時以外は以下は実行しない（GoTo Jump:を有効にする）
   Sheets("Check").Select
   Cells(2, 2) = "Check data": Cells(2, 4) = nDose
   For i = 1 To nDose
      Cells(2 + i, 2) = vv(i, 1): Cells(2 + i, 3) = vv(i, 2)
      Cells(2 + i, 4) = vv(i, 3): Cells(2 + i, 5) = vv(i, 4)
      Cells(2 + i, 6) = vv(i, 5): Cells(2 + i, 7) = vv(i, 7)
   Next i
   For i = 1 To recN9
      Cells(2 + i, 9) = tx(i): Cells(2 + i, 10) = cy(i)
      Cells(2 + i, 11) = recN9
   Next i
   For i = 1 To recN8
      Cells(2 + recN9 + i, 9) = tx8(i): Cells(2 + recN9 + i, 10) = cy8(i)
      Cells(2 + recN9 + i, 11) = recN8
   Next i
   
Jump:
   flagCheckRecom = 0: ' 推奨用量設定フラグ、9 で未設定
End Sub

Sub MoveTop()
   ' 画面移動
   Cells(3, 3).Select
   Exit Sub

End Sub

Sub calc_Time(i)
   Dim yy0 As Single, mn0 As Single, dd0 As Single
   Dim yy1 As Single, mn1 As Single, dd1 As Single
   Dim hh0 As Single, mi0 As Single, hh1 As Single, mi1 As Single
   Dim Ddiff As Single, Tdiff As Single
   
   ' 初回からの投与時間の算出
   recT(1) = 0
   If i = 1 Then
      recT(i) = 0
   Else
      yy0 = Year(recDate(1)): mn0 = Month(recDate(1)): dd0 = Day(recDate(1))
      hh0 = Int(Val(recTime(1)) * 24): mi0 = (Val(recTime(1)) * 24 - hh0) * 60
      yy1 = Year(recDate(i)): mn1 = Month(recDate(i)): dd1 = Day(recDate(i))
      hh1 = Int(Val(recTime(i)) * 24): mi1 = (Val(recTime(i)) * 24 - hh1) * 60
            
      Ddiff = DateSerial(yy1, mn1, dd1) - DateSerial(yy0, mn0, dd0)
      Tdiff = TimeSerial(hh1, mi1, 0) - TimeSerial(hh0, mi0, 0)
      Tdiff = (Ddiff + Tdiff) * 24
      
      ' エラーチェック
      If Tdiff < 0 Then
         Temp2 = MsgBox("経過時間が負の数となります", vbOKOnly + vbExclamation)
         flagCheckData = 1: Exit Sub
      End If
      recT(i) = recT(1) + Tdiff
   End If
End Sub

Sub set_PopParam(BWT, CLcr):  ' Rodvold
   ' EXCEL セルに PPK パラメータが正しく定義されている必要がある
   ' 母集団パラメータ値を追加するときは PopParamSetによる場合分けを増やす
   
   If PopParamSet = 1 Then
      ' PPK by Rodvold, Two-Comp. Inf
      THETA(1) = Cells(53, 4).Value ' Slope for WT
      THETA(2) = Cells(54, 4).Value ' Slope for CLcr
      THETA(3) = Cells(55, 4).Value ' K12
      THETA(4) = Cells(56, 4).Value ' K21
      THETA(5) = Cells(57, 4).Value ' Vc
      OMEGA(1) = Cells(53, 6).Value ' CV% for CL
      OMEGA(2) = Cells(54, 6).Value ' CV% for K12
      OMEGA(3) = Cells(55, 6).Value ' CV% for K21
      OMEGA(4) = Cells(56, 6).Value ' CV% for Vc
      SIGMA(1) = Cells(53, 8).Value ' Absolute error
      SIGMA(2) = Cells(54, 8).Value ' Relative error CV%
      ' set population mean, OM, SG
      PM(1) = THETA(1) * BWT + THETA(2) * CLcr ' CL
      PM(2) = THETA(3)                         ' K12
      PM(3) = THETA(4)                         ' K12
      PM(4) = THETA(5) * BWT                   ' Vc
      OM(1) = OMEGA(1)                         ' CL
      OM(2) = OMEGA(2)                         ' K12
      OM(3) = OMEGA(3)                         ' K21
      OM(4) = OMEGA(4)                         ' Vc
      SG(1) = 0 ' 個体内変動分散は obj_func ルーチン内で計算される
   
   ElseIf PopParamSet = 2 Then
      ' PPK by Yasuhara, Two-Comp. Inf
      THETA(1) = Cells(53, 4).Value ' Slope for WT
      THETA(2) = Cells(54, 4).Value ' Slope for CLcr
      THETA(3) = Cells(55, 4).Value ' K12
      THETA(4) = Cells(56, 4).Value ' K21
      THETA(5) = Cells(57, 4).Value ' Vc
      OMEGA(1) = Cells(53, 6).Value ' CV% for CL
      OMEGA(2) = Cells(54, 6).Value ' CV% for K12
      OMEGA(3) = Cells(55, 6).Value ' CV% for K21
      OMEGA(4) = Cells(56, 6).Value ' CV% for Vc
      SIGMA(1) = Cells(53, 8).Value ' Absolute error
      SIGMA(2) = Cells(54, 8).Value ' Relative error CV%
      ' set population mean, OM, SG
      PM(1) = THETA(1) * CLcr                  ' CL
      PM(2) = THETA(3)                         ' K12
      PM(3) = THETA(4)                         ' K12
      PM(4) = THETA(5)     ' * Vss             ' Vc
      OM(1) = OMEGA(1)                         ' CL
      OM(2) = OMEGA(2)                         ' K12
      OM(3) = OMEGA(3)                         ' K21
      OM(4) = OMEGA(4)     ' * Vss             ' Vc
      SG(1) = 0 ' 個体内変動分散は obj_func ルーチン内で計算される
   End If
   
   Call ShowPmean
End Sub

Sub ShowPmean()
   Dim i As Integer
   ' 母集団平均値（THETA、OMEGA、SIGMA）の表示
   For i = 1 To nParam: Cells(52 + i, 10).Value = PM(i): Next i
End Sub

Sub CopyTo()
    
    ' 投与・採血データの別シートへのコピー（保存）
    Range("C66:M113").Select
    Selection.Copy
    Sheets.Add After:=ActiveSheet
    Range("C66").Select
    ActiveSheet.Paste

    Columns("C:M").Select
    Range("C56").Activate
    Selection.ColumnWidth = 12

End Sub

Sub CopyFrom()
    Dim SheetName As String
    
    ' 指定したシートから投与・採血データをコピー（読み込み）
    SheetName = InputBox("読み込みたいシート名 =")
    On Error GoTo Ertrap:
    
    Sheets(SheetName).Select
    Range("C66:M113").Select
    Selection.Copy
    Sheets("Main").Select
    Range("C66").Select
    ActiveSheet.Paste
    Exit Sub
    
Ertrap:
    temp = MsgBox("読み込み時エラー（シート名の確認等が必要です）", vbOKOnly + vbExclamation)
    Exit Sub
End Sub

Sub DataErase()
   
   ' 投与・採血履歴をクリアする
   temp = MsgBox("データを消去します。元に戻せませんがよろしいですか？", vbOKCancel + vbExclamation)
   If temp = vbOK Then
       Range("C68:M113").Select
       Selection.ClearContents
   End If
   
End Sub

'=== ベイズ最小二乗法計算 ===
Sub goMULTI()
   Dim i As Integer, j As Integer
   
   flagCheckCalc = 0: ' あてはめ計算済フラグを一旦リセット
   ' するので
   ' 一旦推奨用法の計算をした場合、再度データ設定を行ってからあてはめ計算を行なう必要がある
   If flagCheckRecom = 1 Then
      Temp2 = MsgBox("再度、投与・採血履歴を設定してから実行してください", vbOKOnly + vbExclamation)
      Exit Sub
   End If
   '
   If flagCheckData = 1 Then
      Temp2 = MsgBox("投与・採血履歴が正しく設定されていません", vbOKOnly)
      Exit Sub
   End If
   '
   If (nMethod = 1 And recN9 < nParam) Then
      Temp2 = MsgBox("OLS 計算ではパラメータ数以上のデータ数が必要です", vbOKOnly)
      flagCheckData = 1: Exit Sub
   End If
   
   Range("U5", "AB2004").Value = "": ' グラフ描画のためのデータをクリア
   
   ' 初期値の設定
   Call set_InitParam:
   If flagCheckCalc = 1 Then
      Temp2 = MsgBox("初期値設定にエラーがあります", vbOKOnly + vbExclamation)
      Exit Sub
   End If
   
   ' 初期値の表示
   For i = 1 To nParam: Cells(7 + i, 18).Value = P(i): Next i
    
   'Simplex 法
   Call simplex
   If flagCheckCalc = 1 Then
      Temp2 = MsgBox("ベイズ最小二乗法の過程でエラーがありました", vbOKOnly + vbExclamation)
      Exit Sub
   End If
    
   ' 二次パラメータの計算と表示
   Call calc_refine:
   If flagCheckCalc = 1 Then
      Temp2 = MsgBox("二次パラメータの計算課程でエラーがありました", vbOKOnly + vbExclamation)
      Exit Sub
   End If
   
   Call print_output: ' 結果の表示
   Call plot_output: ' 結果のグラフ表示
   Temp2 = MsgBox("ベイズ最小二乗法計算が正しく終了しました", vbOKOnly + vbInformation)
   
End Sub

Sub set_InitParam()
   Dim i As Integer
   
   ' 初期値の読み込み
   If Cells(61, 11) = 1 Then
      ' 母集団パラメータ値を利用
      For i = 1 To nParam: P(i) = PM(i): Cells(52 + i, 10) = P(i): Next i
   Else
      ' 入力値を利用
      For i = 1 To nParam: P(i) = Cells(52 + i, 10).Value: Next i
   End If
   
   ' エラーチェック（負の数の場合エラーとする）
   For i = 1 To nParam
      If P(i) <= 0 Then flagCheckCalc = 1: Exit Sub
   Next i
End Sub

Sub simplex()
   ' Nelder-Mead method (Simplex method)
   ' Coded originally by K. Yamaoka, J.P.Dyn. 1985,8,246-256.
   ' このサブルーチンは山岡先生のプログラムリストをほぼそのまま転用した
   ' LINE XX in MULTI2(BAYES) とは、そのプログラムの行番号を示す。
   Dim i As Integer, j As Integer, kk As Integer, M As Integer
   Dim A1 As Integer, B1 As Integer, C1 As Integer
   Dim JH As Integer, JL As Integer, JS As Integer
   Dim SR As Single, SL As Single, SG As Single, PC As Single
   Dim A(20, 21) As Single, X0(20) As Single
   Dim XR(20) As Single, X1(20) As Single
   
   flagCheckCalc = 0
   Randomize Val(Right$(Time$, 2)): ' 乱数初期化
   A1 = 1: B1 = 0.5: C1 = 2: ' Line 3000 in MULTI2(BAYES)
   SG = 10000000000#: SL = 0: PC = 0.0001
   M = nParam: ' パラメータ数
   
   ' 初期値を配列に代入
   For i = 1 To M: A(i, 1) = P(i): Next i
   
   ' 乱数によりシンプレックス面を作成
   For j = 2 To M + 1
      For i = 1 To M
Label_L1:
         A(i, j) = 2 * Rnd(1) * A(i, 1) + 0.01 * (Rnd(1) - 0.5)
         If A(i, j) <= 0 Then GoTo Label_L1:
      Next i
   Next j
   
   kk = 0
   ' 初期値とその時の残差平方和を表示・確認
   Cells(6, 17).Value = ""
   ' シンプレックス面でのパラメータ値と残差平方和
   For j = 1 To M + 1
      For i = 1 To M
         P(i) = A(i, j): Cells(31 + i, 15 + j).Value = A(i, j)
      Next i
      A(0, j) = obj_func: If flagCheckCalc = 1 Then Exit Sub
      Cells(7, 18).Value = A(0, j)
      Cells(31, 15 + j).Value = A(0, j)
   Next j
   
   ' シンプレクス法ループ開始 - Line 3070 in MULTI2(BAYES)
Label_3070:
   kk = kk + 1
   If kk >= 500 Then
      Temp2 = MsgBox("反復回数が 500 を超えました。再計算してください", vbOKOnly)
      flagCheckCalc = 1: Exit Sub
   End If
   
   Cells(30, 19).Value = kk
   For j = 1 To M + 1
      Cells(31, 15 + j).Value = A(0, i)
      For i = 1 To M: Cells(31 + i, 15 + j).Value = A(i, 1): Next i
   Next j
   
   ' 収束したか確認, Line 3077 in MULTI2(BAYES)
   GoTo Label_5000:

Label_3080:
   SR = 0: SL = 10000000000#: ' 残差平方和の最大値、最小値を得る
   For j = 1 To M + 1
      If SR < A(0, j) Then JH = j: SR = A(0, j)
      If SL > A(0, j) Then JL = j: SL = A(0, j)
   Next j
   
   ' 2番目に大きい残差平方和を得る, Line 3100 in MULTI2(BAYES)
   SR = 0:
   For j = 1 To M + 1
      If (j <> JH) And (SR < A(0, j)) Then
         JS = j: SR = A(0, j)
      End If
   Next j
   
   ' 残差平方和が最大値となるセットを除く、各パラメータ平均値を求める
   For i = 1 To M
      X0(i) = 0
      For j = 1 To M + 1
         If j <> JH Then X0(i) = X0(i) + A(i, j)
      Next j
      X0(i) = X0(i) / M
   Next i
   
   ' Line 3120 in MUTI2(BAYES)
   For i = 1 To M
      A(i, 0) = (1 + A1) * X0(i) - A1 * A(i, JH): P(i) = A(i, 0)
   Next i
   SR = obj_func: If flagCheckCalc = 1 Then Exit Sub

   ' Line 3130 in MULTI2(BAYES)
   If SR <= A(0, JS) Then GoTo Label_3300:
   If SR <= A(0, JH) Then
      For i = 1 To M: A(i, JH) = A(i, 0): Next i
      A(0, JH) = SR
   End If
   
   ' Line 3170 in MULTI2(BAYES)
   For i = 1 To M
      A(i, 0) = B1 * A(i, JH) + (1 - B1) * X0(i): P(i) = A(i, 0)
   Next i
   SR = obj_func: If flagCheckCalc = 1 Then Exit Sub
   
   If SR <= A(0, JH) Then
      For i = 1 To M: A(i, JH) = A(i, 0): Next i
      A(0, JH) = SR: GoTo Label_3070:
   End If
   
   ' Line 3200 in MULTI2(BAYES)
   For j = 1 To M + 1
      For i = 1 To M
         A(i, j) = (A(i, j) + A(i, JL)) / 2: P(i) = A(i, j)
      Next i
      A(0, j) = obj_func: If flagCheckCalc = 1 Then Exit Sub
   Next j
   GoTo Label_3070:
      
Label_3300:
   If SR < A(0, JL) Then GoTo Label_3500:
Label_3320:
   For i = 1 To M: A(i, JH) = A(i, 0): Next i
   A(0, JH) = SR: GoTo Label_3070
      
Label_3500:
   ' Line 3500 in MULTI2(BAYES)
   For i = 1 To M
      X1(i) = C1 * A(i, 0) + (1 - C1) * X0(i): P(i) = X1(i)
   Next i
   SL = obj_func: If flagCheckCalc = 1 Then Exit Sub
   
   If SL < SR Then
      For i = 1 To M: A(i, JH) = X1(i): Next i
      A(0, JH) = SL: GoTo Label_3070
   End If
   GoTo Label_3320:

Label_5000:
   ' 収束したか判定 Line 5000 in MULTI2(BAYES)
   SR = 0
   For i = 1 To M + 1: SR = SR + A(0, i): Next i
   If Abs(SR - SG) > PC * SG Then
      SG = SR
      GoTo Label_3080
   Else
      For i = 1 To M
         P(i) = A(i, JL): Cells(7 + i, 17).Value = P(i)
      Next i
      Cells(6, 17).Value = kk
      Cells(7, 17).Value = A(0, JL)
      Exit Sub
   End If

End Sub

Function obj_func() As Single
   Dim i As Integer, tt As Single, Cp As Single, IW As Integer, S2 As Single
   ' 残差平方和（目的関数の値）の計算
   ' 相対誤差モデルであるが分母を CV値×平均値 = SD で近似している
   
   SS = 0
   
   ' OLS（通常の最小二乗法）
   If nMethod = 1 Then
      IW = Cells(5, 9).Value
      For i = 1 To recN9
         If cy(i) = 0 Then
            Temp2 = MsgBox("血中濃度データが 0 です", vbOKOnly)
            flagCheckCalc = 1: Exit Function
         End If
         
         Cp = calc_Cp(tx(i))
         If flagCheckCalc = 1 Then Exit Function
         
         If IW = 0 Then SS = SS + (cy(i) - Cp) ^ 2 / 1
         If IW = 1 Then SS = SS + (cy(i) - Cp) ^ 2 / cy(i)
         If IW = 2 Then SS = SS + (cy(i) - Cp) ^ 2 / cy(i) / cy(i)
      Next i
   Else
   
   ' BLS（ベイズ最小二乗法）
      For i = 1 To recN9
         Cp = calc_Cp(tx(i))
         If flagCheckCalc = 1 Then Exit Function
         
         SG(1) = SIGMA(1) + SIGMA(2) * Cp
         SS = SS + (cy(i) - Cp) ^ 2 / SG(1) / SG(1)
      Next i
      S2 = 0
      For i = 1 To nParam
         S2 = S2 + ((PM(i) - P(i)) ^ 2) / ((OM(i) * PM(i)) ^ 2)
      Next i
      SS = SS + S2
   End If
   obj_func = SS
End Function

Function calc_Cp(tt)
   ' 2-Comp. 点滴モデル、不定用量・不定間隔での濃度算出
   ' tt を採血時刻として引数としている
   Dim i As Integer, t As Single
   Dim E1 As Single, E2 As Single, E3 As Single, E4 As Single, E5 As Single
   Dim AL As Single, BE As Single
      
   flagCheckCalc = 0
   On Error GoTo ErrTrap:
   
   If nVss = 1 Then Vc = P(4)
   If nVss = 2 Then Vc = P(4) / (1 + P(2) / P(3))
   Ke = P(1) / Vc: K12 = P(2): K21 = P(3)
   E1 = Ke + K12 + K21: E2 = Ke * K21
   
   If (E1 * E1 - 4 * E2) <= 0 Then
      Temp2 = MsgBox("血中濃度計算時にエラーが出ました", vbOKOnly + vbExclamation)
      flagCheckCalc = 1: Exit Function
   End If
   
   AL = 0.5 * (E1 + Sqr(E1 * E1 - 4 * E2))
   BE = 0.5 * (E1 - Sqr(E1 * E1 - 4 * E2))
           
   ' set tt
   tn = 0
   For i = 1 To nDose
      If tt >= vv(i, 1) Then
         tn = i
      Else
         t = tt - vv(tn, 1)
      End If
   Next i
   'tt = t
   ' == CHECK ==
   Cells(15, 17) = tt: Cells(16, 17) = t

   ' Cp(tt) の計算
   Cp = 0
   For i = 1 To tn
      If (tt - vv(i, 1)) < vv(i, 3) Then
      ' 点滴中
         E4 = (AL - K21) / AL / (AL - BE) * (1 - Exp(-AL * (tt - vv(i, 1))))
         E5 = (K21 - BE) / BE / (AL - BE) * (1 - Exp(-BE * (tt - vv(i, 1))))
         Cp = Cp + vv(i, 4) / Vc * (E4 + E5)
      Else:
      ' 点滴終了後
         E4 = (AL - K21) / AL / (AL - BE)
         E4 = E4 * (Exp(AL * vv(i, 3)) - 1) * Exp(-AL * (tt - vv(i, 1)))
         E5 = (K21 - BE) / BE / (AL - BE)
         E5 = E5 * (Exp(BE * vv(i, 3)) - 1) * Exp(-BE * (tt - vv(i, 1)))
         Cp = Cp + vv(i, 4) / Vc * (E4 + E5)
      End If
   Next i
   calc_Cp = Cp
   Cells(17, 17) = Cp: ' == CHECK ==
   Exit Function

ErrTrap:
   ' エラーの種類の応じてメッセージ表示
   If Err.Number = 6 Then
      Temp2 = MsgBox("オーバフローしました", vbOKOnly + vbExclamation)
   ElseIf Err.Number = 11 Then
      Temp2 = MsgBox("ゼロで除算しました", vbOKOnly + vbExclamation)
   Else
      Temp2 = MsgBox("血中濃度計算時にエラーがでました", vbOKOnly + vbExclamation)
   End If
   
   flagCheckCalc = 1: Exit Function
End Function

Sub calc_refine()
   Dim E1 As Single, E2 As Single
   
   ' 二次パラメータの計算と表示
   If (nVss = 1) Then
      Pref(3) = P(4) * (1 + P(2) / P(3)):    ' Vss
      Pref(1) = P(1) / P(4):                 ' ke
   End If
   If (nVss = 2) Then
      Pref(3) = P(4) / (1 + P(2) / P(3)):    ' Vc
      Pref(1) = P(1) / Pref(3):              ' ke
   End If
   '
   E1 = Pref(1) + P(2) + P(3): E2 = Pref(1) * P(3)
   If (E1 * E1 - 4 * E2) <= 0 Then
      flagCheckCalc = 1: Exit Sub
   End If
   AL = 0.5 * (E1 + Sqr(E1 * E1 - 4 * E2))
   BE = 0.5 * (E1 - Sqr(E1 * E1 - 4 * E2))
   Pref(2) = Log(2) / BE:                    ' t1/2(b)

End Sub

Sub print_output()
   Dim i As Integer
   
   ' 推定パラメータ値の表示
   For i = 1 To nParam: Cells(52 + i, 11) = P(i): Next i
   For i = 1 To nRef: Cells(52 + i, 13) = Pref(i): Next i
   
   ' 予測値の表示
   For i = 1 To maxrec
      If Cells(67 + i, 6) = 8 Or Cells(67 + i, 6) = 9 Then
         Cells(67 + i, 13) = calc_Cp(Cells(67 + i, 5))
      End If
   Next i
   
   ' 実測値、予測値、残差の表示と確認
   For i = 1 To recN9
      Cells(20 + i, 16) = tx(i):                                 ' time
      Cells(20 + i, 17) = calc_Cp(tx(i))                         ' pred
      Cells(20 + i, 18) = cy(i)                                  ' obs.
      Cells(20 + i, 19) = cy(i) - calc_Cp(tx(i))                 ' res.
   Next i
   For i = 1 To recN8
      Cells(20 + recN9 + 1 + i, 16) = tx8(i):                    ' time
      Cells(20 + recN9 + 1 + i, 17) = calc_Cp(tx8(i)):           ' pred
      Cells(20 + recN9 + 1 + i, 18) = cy8(i):                    ' obs.
      Cells(20 + recN9 + 1 + i, 19) = cy8(i) - calc_Cp(tx8(i)):  ' res.
   Next i
End Sub

Sub plot_output()
   Dim i As Integer, showgrid As Integer, nSiml As Integer
   Dim xmin As Single, xmax As Single
   Dim xx(2000) As Single, yy(2000) As Single, dayBase As Single
   
   ' モデルあてはめ結果のグラフ表示（グラフデータを指定セルに設定）
   Call xAxisChange
   showgrid = Cells(43, 13).Value: nSiml = Int(Cells(44, 11).Value)
   xmin = Cells(45, 11).Value: xmax = Cells(46, 11).Value
   
   dayBase = DateSerial(Year(recDate(1)), Month(recDate(1)), Day(recDate(1)))
   
   For i = 1 To nSiml: yy(i) = 0: Next i
   For i = 1 To nSiml
      xx(i) = xmin + (xmax - xmin) * (i - 1) / (nSiml - 1)
      yy(i) = calc_Cp(xx(i))
      Cells(4 + i, 21).Value = xx(i) / 24 + recTime(1) + dayBase
      Cells(4 + i, 22).Value = yy(i)
   Next i
   
   For i = 1 To recN9
      Cells(4 + i, 23).Value = tx(i) / 24 + recTime(1) + dayBase
      Cells(4 + i, 24).Value = cy(i)
   Next i

   For i = 1 To recN8
      Cells(4 + i, 25).Value = tx8(i) / 24 + recTime(1) + dayBase
      Cells(4 + i, 26).Value = cy8(i)
   Next i
End Sub

Sub plot_output2()
   Dim i As Integer, showgrid As Integer, nSiml As Integer
   Dim xmin As Single, xmax As Single
   Dim xx(2000) As Single, yy(2000) As Single, dayBase As Single
   
   ' 推奨計算した濃度推移のグラフ表示（グラフデータを指定セルに設定）
   showgrid = Cells(43, 13).Value: nSiml = Int(Cells(44, 11).Value)
   xmin = Cells(45, 11).Value
   xmax = Cells(46, 11).Value
   dayBase = DateSerial(Year(recDate(1)), Month(recDate(1)), Day(recDate(1)))
   
   For i = 1 To nSiml: yy(i) = 0: Next i
   For i = 1 To nSiml
      xx(i) = xmin + (xmax - xmin) * (i - 1) / (nSiml - 1)
      yy(i) = calc_Cp(xx(i))
      Cells(4 + i, 27).Value = xx(i) / 24 + recTime(1) + dayBase
      Cells(4 + i, 28).Value = yy(i)
   Next i
End Sub

Sub xAxisChange()
   Dim chart As ChartObject, dayBase As Single
  
   dayBase = DateSerial(Year(recDate(1)), Month(recDate(1)), Day(recDate(1)))
   For Each chart In ActiveSheet.ChartObjects
      chart.Select
      ActiveChart.Axes(xlCategory).Select
      ActiveChart.Axes(xlCategory).MinimumScale = dayBase
   Next chart
End Sub

Sub recom1()
   Dim i As Integer
   
   ' 推奨用法・用量の計算へ移動（事前チェック）
   For i = 1 To 4: Cells(35 + i, 7).Value = "": Next i
   If flagCheckCalc = 1 Then
      Temp2 = MsgBox("パラメータの推定ができていません", vbOKOnly + vbExclamation)
      Exit Sub
   End If
   Cells(34, 3).Select
End Sub

Sub recom2()
   Dim rMethod As Integer, rTau As Single, rTinf As Single
   Dim rAUC As Single, rPeak As Single, rTrough As Single
   Dim rDose As Single, rRate As Single, A1 As Single, B1 As Single, C1 As Single
   Dim i As Integer
   
   ' 推奨用法・用量の計算へ移動
   For i = 1 To 6: Cells(32 + i, 11).Value = "": Next i
   Cells(34, 3).Select: ' 移動
   
   flagCheckRecom = 0
   rMethod = Cells(33, 4): ' 目標設定が AUC か Cpeak / Ctrough か
   If rMethod <> 1 And rMethod <> 2 Then
      Temp2 = MsgBox("Incorrect value in cell(33, 4).", vbOKOnly)
      flagCheckRecom = 1: Exit Sub
   End If
   
   ' 設定条件の読み取り
   rTau = Cells(35, 4): rTinf = Cells(36, 4)
   rAUC = Cells(37, 4): rPeak = Cells(38, 4): rTrough = Cells(39, 4)
   
   ' 推奨値表示を一旦リセット
   Cells(36, 7) = "": Cells(36, 8) = "": ' interval
   Cells(37, 7) = "": Cells(37, 8) = "": ' inf time
   Cells(38, 7) = "": Cells(38, 8) = "": ' dose
   Cells(39, 7) = "": Cells(39, 8) = "": ' rate
   
   If rMethod = 1 Then
   ' AUC を目標値とする ... 投与間隔、点滴時間は固定、点滴速度を計算
      Cells(36, 7).Value = rTau: Cells(36, 8) = rTau
      Cells(37, 7).Value = rTinf: Cells(37, 8) = rTinf
      rDose = rAUC * Vc * Ke * rTau / 24: Cells(38, 7) = rDose
      rRate = rDose / rTinf: Cells(39, 7).Value = rRate
   Else
   ' Peak / trough を目標値とする ... 点滴時間固定、投与間隔、点滴速度を計算
      Cells(37, 7).Value = rTinf
      A1 = (AL - K21) / AL * (1 - Exp(-AL * rTinf))
      B1 = (K21 - BE) / BE * (1 - Exp(-BE * rTinf))
      rTau = (A1 + rPeak / rTrough * B1 * Exp(BE * rTinf)) / (A1 + B1)
      rTau = Log(rTau) / BE
      rDose = rPeak * Exp(BE * rTinf) - rTrough
      rDose = rDose / (A1 + B1) / Exp(BE * rTinf)
      rDose = rDose * rTinf * Vc * (AL - BE)
      Cells(36, 7).Value = rTau: Cells(37, 7).Value = rTinf
      Cells(38, 7).Value = rDose
      rRate = rDose / rTinf: Cells(39, 7).Value = rRate
   End If
End Sub

Sub recom3()
   Dim i As Integer, j As Integer
   Dim rTau As Single, rTinf As Single
   Dim rAUC As Single, rPeak As Single, rTrough As Single
   Dim rDose As Single, rRate As Single, rC1h As Single, rC2h As Single
   Dim A1 As Single, B1 As Single
   
   ' 計算した推奨用量・用法から血中濃度を算出
   ' エラーチェック
   If flagCheckCalc = 1 Then
      Temp2 = MsgBox("パラメータが推定されていません", vbOKOnly)
      Exit Sub
   End If
   If flagCheckRecom = 1 Then
      Temp2 = MsgBox("推奨条件が設定できていません", vbOKOnly)
      flagCheckRecom = 0
      Exit Sub
   End If
     
   rTau = Cells(36, 8): rTinf = Cells(37, 8): rDose = Cells(38, 8): rRate = Cells(39, 8)
   
   If (rTau <= 0 Or rTinf <= 0 Or rDose <= 0 Or rRate <= 0) Then
      Temp2 = MsgBox("推奨条件が不適切です", vbOKOnly)
      flagCheckRecom = 1: Exit Sub
   End If
   
   ' 血中濃度予測値の計算と表示. (AUC, Peak, Trough, Css, C1h, C2h)
   rAUC = rDose / Vc / Ke * 24 / rTau
   ' Peak
   A1 = (AL - K21) / AL * (1 - Exp(-AL * rTinf)) / (1 - Exp(-AL * rTau))
   B1 = (K21 - BE) / BE * (1 - Exp(-BE * rTinf)) / (1 - Exp(-BE * rTau))
   rPeak = rDose / Vc / rTinf / (AL - BE) * (A1 + B1)
   ' Trough
   A1 = (AL - K21) / AL * (1 - Exp(-AL * rTinf))
   A1 = A1 * Exp(-AL * (rTau - rTinf)) / (1 - Exp(-AL * rTau))
   B1 = (K21 - BE) / BE * (1 - Exp(-BE * rTinf))
   B1 = B1 * Exp(-BE * (rTau - rTinf)) / (1 - Exp(-BE * rTau))
   rTrough = rDose / Vc / rTinf / (AL - BE) * (A1 + B1)
   ' C1h
   A1 = (AL - K21) / AL * (1 - Exp(-AL * rTinf)) * Exp(-AL * 1) / (1 - Exp(-AL * rTau))
   B1 = (K21 - BE) / BE * (1 - Exp(-BE * rTinf)) * Exp(-BE * 1) / (1 - Exp(-BE * rTau))
   rC1h = rDose / Vc / rTinf / (AL - BE) * (A1 + B1)
   ' C2h
   A1 = (AL - K21) / AL * (1 - Exp(-AL * rTinf)) * Exp(-AL * 2) / (1 - Exp(-AL * rTau))
   B1 = (K21 - BE) / BE * (1 - Exp(-BE * rTinf)) * Exp(-BE * 2) / (1 - Exp(-BE * rTau))
   rC2h = rDose / Vc / rTinf / (AL - BE) * (A1 + B1)
   
   Cells(33, 11).Value = rAUC
   Cells(34, 11).Value = rPeak: Cells(35, 11).Value = rTrough
   Cells(36, 11).Value = rAUC / 24
   Cells(37, 11).Value = rC1h: Cells(38, 11).Value = rC2h
   
   ' 推奨用法・用量部分の投与記録を一旦リセット
   For i = 1 To 200
      If Cells(67 + i, 13).Value = "*" Then
         maxrec = maxrec - 1
         For j = 3 To 13: Cells(67 + i, j) = "": Next j
      End If
   Next i
   
  ' 推奨用法・用量部分の投与記録を設定
   Cells(67 + maxrec + 1, 3) = Cells(33, 8):  ' day
   Cells(67 + maxrec + 1, 4) = Cells(33, 9):  ' Time
   Cells(67 + maxrec + 1, 6) = 1:             ' EVID
   Cells(67 + maxrec + 1, 7) = 4:             ' Route (inf)
   Cells(67 + maxrec + 1, 8) = Cells(38, 8):  ' Dose
   Cells(67 + maxrec + 1, 9) = Cells(37, 8):  ' Tinf
   Cells(67 + maxrec + 1, 10) = Cells(34, 8): ' nSiml
   Cells(67 + maxrec + 1, 11) = Cells(36, 8): ' Tau
   Cells(67 + maxrec + 1, 13) = "*": ' Check Mark

   flagR = 1
   Call set_Dosage: ' 投与履歴の再読み込み
   flagR = 0
   Call plot_output2: ' 推奨推移のグラフ化
   
   flagCheckRecom = 1
End Sub