close all;
clc;
clear;

filename0 = 'A.txt';
path0 = '/Users/sergeymakovkin/Documents/Kuramoto_network/Solver4_experiment/Solver/Solver';
Data = importdata( [ path0 '/' filename0 ] );

%figure;
figure('units', 'normalized', 'outerposition', [0 0 1 1]);
imagesc(Data, [0 1]);
%axis equal;
daspect([1 1 1])
%%axis off;
colormap(gray(2));
caxis([0 1]);
cbh = colorbar('v');
set(cbh, 'YTick', [0:1:1]);
set(gca, 'XTick', [1 size(Data, 1)/2 size(Data, 1)]);
set(gca, 'YTick', [1 size(Data, 2)/2 size(Data, 2)]);
ax = gca;
ax.XAxisLocation = 'top';
title('Adjacency matrix');
set(gca, 'FontSize', 40);
saveas(gcf, 'A', 'epsc');

filename3 = 'links.txt';
path3 = '/Users/sergeymakovkin/Documents/Kuramoto_network/Solver4_experiment/Solver/Solver';
Data = importdata( [ path3 '/' filename3 ] );

% f_n_omega = figure('units', 'normalized', 'outerposition', [0 0 1 1]);
% subplot(3, 1, 1);
figure;
oscillator_number = linspace(size(Data, 1)+1, 2*size(Data, 1), size(Data, 1));
bar(oscillator_number, Data(:, 1), 'r');
hold on;
plot([size(Data, 1)+0.5 2*size(Data, 1)+0.5], [mean(Data(:, 1)) mean(Data(:, 1))], 'b', 'LineWidth', 3);
hold on;
xlim([size(Data, 1)+0.5, 2*size(Data, 1)+0.5]);
set(gca, 'XTick', [size(Data, 1)+1, 2*size(Data, 1)]);
title('Links counter diagram');
xlabel('Neuron layer oscillator number');
ylabel('Links count');
saveas(gcf, 'Links', 'epsc');

filename1 = 'omega_init.txt';
path1 = '/Users/sergeymakovkin/Documents/Kuramoto_network/Solver4_experiment/Solver/Solver';
omega_init = importdata( [ path1 '/' filename1 ] );
oscillator_number = linspace(1, size(omega_init, 1), size(omega_init, 1));

f_omega_1 = figure('units', 'normalized', 'outerposition', [0 0 1 1]);
subplot(1, 2, 1);
plot(oscillator_number(1:end/2),  omega_init(1:end/2)', 'r.', 'MarkerSize', 21);
hold on;
plot(oscillator_number(end/2+1:end),  omega_init(end/2+1:end)', 'b.', 'MarkerSize', 21);
grid on;
xlim([0 size(oscillator_number, 2)]);
ylim([0 4]);
set(gca, 'XTick', 0:size(oscillator_number, 2)/2:size(oscillator_number, 2));
set(gca, 'FontSize', 28);
title('\omega_{init}, natural frequencies of i^{th} oscillator', 'FontSize', 28);
xlabel('Oscillator number', 'FontSize', 50);
ylabel('\omega, Hz', 'FontSize', 50);
legend('\omega_{init wired glial}', '\omega_{init random neuron}', 'Location', 'northwest', 'FontSize', 28);

% figure(f_n_omega);
% subplot(3, 1, 2);
% plot(oscillator_number(end/2+1:end)+1, omega_init(end/2+1:end), 'b.', 'MarkerSize', 21);
% grid on;
% grid minor;
% xlim([oscillator_number(end/2+1)+1 oscillator_number(end)+1]);
% ylim([3 4]);
% set(gca, 'XTick', [oscillator_number(end/2+1)+1, oscillator_number(end)+1]);
% xlabel('Neuron layer oscillator number');
% ylabel('\omega_{init random neuron}, Hz');

filename2 = 'results.txt';
path2 = '/Users/sergeymakovkin/Documents/Kuramoto_network/Solver4_experiment/Solver/Solver';
Data = importdata( [ path2 '/' filename2 ] );
time = Data(:, 1);
theta = Data(:, 2:end);

figure('units', 'normalized', 'outerposition', [0 0 1 1]);
plot(time, theta);
%plot(time, mod(theta, 2*pi));
%grid on;
%title('Kuramoto C results on MacBook Pro');
set(gca, 'FontSize', 24);
xlabel('Time, msec', 'FontSize', 50);
ylabel('\theta, rad', 'FontSize', 50);
%saveas(gcf, 'results.png');
saveas(gcf, 'results', 'epsc');

%%%%%%%%%%%%%%%%%%%%%%%%%%%
% omega calculation HARD!!!
% for i = 0.5:0.1:0.6
%     omega_end_distribution(i, :) = ( bsxfun(@minus, theta(i, :), theta(size(theta(round(end/2):end, :), 1), :)) ) / ( time(i) - time(size(theta(round(end/2):end, :), 1)) );
% end
%
% figure('units', 'normalized', 'outerposition', [0 0 1 1]);
% plot(time(:, :), omega_end_distribution(:, :));
% grid on;
% hold on;
%%%%%%%%%%%%%%%%%%%%%%%%%%%

omega_end(:) = ( bsxfun(@minus, theta(end, :), theta(round(end/2), :)) ) / ( time(end) - time(round(end/2)) );
% plot(time(end), omega_end(:), '.', 'MarkerSize', 20);
% grid on;

%%%%%%%%%%%
omega_end_055(:) = ( bsxfun(@minus, theta(round(0.55*end), :), theta(round(end/2), :)) ) / ( time(round(0.55*end)) - time(round(end/2)) );
omega_end_060(:) = ( bsxfun(@minus, theta(round(0.60*end), :), theta(round(end/2), :)) ) / ( time(round(0.60*end)) - time(round(end/2)) );
omega_end_065(:) = ( bsxfun(@minus, theta(round(0.65*end), :), theta(round(end/2), :)) ) / ( time(round(0.65*end)) - time(round(end/2)) );
omega_end_070(:) = ( bsxfun(@minus, theta(round(0.70*end), :), theta(round(end/2), :)) ) / ( time(round(0.70*end)) - time(round(end/2)) );
omega_end_075(:) = ( bsxfun(@minus, theta(round(0.75*end), :), theta(round(end/2), :)) ) / ( time(round(0.75*end)) - time(round(end/2)) );
omega_end_080(:) = ( bsxfun(@minus, theta(round(0.80*end), :), theta(round(end/2), :)) ) / ( time(round(0.80*end)) - time(round(end/2)) );
omega_end_085(:) = ( bsxfun(@minus, theta(round(0.85*end), :), theta(round(end/2), :)) ) / ( time(round(0.85*end)) - time(round(end/2)) );
omega_end_090(:) = ( bsxfun(@minus, theta(round(0.90*end), :), theta(round(end/2), :)) ) / ( time(round(0.90*end)) - time(round(end/2)) );
omega_end_095(:) = ( bsxfun(@minus, theta(round(0.95*end), :), theta(round(end/2), :)) ) / ( time(round(0.95*end)) - time(round(end/2)) );
omega_end_100(:) = ( bsxfun(@minus, theta(round(1.00*end), :), theta(round(end/2), :)) ) / ( time(round(1.00*end)) - time(round(end/2)) );
omega_end_time_distribution = [omega_end_055(:) omega_end_060(:) omega_end_065(:) omega_end_070(:) omega_end_075(:) omega_end_080(:) omega_end_085(:) omega_end_090(:) omega_end_095(:) omega_end_100(:)];
time_omega_end_time_distribution = [time(round(0.55*end)) time(round(0.60*end)) time(round(0.65*end)) time(round(0.70*end)) time(round(0.75*end)) time(round(0.80*end)) time(round(0.85*end)) time(round(0.90*end)) time(round(0.95*end)) time(round(1.00*end))];

figure('units', 'normalized', 'outerposition', [0 0 1 1]);
plot(time_omega_end_time_distribution, omega_end_time_distribution, '.-');
grid on;
ylim([0 4]);
xlabel('Time, msec');
ylabel('\omega_{end random neuron}, Hz');
title('\omega by time distribution');
%%%%%%%%%%%

% figure(f_n_omega);
% subplot(3, 1, 3);
% plot(oscillator_number(end/2+1:end)+1, omega_end(end/2+1:end), 'b.', 'MarkerSize', 21);
% grid on;
% grid minor;
% xlim([oscillator_number(end/2+1)+1 oscillator_number(end)+1]);
% ylim([3 4]);
% set(gca, 'XTick', [oscillator_number(end/2+1)+1, oscillator_number(end)+1]);
% xlabel('Neuron layer oscillator number');
% ylabel('\omega_{end random neuron}, Hz');

figure(f_omega_1);
subplot(1, 2, 2);
plot(oscillator_number(1:end/2),  omega_end(1:end/2), 'r.', 'MarkerSize', 21);
hold on;
plot(oscillator_number(end/2+1:end),  omega_end(end/2+1:end), 'b.', 'MarkerSize', 21);
grid on;
xlim([0 size(oscillator_number, 2)]);
ylim([0 4]);
set(gca, 'XTick', 0:size(oscillator_number, 2)/2:size(oscillator_number, 2));
set(gca, 'FontSize', 28);
title('\omega_{end average} in the last time point, natural frequencies of i^{th} oscillator', 'FontSize', 28);
xlabel('Oscillator number', 'FontSize', 50);
ylabel('\omega, Hz', 'FontSize', 50);
legend('\omega_{end wired glial}', '\omega_{end random neuron}', 'Location', 'northwest', 'FontSize', 28);
saveas(gcf, 'omega_disrtibution_1', 'epsc');

% Dispersion: mean((x - mean(x)).^2) or var(x, 1)
Delta_total = sqrt(var(omega_end(:), 1));
Delta_wired_glial = sqrt(var(omega_end(1:end/2), 1));
Delta_random_neuron = sqrt(var(omega_end(end/2+1:end), 1));

% fid = fopen('Delta_N25_s1.000.txt', 'at');
% fprintf(fid, '%f\t %f\t %f \n', [Delta_total Delta_wired_glial Delta_random_neuron]');
% fclose(fid);
%%%%%%%%%

% Answer here: http://stackoverflow.com/questions/15754459/how-can-i-make-a-color-map-plot-in-matlab
% surf(i, j, B);
% colorbar;
% view(90, 90);
omega_init_wired_glial = omega_init(1:end/2);
omega_init_random_neuron = omega_init(end/2+1:end);

image_init_wired_glial = reshape(omega_init_wired_glial, [sqrt(size(omega_init_wired_glial, 1)), sqrt(size(omega_init_wired_glial, 1))]);
image_init_random_neuron = reshape(omega_init_random_neuron, [sqrt(size(omega_init_random_neuron, 1)), sqrt(size(omega_init_wired_glial, 1))]);

f_omega_2 = figure('units', 'normalized', 'outerposition', [0 0 1 1]);
subplot(1, 2, 1);
imagesc(image_init_wired_glial');
axis equal;
axis off;
colormap jet;
caxis([0 4]);
colorbar;
title(['\omega_{init wired glial}. Layer size = ' num2str(sqrt(size(omega_init_wired_glial, 1))) 'x' num2str(sqrt(size(omega_init_wired_glial, 1)))]);
set(gca, 'FontSize', 40);

figure(f_omega_2);
subplot(1, 2, 2);
imagesc(image_init_random_neuron');
axis equal;
axis off;
colormap jet;
caxis([0 4]);
colorbar;
title(['\omega_{init random neuron}. Layer size = ' num2str(sqrt(size(omega_init_random_neuron, 1))) 'x' num2str(sqrt(size(omega_init_random_neuron, 1)))]);
set(gca, 'FontSize', 40);
saveas(gcf, 'omega_disrtibution_2', 'epsc');

omega_end_wired_glial = omega_end(1:end/2);
omega_end_random_neuron = omega_end(end/2+1:end);

image_end_wired_glial = reshape(omega_end_wired_glial, [sqrt(size(omega_end_wired_glial, 2)), sqrt(size(omega_end_wired_glial, 2))]);
image_end_random_neuron = reshape(omega_end_random_neuron, [sqrt(size(omega_end_random_neuron, 2)), sqrt(size(omega_end_wired_glial, 2))]);

f_omega_3 = figure('units', 'normalized', 'outerposition', [0 0 1 1]);
subplot(1, 2, 1);
imagesc(image_end_wired_glial');
axis equal;
axis off;
colormap jet;
caxis([0 4]);
colorbar;
title(['\omega_{end wired glial}. Layer size = ' num2str(sqrt(size(omega_end_wired_glial, 2))) 'x' num2str(sqrt(size(omega_end_wired_glial, 2)))]);
set(gca, 'FontSize', 40);

figure(f_omega_3);
subplot(1, 2, 2);
imagesc(image_end_random_neuron');
axis equal;
axis off;
colormap jet;
caxis([0 4]);
colorbar;
title(['\omega_{end random neuron}. Layer size = ' num2str(sqrt(size(omega_end_random_neuron, 2))) 'x' num2str(sqrt(size(omega_end_random_neuron, 2)))]);
set(gca, 'FontSize', 40);
%saveas(gcf, 'omega_disrtibution.png');
saveas(gcf, 'omega_disrtibution_3', 'epsc');

filename4 = 'results.txt';
path4 = '/Users/sergeymakovkin/Documents/Kuramoto_network/Solver4_experiment/Solver/Solver';
Data = importdata( [ path4 '/' filename4 ] );
time = Data(:, 1);
theta = Data(:, 2:end);

rx_total = (1/(size(theta, 2)))*sum(cos(theta(:, :)')); % Sum of mean x-part of theta
ry_total = (1/(size(theta, 2)))*sum(sin(theta(:, :)')); % Sum of mean y-part of theta
r_total = sqrt(rx_total.*rx_total + ry_total.*ry_total);
r_total_mean = mean(r_total(:, round(size(r_total, 2)*0.5):end));

rx_wired_glial = (1/(size(theta(:, 1:end/2), 2)))*sum(cos(theta(:, 1:end/2)')); % Sum of mean x-part of theta
ry_wired_glial = (1/(size(theta(:, 1:end/2), 2)))*sum(sin(theta(:, 1:end/2)')); % Sum of mean y-part of theta
r_wired_glial = sqrt(rx_wired_glial.*rx_wired_glial + ry_wired_glial.*ry_wired_glial);
r_wired_glial_mean = mean(r_wired_glial(:, round(size(r_wired_glial, 2)*0.5):end));

rx_random_neuron = (1/(size(theta(:, end/2+1:end), 2)))*sum(cos(theta(:, end/2+1:end)')); % Sum of mean x-part of theta
ry_random_neuron = (1/(size(theta(:, end/2+1:end), 2)))*sum(sin(theta(:, end/2+1:end)')); % Sum of mean y-part of theta
r_random_neuron = sqrt(rx_random_neuron.*rx_random_neuron + ry_random_neuron.*ry_random_neuron);
r_random_neuron_mean = mean(r_random_neuron(:, round(size(r_random_neuron, 2)*0.5):end));

% fid = fopen('r_sigma_N50_sG0.500_sGN0.370_sN0.500.txt', 'at');
% fprintf(fid, '%f\t %f\t %f \n', [r_total_mean r_wired_glial_mean r_random_neuron_mean]');
% %dlmwrite('r_data.txt', [r_total_mean r_wired_glial_mean r_random_neuron_mean], '-append');
% fclose(fid);

figure('units', 'normalized', 'outerposition', [0 0 1 1]);
%plot(time, rx_total, 'g');
%hold on;
%plot(time, ry_total, 'b');
%hold on;
plot(time, r_total, 'k', time, r_wired_glial, 'r', time, r_random_neuron, 'b');
%plot(time, r_wired_glial, 'r', time, r_random_neuron, 'b');
%title('Kuramoto C results on MacBook Pro');
hold on;
line([time(round(size(r_total, 2)*0.5)) time(end)], [r_total_mean r_total_mean], 'Color', 'k', 'LineWidth', 2);
hold on;
line([time(round(size(r_wired_glial, 2)*0.5)) time(size(r_wired_glial, 2))], [r_wired_glial_mean r_wired_glial_mean], 'Color', 'r', 'LineWidth', 2);
hold on;
line([time(round(size(r_random_neuron, 2)*0.5)) time(size(r_random_neuron, 2))], [r_random_neuron_mean r_random_neuron_mean], 'Color', 'b', 'LineWidth', 2);
hold on;
%grid on;
legend('r_{total}', 'r_{wired glial}', 'r_{random neuron}', 'r_{total mean}', 'r_{wired glial mean}', 'r_{random neuron mean}', 'Location', 'east');
%legend('r_{wired glial}', 'r_{random neuron}', 'r_{wired glial mean}', 'r_{random neuron mean}', 'Location', 'east');
set(gca, 'FontSize', 28);
xlabel('Time, msec', 'FontSize', 50);
ylabel('Mean field', 'FontSize', 50);
ylim([0 1]);
%saveas(gcf, 'r_order.png');
saveas(gcf, 'r_order', 'epsc');