%% CR Model: Abundance Correlation vs D for S = 2, Multiple Qs
clear; clc;

%% Parameters
alpha = 0.5; sigma = 0.0; sigma_e = 0.05;
steps = 100; dt = 1/steps; T = 5000;
epsilon = 1e-8; penalty_weight = 100;
S = 2; Q_vals = [3,4,5,6]; runs = 900;
N_tau = 10; % Number of resource updates per consumer step
Alpha_Matrix = [1 alpha; alpha 1];
colors = lines(length(Q_vals));
all_D = {}; all_corr = {}; labels = {};

%% Simulation Loop
for q_idx = 1:length(Q_vals)
    Q = Q_vals(q_idx);
    D_list = []; C_list = [];

    for run = 1:runs
        fprintf('S=%d, Q=%d, run %d/%d\n', S, Q, run, runs);
        [gamma, lambda, D] = optimize_gamma_lambda_Angle(Alpha_Matrix, Q, penalty_weight);
        
        x = rand(1, S); 
        R = rand(1, Q);
        history = zeros(S, T); 
        growth_rate = zeros(S, T);

        for t = 1:T
            x_prev = x;

            for step = 1:steps
                % Consumer update
                x = x + dt * ((-1 + gamma * R') .* x')' + epsilon;

                % Fast resource updates
                for tau = 1:N_tau
                    R = R + dt * (R - R.^2 - R .* (x * lambda') + R * sigma_e^2 / 2) ...
                        + randn(1, Q) .* R * sqrt(dt) * sigma_e;
                end
            end

            history(:,t) = x'; 
            growth_rate(:,t) = log(x ./ x_prev)';
        end

        % Correlation and storage
        C = corr(history(:,T/10:end)');
        D_list = [D_list; D];
        C_list = [C_list; C(1,2)];
    end

    all_D{q_idx} = D_list;
    all_corr{q_idx} = C_list;
    labels{q_idx} = sprintf('Q = %d', Q);
end

%% Plot Results
figure; hold on;
for k = 1:length(Q_vals)
    scatter(all_D{k}, all_corr{k}, 50, 'filled', 'MarkerFaceColor', colors(k,:));
end
xlabel('D '); ylabel('Growth Rate Correlation');
title('Growth Correlation vs D ');
legend(labels, 'Location', 'best'); grid on;

%% Final Run Dynamics
figure; semilogy(history', 'LineWidth', 1.5);
xlabel('Time'); ylabel('Consumer Population');
title(sprintf('Final Simulation: S = 2, Q = %d', Q_vals(end)));
legend('Species 1', 'Species 2'); grid on;