!-----------------------------------------------------------------------
! シグモイド偏微分方程式近似解法プログラム
! 
! 超次元空間における入力ベクトルと出力ベクトル間の
! シグモイド偏微分方程式の数値近似解法を実装します
!-----------------------------------------------------------------------

MODULE sigmoid_pde_module
  IMPLICIT NONE
  
  ! パラメータ設定
  INTEGER, PARAMETER :: DP = SELECTED_REAL_KIND(15, 307)
  REAL(DP), PARAMETER :: EPSILON = 1.0E-6_DP
  REAL(DP), PARAMETER :: ALPHA = 0.5_DP    ! 学習率
  INTEGER, PARAMETER :: MAX_ITER = 1000    ! 最大反復回数
  
CONTAINS

  !-----------------------------------------------------------------------
  ! シグモイド関数
  !-----------------------------------------------------------------------
  FUNCTION sigmoid(x)
    REAL(DP), INTENT(IN) :: x
    REAL(DP) :: sigmoid
    
    sigmoid = 1.0_DP / (1.0_DP + EXP(-x))
  END FUNCTION sigmoid
  
  !-----------------------------------------------------------------------
  ! シグモイド関数の導関数
  !-----------------------------------------------------------------------
  FUNCTION sigmoid_derivative(x)
    REAL(DP), INTENT(IN) :: x
    REAL(DP) :: sigmoid_derivative
    REAL(DP) :: sig
    
    sig = sigmoid(x)
    sigmoid_derivative = sig * (1.0_DP - sig)
  END FUNCTION sigmoid_derivative
  
  !-----------------------------------------------------------------------
  ! 超次元空間における勾配計算
  !-----------------------------------------------------------------------
  SUBROUTINE compute_gradient(input_vec, output_vec, target_vec, grad, dim_in, dim_out)
    INTEGER, INTENT(IN) :: dim_in, dim_out
    REAL(DP), INTENT(IN) :: input_vec(dim_in)
    REAL(DP), INTENT(IN) :: output_vec(dim_out)
    REAL(DP), INTENT(IN) :: target_vec(dim_out)
    REAL(DP), INTENT(OUT) :: grad(dim_in, dim_out)
    
    INTEGER :: i, j
    REAL(DP) :: error(dim_out)
    REAL(DP) :: deriv
    
    ! 誤差計算
    error = target_vec - output_vec
    
    ! 超次元空間での勾配計算
    DO j = 1, dim_out
      DO i = 1, dim_in
        deriv = sigmoid_derivative(output_vec(j))
        grad(i, j) = error(j) * deriv * input_vec(i)
      END DO
    END DO
  END SUBROUTINE compute_gradient
  
  !-----------------------------------------------------------------------
  ! 偏微分方程式のヤコビアン行列計算
  !-----------------------------------------------------------------------
  SUBROUTINE compute_jacobian(weight_matrix, input_vec, jacobian, dim_in, dim_out)
    INTEGER, INTENT(IN) :: dim_in, dim_out
    REAL(DP), INTENT(IN) :: weight_matrix(dim_in, dim_out)
    REAL(DP), INTENT(IN) :: input_vec(dim_in)
    REAL(DP), INTENT(OUT) :: jacobian(dim_out, dim_out)
    
    INTEGER :: i, j, k
    REAL(DP) :: temp
    
    jacobian = 0.0_DP
    
    ! ヤコビアン行列の計算
    DO i = 1, dim_out
      DO j = 1, dim_out
        temp = 0.0_DP
        DO k = 1, dim_in
          temp = temp + weight_matrix(k,i) * weight_matrix(k,j) * input_vec(k)
        END DO
        jacobian(i,j) = temp
      END DO
    END DO
  END SUBROUTINE compute_jacobian
  
  !-----------------------------------------------------------------------
  ! 前進オイラー法による数値積分
  !-----------------------------------------------------------------------
  SUBROUTINE forward_euler_step(weight_matrix, grad, dim_in, dim_out, dt)
    INTEGER, INTENT(IN) :: dim_in, dim_out
    REAL(DP), INTENT(INOUT) :: weight_matrix(dim_in, dim_out)
    REAL(DP), INTENT(IN) :: grad(dim_in, dim_out)
    REAL(DP), INTENT(IN) :: dt
    
    weight_matrix = weight_matrix + dt * grad
  END SUBROUTINE forward_euler_step
  
  !-----------------------------------------------------------------------
  ! ルンゲ・クッタ法による数値積分（高精度）
  !-----------------------------------------------------------------------
  SUBROUTINE runge_kutta_step(weight_matrix, input_vec, target_vec, dim_in, dim_out, dt)
    INTEGER, INTENT(IN) :: dim_in, dim_out
    REAL(DP), INTENT(INOUT) :: weight_matrix(dim_in, dim_out)
    REAL(DP), INTENT(IN) :: input_vec(dim_in)
    REAL(DP), INTENT(IN) :: target_vec(dim_out)
    REAL(DP), INTENT(IN) :: dt
    
    REAL(DP) :: k1(dim_in, dim_out), k2(dim_in, dim_out)
    REAL(DP) :: k3(dim_in, dim_out), k4(dim_in, dim_out)
    REAL(DP) :: temp_weights(dim_in, dim_out)
    REAL(DP) :: output_vec(dim_out)
    
    ! K1 計算
    CALL forward_propagation(weight_matrix, input_vec, output_vec, dim_in, dim_out)
    CALL compute_gradient(input_vec, output_vec, target_vec, k1, dim_in, dim_out)
    
    ! K2 計算
    temp_weights = weight_matrix + 0.5_DP * dt * k1
    CALL forward_propagation(temp_weights, input_vec, output_vec, dim_in, dim_out)
    CALL compute_gradient(input_vec, output_vec, target_vec, k2, dim_in, dim_out)
    
    ! K3 計算
    temp_weights = weight_matrix + 0.5_DP * dt * k2
    CALL forward_propagation(temp_weights, input_vec, output_vec, dim_in, dim_out)
    CALL compute_gradient(input_vec, output_vec, target_vec, k3, dim_in, dim_out)
    
    ! K4 計算
    temp_weights = weight_matrix + dt * k3
    CALL forward_propagation(temp_weights, input_vec, output_vec, dim_in, dim_out)
    CALL compute_gradient(input_vec, output_vec, target_vec, k4, dim_in, dim_out)
    
    ! 最終更新
    weight_matrix = weight_matrix + dt * (k1 + 2.0_DP*k2 + 2.0_DP*k3 + k4) / 6.0_DP
  END SUBROUTINE runge_kutta_step
  
  !-----------------------------------------------------------------------
  ! 順伝播計算
  !-----------------------------------------------------------------------
  SUBROUTINE forward_propagation(weight_matrix, input_vec, output_vec, dim_in, dim_out)
    INTEGER, INTENT(IN) :: dim_in, dim_out
    REAL(DP), INTENT(IN) :: weight_matrix(dim_in, dim_out)
    REAL(DP), INTENT(IN) :: input_vec(dim_in)
    REAL(DP), INTENT(OUT) :: output_vec(dim_out)
    
    INTEGER :: i, j
    REAL(DP) :: activation
    
    DO j = 1, dim_out
      activation = 0.0_DP
      DO i = 1, dim_in
        activation = activation + weight_matrix(i,j) * input_vec(i)
      END DO
      output_vec(j) = sigmoid(activation)
    END DO
  END SUBROUTINE forward_propagation
  
  !-----------------------------------------------------------------------
  ! 平均二乗誤差計算
  !-----------------------------------------------------------------------
  FUNCTION mean_squared_error(output_vec, target_vec, dim_out)
    INTEGER, INTENT(IN) :: dim_out
    REAL(DP), INTENT(IN) :: output_vec(dim_out)
    REAL(DP), INTENT(IN) :: target_vec(dim_out)
    REAL(DP) :: mean_squared_error
    
    INTEGER :: i
    REAL(DP) :: sum_error
    
    sum_error = 0.0_DP
    DO i = 1, dim_out
      sum_error = sum_error + (output_vec(i) - target_vec(i))**2
    END DO
    
    mean_squared_error = sum_error / REAL(dim_out, DP)
  END FUNCTION mean_squared_error
  
  !-----------------------------------------------------------------------
  ! シグモイド偏微分方程式の近似解法
  !-----------------------------------------------------------------------
  SUBROUTINE solve_sigmoid_pde(input_data, target_data, weight_matrix, &
                               n_samples, dim_in, dim_out)
    INTEGER, INTENT(IN) :: n_samples, dim_in, dim_out
    REAL(DP), INTENT(IN) :: input_data(n_samples, dim_in)
    REAL(DP), INTENT(IN) :: target_data(n_samples, dim_out)
    REAL(DP), INTENT(OUT) :: weight_matrix(dim_in, dim_out)
    
    INTEGER :: iter, sample_idx
    REAL(DP) :: output_vec(dim_out)
    REAL(DP) :: grad(dim_in, dim_out)
    REAL(DP) :: error, prev_error, dt
    REAL(DP) :: input_vec(dim_in)
    REAL(DP) :: target_vec(dim_out)
    
    ! 重み行列の初期化 (小さなランダム値)
    CALL random_seed()
    CALL random_number(weight_matrix)
    weight_matrix = 0.1_DP * (weight_matrix - 0.5_DP)
    
    prev_error = HUGE(0.0_DP)
    dt = ALPHA  ! 時間ステップ = 学習率
    
    ! 反復計算
    DO iter = 1, MAX_ITER
      error = 0.0_DP
      
      ! 各サンプルに対して処理
      DO sample_idx = 1, n_samples
        input_vec = input_data(sample_idx, :)
        target_vec = target_data(sample_idx, :)
        
        ! 順伝播
        CALL forward_propagation(weight_matrix, input_vec, output_vec, dim_in, dim_out)
        
        ! 誤差計算
        error = error + mean_squared_error(output_vec, target_vec, dim_out)
        
        ! 勾配計算
        CALL compute_gradient(input_vec, output_vec, target_vec, grad, dim_in, dim_out)
        
        ! ルンゲ・クッタ法による重み更新（より高精度な数値解法）
        CALL runge_kutta_step(weight_matrix, input_vec, target_vec, dim_in, dim_out, dt)
      END DO
      
      ! 平均誤差
      error = error / REAL(n_samples, DP)
      
      ! 収束判定
      IF (ABS(error - prev_error) < EPSILON) THEN
        PRINT *, "収束しました。反復回数: ", iter
        EXIT
      END IF
      
      IF (MOD(iter, 100) == 0) THEN
        PRINT *, "反復: ", iter, " 誤差: ", error
      END IF
      
      prev_error = error
    END DO
    
    IF (iter > MAX_ITER) THEN
      PRINT *, "最大反復回数に達しました。最終誤差: ", error
    END IF
  END SUBROUTINE solve_sigmoid_pde

END MODULE sigmoid_pde_module

!-----------------------------------------------------------------------
! メインプログラム
!-----------------------------------------------------------------------
PROGRAM main
  USE sigmoid_pde_module
  IMPLICIT NONE
  
  INTEGER, PARAMETER :: DP = SELECTED_REAL_KIND(15, 307)
  INTEGER :: dim_in, dim_out, n_samples, i, j
  REAL(DP), ALLOCATABLE :: input_data(:,:), target_data(:,:), weight_matrix(:,:)
  REAL(DP), ALLOCATABLE :: test_input(:), test_output(:)
  
  ! 問題の次元とサンプル数の設定
  dim_in = 5      ! 入力ベクトルの次元
  dim_out = 3     ! 出力ベクトルの次元
  n_samples = 100 ! 訓練サンプル数
  
  ! メモリ割り当て
  ALLOCATE(input_data(n_samples, dim_in))
  ALLOCATE(target_data(n_samples, dim_out))
  ALLOCATE(weight_matrix(dim_in, dim_out))
  ALLOCATE(test_input(dim_in))
  ALLOCATE(test_output(dim_out))
  
  ! サンプルデータの生成（実際のアプリケーションではデータ読み込みが必要）
  CALL random_seed()
  CALL random_number(input_data)
  
  ! ターゲットデータの生成（単純化のためのダミーデータ）
  DO i = 1, n_samples
    DO j = 1, dim_out
      ! 単純な非線形関係を仮定
      target_data(i, j) = 0.5_DP * SIN(SUM(input_data(i, :)) * j) + 0.5_DP
    END DO
  END DO
  
  ! シグモイド偏微分方程式の近似解法の実行
  PRINT *, "超次元空間におけるシグモイド偏微分方程式の近似解法を開始します..."
  CALL solve_sigmoid_pde(input_data, target_data, weight_matrix, n_samples, dim_in, dim_out)
  
  ! 結果のテスト
  PRINT *, "訓練後の重み行列:"
  DO i = 1, dim_in
    WRITE(*, '(*(F8.4, 1X))') (weight_matrix(i, j), j = 1, dim_out)
  END DO
  
  ! テストデータでの予測
  PRINT *, "テストデータでの予測:"
  CALL random_number(test_input)
  CALL forward_propagation(weight_matrix, test_input, test_output, dim_in, dim_out)
  
  PRINT *, "入力ベクトル:"
  WRITE(*, '(*(F8.4, 1X))') test_input
  
  PRINT *, "予測出力ベクトル:"
  WRITE(*, '(*(F8.4, 1X))') test_output
  
  ! メモリ解放
  DEALLOCATE(input_data, target_data, weight_matrix, test_input, test_output)
  
  PRINT *, "プログラム終了"
END PROGRAM main