% transform_to_corr_mat: whether to trasnform the input and output covariance matrix
% to the correlation matrix.
%   0: no transform
%   1: transform (default)
% Input data can be either a covariance matrix or a correlation matrix

clear
curr_dir = '../data'
% curr_dir = '.' % should be modified to be consistent with your folder name

% We provide only motivation.txt (as described in Masuda, Kojaku & Sano, Physical Review E, 2018)
% Matrix C can be either a covariance or correlation matrix
C = load([curr_dir '/motivation.txt']);
% C = load([curr_dir '/cov105115_all.txt']);
% C = load([curr_dir '/cov118932_all.txt']);
% C = load([curr_dir '/JapanCov.txt']);
% C = load([curr_dir '/USCov.txt']);

N = size(C,1); % size of the matrix


min_eig_C = min(eig(C));
if min_eig_C < 1e-6 * max(diag(C))
    fprintf('min eig = %f\n', min_eig_C);
    error('Input correlation/covariance matrix must be of full rank');
end

tolerance = 1e-5; % to judge whether the algorithm has converged. In the paper = 1e-5

r = 1e-4; % learning rate. If r is too large, the algorithm would not converge
% default r = 1e-4

transform_to_corr_mat = 1; 
% If transform_to_corr_mat == 1, transform the input covariance matrix to the correlation matrix before running the gradient descent method
% Default value is 1
[C_con, alpha, beta, it] = max_ent_config_naive_gradient_descent(C, tolerance, r, transform_to_corr_mat);
% C_con is a covariance matrix but not a correlation matrix

if transform_to_corr_mat == 1
    C_con = diag(diag(C_con))^(-1/2) * C_con * diag(diag(C_con))^(-1/2); % corr matrix
end
% C_con is either a covariance or correlation matrix depending on the value of transform_to_corr_mat.

K_con = inv(C_con); % precision matrix

L = N; % L should be set to the length of the original data for which the Pearson correlation is calculated (e.g., length of the time series)

X = mvnrnd(zeros(N,1), C_con, L);
C_sam = X' * X / L; % sample covariance matrix generated by the estimated configuration model
if transform_to_corr_mat == 1 % transform the generated sample covariance matrix to the correlation matrix
    C_sam = diag(diag(C_sam))^(-1/2) * C_sam * diag(diag(C_sam))^(-1/2);
end
K_sam = inv(C_sam); % sample precision matrix