Program main
  Use :: nag_library, Only: g02daf, nag_wp
  Implicit None
  Real (Kind=nag_wp) :: mse_tolerance_percent = 1.0d0
  Real (Kind=nag_wp) :: train_data_ratio_percent = 80.0D+0   
  Integer, Parameter :: number_of_folds = 5
  Integer, Parameter :: number_of_alpha_values = 50
  Integer, Parameter :: number_of_ensembles = 5
  Type :: metrics
    Real (Kind=nag_wp) :: mse      ! Mean Squared Error
    Real (Kind=nag_wp) :: mae      ! Mean Absolute Error
    Real (Kind=nag_wp) :: r_squared ! R-squared value
    Real (Kind=nag_wp) :: num_nonzero_coef ! 非ゼロ係数の平均数
  End Type metrics
  Type (metrics) :: model_metrics_lin, model_metrics_lin_org
  Type (metrics) :: model_metrics_lasso, model_metrics_lasso_org
  Type (metrics) :: test_metrics_lin, test_metrics_lin_org
  Type (metrics) :: test_metrics_lasso, test_metrics_lasso_org
  Integer :: n, m, m_plus_one, train_size
  Real (nag_wp), Allocatable :: data(:, :), train_data(:, :), test_data(:, :), means(:), stdvs(:)
  Real (nag_wp), Allocatable :: beta_lin(:), beta_lin_org(:), beta_lasso(:), beta_lasso_org(:)
  Real (nag_wp) :: intercept_lin, intercept_lin_org, intercept_lasso, intercept_lasso_org
  Real (Kind=nag_wp) :: best_alpha
  Character (Len=100) :: filename = 'semi.csv'

  ! データの読み込み
  Call read_data(data, filename)

  ! 最初に全体をシャフルする
  Call shuffle_data(data)

  n = size(data, 1)
  m_plus_one = size(data, 2)
  m = m_plus_one - 1
  train_size = int(n*train_data_ratio_percent/100d0)

  ! Split train_data:test_data
  Allocate (train_data(train_size,m_plus_one), test_data(n-train_size,m_plus_one))
  train_data = data(1:train_size, :)
  test_data = data(train_size+1:n, :)

  Print *, "Performing cross_validaion..."
  Print '(A,F8.3,A)', "   MSE-Tolerance =", mse_tolerance_percent, "%"
  ! cross validation
  Call perform_cross_validation_lasso(train_data, number_of_folds, number_of_alpha_values, number_of_ensembles, mse_tolerance_percent/100d0, &
    best_alpha)
  Print '(A,F8.4)', "Best alpha =", best_alpha

  ! compute means and stdvs
  Allocate (means(m_plus_one), stdvs(m_plus_one))
  Call find_mean_and_std(train_data, means, stdvs)

  ! build final models
  Allocate (beta_lin(m), beta_lasso(m), beta_lin_org(m), beta_lasso_org(m))
  Call standardize(train_data, means, stdvs)
  Call perform_lin_reg(train_data, beta_lin, intercept_lin, model_metrics_lin)
  Call perform_lasso(train_data, best_alpha, beta_lasso, intercept_lasso)
  Call evaluate(train_data, beta_lasso, intercept_lasso, model_metrics_lasso)

  ! desclae model metrics to original scale
  Call descale_metrics(model_metrics_lin, stdvs(m_plus_one), model_metrics_lin_org)
  Call descale_metrics(model_metrics_lasso, stdvs(m_plus_one), model_metrics_lasso_org)

  Print *, "======== Model metrics ========"
  Call print_metrics('  Linear Regression', model_metrics_lin_org)
  Call print_metrics('  Lasso Regression', model_metrics_lasso_org)

 ! モデルの切片と係数を元のスケールに
  Call descale_beta_and_intercept(beta_lin, intercept_lin, stdvs(1:m), stdvs(m_plus_one), means, beta_lin_org, intercept_lin_org)
  Call descale_beta_and_intercept(beta_lasso, intercept_lasso, stdvs(1:m), stdvs(m_plus_one), means, beta_lasso_org, intercept_lasso_org)
  ! モデル係数と切片を表示
  Call print_beta_and_intercept('  Linear Regression', beta_lin_org, intercept_lin_org)
  Call print_beta_and_intercept('  Lasso Regression', beta_lasso_org, intercept_lasso_org)
  
  ! テストデータでの評価
  Call standardize(test_data, means, stdvs)
  Call evaluate(test_data, beta_lin, intercept_lin, test_metrics_lin)
  Call evaluate(test_data, beta_lasso, intercept_lasso, test_metrics_lasso)
  ! 指標をもとのスケールに変換
  Call descale_metrics(test_metrics_lin, stdvs(m_plus_one), test_metrics_lin_org)
  Call descale_metrics(test_metrics_lasso, stdvs(m_plus_one), test_metrics_lasso_org)

  ! 指標を表示
  Print *, "======== Validation Results ========"
  Call print_metrics('Linear Regression', test_metrics_lin_org)
  Call print_metrics('Lasso Regression', test_metrics_lasso_org)
  Print '(A,F8.3)', "    Alpha:", best_alpha

Contains
  Subroutine perform_cross_validation_lasso(data, n_folds, nalpha, n_ensemble, mse_tol, best_alpha_value)
    Use, Intrinsic :: iso_fortran_env, Only: nag_wp => real64
    Implicit None
    Real (Kind=nag_wp), Intent (In) :: data(:, :)
    Integer, Intent (In) :: n_folds, nalpha, n_ensemble
    Real (Kind=nag_wp), Intent (In) :: mse_tol
    Real (Kind=nag_wp), Intent (Out) :: best_alpha_value
    Integer :: i, j, k, n, m, fold_size, start_idx, end_idx
    Real (Kind=nag_wp), Allocatable :: train_data(:, :), val_data(:, :), shuffled_data(:, :)
    Real (Kind=nag_wp), Allocatable :: beta(:), means(:), stdvs(:)
    Real (Kind=nag_wp) :: intercept, alpha
    Type (metrics), Allocatable :: fold_metrics(:, :)
    Real (Kind=nag_wp), Allocatable :: avg_mse_values(:), ensemble_best_alphas(:)
    Real (Kind=nag_wp) :: min_mse, mse_threshold
    Real (Kind=nag_wp) :: alpha_min = 0.0_nag_wp, alpha_max = 1.0_nag_wp

    n = size(data, 1)
    m = size(data, 2)
    fold_size = n/n_folds
    Allocate (fold_metrics(n_folds,nalpha))
    Allocate (beta(m-1), means(m), stdvs(m))
    Allocate (avg_mse_values(nalpha))
    Allocate (ensemble_best_alphas(n_ensemble))
    Allocate (shuffled_data(n,m))
    Do k = 1, n_ensemble
      Do i = 1, n_folds
        Call make_shuffled_copy(data, shuffled_data)
        start_idx = (i-1)*fold_size + 1
        end_idx = i*fold_size
        Allocate (train_data(n-fold_size,m), val_data(fold_size,m))
        train_data(1:start_idx-1, :) = shuffled_data(1:start_idx-1, :)
        train_data(start_idx:n-fold_size, :) = shuffled_data(end_idx+1:n, :)
        val_data = shuffled_data(start_idx:end_idx, :)
        Call find_mean_and_std(train_data, means, stdvs)
        Call standardize(train_data, means, stdvs)
        Call standardize(val_data, means, stdvs)
        Do j = 1, nalpha
          alpha = alpha_min + (alpha_max-alpha_min)*(j-1)/(nalpha-1)
          Call perform_lasso(train_data, alpha, beta, intercept)
          Call evaluate(val_data, beta, intercept, fold_metrics(i,j))
        End Do
        Deallocate (train_data, val_data)
      End Do
      Do j = 1, nalpha
        avg_mse_values(j) = sum(fold_metrics(:,j)%mse)/n_folds
      End Do
      min_mse = minval(avg_mse_values)
      mse_threshold = min_mse*(1.0_nag_wp+mse_tol)
      Do j = nalpha, 1, -1
        If (avg_mse_values(j)<=mse_threshold) Then
          ensemble_best_alphas(k) = alpha_min + (alpha_max-alpha_min)*(j-1)/(nalpha-1)
          Exit
        End If
      End Do
    End Do
    best_alpha_value = sum(ensemble_best_alphas)/n_ensemble
    Deallocate (fold_metrics, beta, means, stdvs, avg_mse_values, ensemble_best_alphas, shuffled_data)
  End Subroutine perform_cross_validation_lasso


  Subroutine calculate_std_metrics(fold_metrics, avg_metrics, std_metrics)
    Type (metrics), Intent (In) :: fold_metrics(:)
    Type (metrics), Intent (In) :: avg_metrics
    Type (metrics), Intent (Out) :: std_metrics
    Integer :: n

    n = size(fold_metrics)
    std_metrics%mse = sqrt(sum((fold_metrics%mse-avg_metrics%mse)**2)/(n-1))
    std_metrics%mae = sqrt(sum((fold_metrics%mae-avg_metrics%mae)**2)/(n-1))
    std_metrics%r_squared = sqrt(sum((fold_metrics%r_squared-avg_metrics%r_squared)**2)/(n-1))
    std_metrics%num_nonzero_coef = sqrt(sum((fold_metrics%num_nonzero_coef-avg_metrics%num_nonzero_coef)**2)/(n-1))
  End Subroutine calculate_std_metrics


  Subroutine perform_cross_validation_linear(data, n_folds, avg_metrics, std_metrics)
    Real (Kind=nag_wp), Intent (In) :: data(:, :)
    Integer, Intent (In) :: n_folds
    Type (metrics), Intent (Out) :: avg_metrics, std_metrics

    Integer :: i, n, m, fold_size, start_idx, end_idx
    Real (Kind=nag_wp), Allocatable :: train_data(:, :), val_data(:, :)
    Real (Kind=nag_wp), Allocatable :: beta(:), means(:), stdvs(:)
    Real (Kind=nag_wp) :: intercept
    Type (metrics), Allocatable :: fold_metrics(:)

    n = size(data, 1)
    m = size(data, 2)
    fold_size = n/n_folds

    Allocate (fold_metrics(n_folds))
    Allocate (beta(m-1), means(m), stdvs(m))

    Do i = 1, n_folds
      start_idx = (i-1)*fold_size + 1
      end_idx = i*fold_size

      ! Split data into training and validation sets
      Allocate (train_data(n-fold_size,m), val_data(fold_size,m))

      ! Correct way to combine arrays
      train_data(1:start_idx-1, :) = data(1:start_idx-1, :)
      train_data(start_idx:n-fold_size, :) = data(end_idx+1:n, :)

      val_data = data(start_idx:end_idx, :)

      ! Standardize data
      Call find_mean_and_std(train_data, means, stdvs)
      Call standardize(train_data, means, stdvs)
      Call standardize(val_data, means, stdvs)

      ! Perform linear regression
      Call perform_lin_reg(train_data, beta, intercept, fold_metrics(i))

      ! Evaluate on validation set
      Call evaluate(val_data, beta, intercept, fold_metrics(i))

      Deallocate (train_data, val_data)
    End Do

    ! Calculate average and standard deviation of metrics
    Call calculate_avg_std_metrics(fold_metrics, avg_metrics, std_metrics)

    Deallocate (fold_metrics, beta, means, stdvs)
  End Subroutine perform_cross_validation_linear


  Subroutine calculate_avg_std_metrics(fold_metrics, avg_metrics, std_metrics)
    Type (metrics), Intent (In) :: fold_metrics(:)
    Type (metrics), Intent (Out) :: avg_metrics, std_metrics
    Integer :: n

    n = size(fold_metrics)

    ! Calculate averages
    avg_metrics%mse = sum( [fold_metrics%mse] )/n
    avg_metrics%mae = sum( [fold_metrics%mae] )/n
    avg_metrics%r_squared = sum( [fold_metrics%r_squared] )/n
    avg_metrics%num_nonzero_coef = sum( [fold_metrics%num_nonzero_coef] )/n

    ! Calculate standard deviations
    std_metrics%mse = sqrt(sum(([fold_metrics%mse]-avg_metrics%mse)**2)/(n-1))
    std_metrics%mae = sqrt(sum(([fold_metrics%mae]-avg_metrics%mae)**2)/(n-1))
    std_metrics%r_squared = sqrt(sum(([fold_metrics%r_squared]-avg_metrics%r_squared)**2)/(n-1))
    std_metrics%num_nonzero_coef = sqrt(sum(([fold_metrics%num_nonzero_coef]-avg_metrics%num_nonzero_coef)**2)/(n-1))

  End Subroutine calculate_avg_std_metrics

  Subroutine print_beta_and_intercept(title, beta, intercept)
    Character (*), Intent (In) :: title
    Real (Kind=nag_wp), Intent (In) :: beta(:)
    Real (Kind=nag_wp), Intent (In) :: intercept
    Integer :: i

    Print *, 'Beta and Intercept for: ', trim(title)
    Print *, '-----------------------------------------------'
    Do i = 1, size(beta)
      Print '(A,I2,A,F12.5)', '    beta(', i, ') =', beta(i)
    End Do
    Print '(A,F12.5)', '    Intercept =', intercept
  End Subroutine print_beta_and_intercept

  Subroutine print_metrics(title, avg_met, std_met)
    Character (*), Intent (In) :: title
    Type (metrics), Intent (In) :: avg_met
    Type (metrics), Intent (In), Optional :: std_met

    If (present(std_met)) Then
      Print *, 'Metrics for: ', trim(title)
      Print '(A,F16.8,A,F16.8)', '   Average R2:  ', avg_met%r_squared, ' Std Dev: ', std_met%r_squared
      Print '(A,F16.8,A,F16.8)', '   Average MSE: ', avg_met%mse, ' Std Dev: ', std_met%mse
      Print '(A,F16.8,A,F16.8)', '   Average MAE: ', avg_met%mae, ' Std Dev: ', std_met%mae
      Print '(A,F10.2)', '   Average Non-zero Coefficients: ', avg_met%num_nonzero_coef
    Else
      Print *, 'Metrics for: ', trim(title)
      Print '(A,F16.8)', '   Average R2:  ', avg_met%r_squared
      Print '(A,F16.8)', '   Average MSE: ', avg_met%mse
      Print '(A,F16.8)', '   Average MAE: ', avg_met%mae
      Print '(A,F10.2)', '   Average Non-zero Coefficients: ', avg_met%num_nonzero_coef
    End If
  End Subroutine print_metrics


  Subroutine perform_lin_reg(data, beta, intercept, met)
    Real (nag_wp), Intent (In) :: data(:, :)
    Real (nag_wp), Intent (Out) :: beta(:)
    Type (metrics), Intent (Out) :: met
    Real (nag_wp) :: rss, intercept
    Integer :: idf
    Real (nag_wp), Allocatable :: x(:, :), y(:), wt(:), b(:), se(:), cov(:), h(:), q(:, :), p(:), wk(:), res(:)
    Integer, Allocatable :: isx(:)
    Real (nag_wp) :: tol
    Integer :: n, m, ldx, ldq, ifail, ip, irank
    Logical :: svd
    Character (1) :: mean, weight

    n = size(data, 1)
    m = size(data, 2) - 1          ! 一番右の列はY
    ldx = n
    ldq = n

    Allocate (x(ldx,m), y(n), isx(m), b(m+1), se(m+1), cov((m+2)*(m+1)/2), h(n))
    Allocate (q(ldq,m+2), p(2*(m+1)+(m+1)*(m+1)), wk(max(2,5*m+m*m)), res(n))
    Allocate (wt(1))               ! 重みなしの場合

    x = data(:, 1:m)
    y = data(:, size(data,2))      ! 一番右の列がY

    mean = 'M'                     ! 平均項を含める
    weight = 'U'                   ! 重みなし
    ip = m + 1                     ! 定数項を含む独立変数の数
    isx = 1                        ! すべての変数を使用
    tol = 0.000001_nag_wp          ! 許容誤差

    ifail = -1
    Call g02daf(mean, weight, n, x, ldx, m, isx, ip, y, wt, rss, idf, b, se, cov, res, h, q, ldq, svd, irank, p, tol, wk, ifail)

    If (ifail==0) Then
      If (svd) Then
        Print *, 'Warning in g02daf.  Model is not full rank', ip, m
      End If
      beta = b(2:)
      intercept = b(1)
      ! print '(A,(99(F7.3)))', "B and Intercept ", b
    Else
      Print *, 'Error in g02daf. ifail =', ifail
    End If

    Call calculate_metrics(rss, res, y, n, idf, met%r_squared, met%mse, met%mae)
    met%num_nonzero_coef = count(abs(beta)>1.0E-6_nag_wp)

    Deallocate (x, y, isx, b, se, cov, h, q, p, wk, wt, res)
  End Subroutine perform_lin_reg

  Subroutine perform_lasso(data, alpha, best_beta, best_intercept)
    Real (Kind=nag_wp), Intent (In) :: data(:, :)
    Real (Kind=nag_wp), Intent (In) :: alpha ! 0 <= alpha <= 1
    Real (Kind=nag_wp), Allocatable, Intent (Out) :: best_beta(:)
    Real (Kind=nag_wp), Intent (Out) :: best_intercept
    Integer :: n, m, mtype, pred, prey, mnstep, lisx, ip, nstep, ldb, ldd, lropt, ifail
    Real (Kind=nag_wp), Allocatable :: b(:, :), fitsum(:, :), ropt(:), y(:), x(:, :)
    Integer, Allocatable :: isx(:)
    Integer :: ktype, lnk, ldnb
    Real (Kind=nag_wp), Allocatable :: nk(:), nb(:, :)
    Real (Kind=nag_wp) :: alpha_ranged

    alpha_ranged = alpha

    ! 入力パラメータのチェック
    If (alpha_ranged<0.0_nag_wp) alpha_ranged = 0.0D+0
    If (alpha_ranged>1.0_nag_wp) alpha_ranged = 1.0D+0

    ! データのサイズ取得
    n = size(data, 1)
    m = size(data, 2) - 1          ! 最後の1列は目的変数

    ! モデルの設定
    mtype = 3                      ! LASSO
    pred = 0                       ! 説明変数の標準化を行わない（すでに標準化済み）
    prey = 0                       ! 目的変数の平均を考慮する
    mnstep = m*10                  ! 最大ステップ数を設定
    lisx = 0                       ! すべての変数を含める

    ! オプション引数（デフォルト値を使用）
    lropt = 5
    Allocate (ropt(lropt))
    ropt = 0.0_nag_wp
    ropt(1) = 1.0E-7_nag_wp        ! 最小ステップサイズ
    ropt(2) = 1.0E-7_nag_wp        ! 一般的な許容誤差
    ropt(3) = 1.0_nag_wp           ! パラメータ推定値をリスケールする
    ropt(4) = 1.0_nag_wp           ! モデルに切片を含める
    ropt(5) = 0.0_nag_wp           ! クロスプロダクト行列を使用

    ! 出力配列の割り当て
    ip = m
    ldb = ip
    Allocate (b(ldb,mnstep+2), fitsum(6,mnstep+1))
    Allocate (isx(0))              ! to avoid error

    ! 説明変数と目的変数の分離
    Allocate (x(n,m), y(n))
    x = data(:, 1:m)
    y = data(:, m+1)
    ldd = n

    ! モデル適合ルーチンの呼び出し
    ifail = 0
    Call g02maf(mtype, pred, prey, n, m, x, ldd, isx, lisx, y, mnstep, ip, nstep, b, ldb, fitsum, ropt, lropt, ifail)

    ! エラーチェック
    If (ifail/=0) Then
      Print *, 'Error in g02maf. ifail =', ifail
      Return
    End If

    ! nkとktypeの設定
    ktype = 3                      ! スケーリングされたL1ノルムの最大値に対する比率
    lnk = 1
    Allocate (nk(lnk))
    nk(1) = 1.0D+0 - alpha_ranged

    ! nbの割り当て
    ldnb = ip
    Allocate (nb(ldnb,lnk))

    ! g02mcfの呼び出し
    ifail = 0
    Call g02mcf(nstep, ip, b, ldb, fitsum, ktype, nk, lnk, nb, ldnb, ifail)

    ! エラーチェック
    If (ifail/=0) Then
      Print *, 'Error in g02mcf. ifail =', ifail
      Return
    End If

    ! パラメータ推定値の取得
    best_beta = nb(:, 1)
    best_intercept = fitsum(1, nstep+1) ! 切片はfitsum(1, nstep+1)に格納されている

    ! メモリの解放
    Deallocate (b, fitsum, x, y, nb, nk, ropt)

  End Subroutine perform_lasso




  Subroutine calculate_metrics(rss, residuals, y, n, df, rsq, mse, mae)
    Real (Kind=nag_wp), Intent (In) :: rss
    Real (Kind=nag_wp), Intent (In) :: residuals(:)
    Real (Kind=nag_wp), Intent (In) :: y(:)
    Integer, Intent (In) :: n, df
    Real (Kind=nag_wp), Intent (Out) :: rsq, mse, mae
    Real (Kind=nag_wp) :: tss, y_mean

    ! 総平方和 (TSS) の計算
    y_mean = sum(y)/n
    tss = sum((y-y_mean)**2)

    ! 決定係数 (R-squared) の計算
    If (tss>0.0_nag_wp) Then
      rsq = 1.0_nag_wp - rss/tss
    Else
      rsq = 0.0_nag_wp
    End If

    ! 平均二乗誤差 (MSE) の計算
    mse = rss/df

    ! 平均絶対誤差 (MAE) の計算
    mae = sum(abs(residuals))/n

  End Subroutine calculate_metrics

  Subroutine read_data(data, filename)
    Real (Kind=nag_wp), Allocatable, Intent (Out) :: data(:, :)
    Character (Len=*), Intent (In) :: filename
    Integer :: io_status, num_rows, num_cols, i
    Character (Len=4000) :: line

    ! カウント行数とカラム数
    num_rows = 0
    Open (Unit=10, File=filename, Status='old', Action='read', Iostat=io_status)
    If (io_status/=0) Then
      Print *, 'Error opening file: ', filename
      Return
    End If

    ! 1行目を読み込んでカラム数を数える
    Read (10, '(A)', Iostat=io_status) line
    num_cols = count( [(line(i:i),i=1,len_trim(line))] ==',') + 1
    num_rows = num_rows + 1
    ! num_cols = count(transfer(header_line,'a',len(header_line))==',') + 1

    ! 行数を数える
    Do
      Read (10, '(A)', Iostat=io_status) line
      If (io_status/=0) Exit
      If (len_trim(line)==0) Exit
      num_rows = num_rows + 1
    End Do

    ! データの読み込み
    Allocate (data(num_rows,num_cols))
    Rewind (10)
    Do i = 1, num_rows
      Read (10, *, Iostat=io_status) data(i, :)
    End Do

    Close (10)
  End Subroutine read_data


  Subroutine find_mean_and_std(data, means, stdvs)
    Real (nag_wp), Intent (In) :: data(:, :)
    Real (nag_wp), Intent (Out) :: means(:), stdvs(:)
    Integer :: n

    n = size(data, 1)
    means = sum(data, dim=1)/n
    stdvs = sqrt(sum((data-spread(means,1,n))**2,dim=1)/(n-1))
  End Subroutine find_mean_and_std

  Subroutine standardize(data, means, stdvs)
    Real (nag_wp), Intent (Inout) :: data(:, :)
    Real (nag_wp), Intent (In) :: means(:), stdvs(:)
    Integer :: n

    n = size(data, 1)
    data = (data-spread(means,1,n))/spread(stdvs, 1, n)
  End Subroutine standardize


  Subroutine descale_beta_and_intercept(beta_lin, intercept_lin, x_std, y_std, means, beta_lin_org, intercept_lin_org)
    Real (Kind=nag_wp), Intent (In) :: beta_lin(:), intercept_lin
    Real (Kind=nag_wp), Intent (In) :: x_std(:), y_std
    Real (Kind=nag_wp), Intent (In) :: means(:)
    Real (Kind=nag_wp), Intent (Out) :: beta_lin_org(:)
    Real (Kind=nag_wp), Intent (Out) :: intercept_lin_org
    Integer :: i, m

    m = size(beta_lin)

    ! Descale beta coefficients
    Do i = 1, m
      beta_lin_org(i) = beta_lin(i)*(y_std/x_std(i))
    End Do

    ! Descale intercept
    intercept_lin_org = means(m+1) + y_std*intercept_lin
    Do i = 1, m
      intercept_lin_org = intercept_lin_org - beta_lin_org(i)*means(i)
    End Do
  End Subroutine descale_beta_and_intercept


  Subroutine descale_metrics(met_in, y_std, met_out)
    Type (metrics), Intent (In) :: met_in
    Real (Kind=nag_wp), Intent (In) :: y_std
    Type (metrics), Intent (Out) :: met_out

    met_out = met_in
    met_out%mse = met_in%mse*(y_std**2)
    met_out%mae = met_in%mae*y_std
  End Subroutine descale_metrics

  Subroutine evaluate(data_in, beta_in, intercept_in, met_out)
    Use :: nag_library, Only: nag_wp
    Implicit None
    Real (Kind=nag_wp), Intent (In) :: data_in(:, :)
    Real (Kind=nag_wp), Intent (In) :: beta_in(:)
    Real (Kind=nag_wp), Intent (In) :: intercept_in
    Type (metrics), Intent (Out) :: met_out
    Integer :: n, m, i, df
    Real (Kind=nag_wp), Allocatable :: y_pred(:), residuals(:)
    Real (Kind=nag_wp) :: rss

    n = size(data_in, 1)
    m = size(data_in, 2) - 1       ! 最後の列は応答変数
    df = n - m - 1                 ! 自由度

    Allocate (y_pred(n), residuals(n))

    ! 予測値とバラつきの計算
    Do i = 1, n
      y_pred(i) = intercept_in + dot_product(beta_in, data_in(i,1:m))
      residuals(i) = data_in(i, m+1) - y_pred(i)
    End Do

    ! RSSの計算
    rss = sum(residuals**2)


    ! メトリクスの計算
    Call calculate_metrics(rss, residuals, data_in(:,m+1), n, df, met_out%r_squared, met_out%mse, met_out%mae)

    ! 非ゼロ係数の数を計算
    met_out%num_nonzero_coef = count(abs(beta_in)>1.0E-6_nag_wp)

    Deallocate (y_pred, residuals)
  End Subroutine evaluate


  Subroutine make_shuffled_copy(indata, outdata)
    Real (Kind=nag_wp), Intent (In) :: indata(:, :)
    Real (Kind=nag_wp), Intent (Out) :: outdata(:, :)
    Integer :: i, n, idx
    Real (Kind=nag_wp) :: temp(size(indata,2)), r

    n = size(indata, 1)
    outdata = indata
    Do i = n, 2, -1
      Call random_number(r)
      idx = int(r*i) + 1
      temp = outdata(i, :)
      outdata(i, :) = outdata(idx, :)
      outdata(idx, :) = temp
    End Do
  End Subroutine make_shuffled_copy

  Subroutine shuffle_data(data)
    Real (Kind=nag_wp), Intent (Inout) :: data(:, :)
    Real (Kind=nag_wp), Allocatable :: tmp_data(:, :)

    Allocate (tmp_data(size(data,1),size(data,2)))
    Call make_shuffled_copy(data, tmp_data)
    data = tmp_data
    Deallocate (tmp_data)
  End Subroutine shuffle_data

End Program main