% This program generates large pool randomized parameters and simulation result
clear
clc

addpath('./Tools/');
% color for pure-producer; partial-producer; pure-cheater
strategy_colors=[0,0.8,0;
    1,0.8,0;
    0.8,0,0];

rng('shuffle');  % Set to new time-based seed
state = rng;      % Get the current (new) state
randomSeed = state.Seed;
%% Run simulation
cyclenum=1000; % how many cycles to generate

N_sid=50; % Number of siderophores
N_spe=50; % Number of species

% parameter set for the ODE
paraset=[];
paraset.e= 10 * ones(N_sid,1); % siderophore conversion efficiency， 10 比1 更快进入limit cycle
paraset.u= 1 * ones(N_sid,1); % Iron-binding affinity， 10 比1慢
paraset.migr=10^-8 *1; % migration rate
paraset.gamma=1; % scale for growth rate
paraset.d=0.1; % dilution rate
paraset.R_sup=3;% 少一点似乎更快，但多了至少也会震,少了容易灭绝
%%%%%%%%%
M_scale= paraset.gamma * paraset.R_sup; % scale for the biomass
R_scale=mean(paraset.e(:))*M_scale; % scale of the iron resource

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% folder to store the data
typestr=[num2str(randomSeed),'_',num2str(paraset.R_sup),'_',num2str(N_spe),'_',num2str(N_sid)];
foldername=['largepool_simulation/v2_largepool_simulation_nosaving_',typestr];
if ~exist(foldername, 'dir')
    % Create folder if it doesn't exist
    mkdir(foldername);
    disp(['Folder "', foldername, '" created successfully.']);
else
    disp(['Folder "', foldername, '" already exists.']);
end
% record the result:
resultfilepath=[foldername,'/','ResultSummary_',typestr,'.csv'];
fid = fopen(resultfilepath, 'a');

CheatingCountlist=(N_spe*(0:1:min(100,N_sid-1))); % Number of Cheating Receptors
p_cheater_list=(0:1:5)/N_spe; % percentage of pure-cheaters
for cyc=1:cyclenum % different repeating cycles
    for cc=1:length(CheatingCountlist)
        CheatingCount=CheatingCountlist(cc);
        for p=1:length(p_cheater_list)
            p_cheater=p_cheater_list(p);
            [alpha_sid_pro, v_sid_rec] = Generate_Random_IronNet(N_spe,N_sid,p_cheater,CheatingCount);
            
            % get the actual r_c and CB
            cheaternum=sum(sum(alpha_sid_pro,2)==0);
            producernum=(N_spe-cheaternum);
            Sid_beingproduced=find(sum(alpha_sid_pro>0)>0);
            totalRec=sum(sum(v_sid_rec(:,Sid_beingproduced)>0));
            CR_num=totalRec-producernum;
            actualCB=CR_num/N_spe;
            
            
            
            filenamestr=[num2str(N_spe),'_',num2str(N_sid),'_',num2str(p_cheater),'_',num2str(CheatingCount),'_',num2str(cyc)];
            parafilepath=[foldername,'/',filenamestr,'.mat'];
            
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            %%% define the ODE s
            f = ODE_siderophore_iron_partition(alpha_sid_pro,v_sid_rec,paraset);
            % randomized initial condition
            randinit=[M_scale*rand(N_spe,1)/N_spe;R_scale*rand(N_sid,1);paraset.R_sup];
            
            %%%%%%%%%%%%%%%%%%%%%%%%%%%
            clear otherp
            % parameters to assess the attractor type
            otherp.roundmax=5;% how many rounds for maximal assesement
            otherp.ssthresh=10^-3;  % threshold for steady-state
            otherp.varthresh=10;
            % threshold for limit cycle
            otherp.dominant_thresh=0.3;
            otherp.entropy_thresh=0.3;
            otherp.harmonic_thresh=10^-3;
            % simulation length for the first round
            otherp.onesimulationspan=3*10^4;
            
            
            %%%%%%%%%%%%%%%%%%%%%%%%%%%
            [meanvalue,t,z,otheroutput] = Simulate2attractor_siderophore(f,randinit,otherp);
            %%%%%%%%%%%%%%%%%%%%%%% record the result of parameter and simulation
            
            fprintf(fid,'%s,%d,%d,%f,%d,',filenamestr,N_spe,N_sid,cheaternum,CR_num);
            fprintf(fid,'%d,%d,%d,',otheroutput.r, otheroutput.issteady,otheroutput.isolimitcycle);
            fprintf(fid,'%f,',otheroutput.firstperiod,otheroutput.dominant_ratio,otheroutput.spectral_entropy,otheroutput.harmonic_dist);
            fprintf(fid,'%f,',meanvalue);
            % record last time values
            fprintf(fid,'%f,',z(end,:));
            fprintf(fid,'\n');
            
            % save the parameter if needed
            if 0%sum(meanvalue(1:N_spe)>10^-3)>8 % surviving species number
                save_odeparameters(parafilepath, alpha_sid_pro, v_sid_rec, paraset);
            end
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            if p==1
                fprintf('%d\t%d\t%d\tround(%d)\tissteady(%d)\tbiomasstotal(%f)\n',cyc,cc,p,otheroutput.r,otheroutput.issteady,sum(meanvalue(1:N_spe)))
            end %for p=1:length(p_cheater_list)
        end
    end
end
fclose(fid);