exec('/Volumes/shared_data/paulb/Scilab/Infected/infected_init_local.sci'); // Load the models and the data
exec('/Volumes/shared_data/paulb/Scilab/Infected/infected_ensemble.sci'); // Load the multi-model
exec('/Volumes/shared_data/paulb/Scilab/Verification/Verification_Dates.sce'); // Load the functions start_dates and end_dates

function [cor]=correction(param,N_y,N_e,XXX)
    cor=[0,0,0,0,0,0,0,0,0,0,0,0];
    I=[]; cas=[];
    for j=1:N_y // Years 80 to 03 (to avoid limits)

        // Define starting (april) en ending dates
        ti=Start_Dates(j,4); //See comment on file Verification_Dates.sce
        tf=End_Dates(j,15); // Here, we mak prediction for one year, starting on the first of April
        if tf-ti+1>365 then // If bisextile year
            tf=tf-1;
        end

        // Define predicted, actual and climatology cases
        EI=infected_ensemble(param,ti,tf,XXX); // Predicted cases for one year
        Cas=Cases(ti:tf,2); // Actual cases for the same period of time
        Clima=Climatology(ti,tf); // Climatology of the actual cases for the same period of time

        // Take the monthly mean of them
        EI_mean=zeros(12,N_e);
        ss=0;
        for i=1:12 // Monthly average
            gg=End_Dates(j,i)-Start_Dates(j,i)+1;
            if ss+gg>365 then gg=gg-1; end // For bisextiles years
            EI_mean(i,:)=mean(EI(ss+1:ss+gg,:),'r'); // monthly average of the predicted cases
            Cas_month(i)=Cas((ss+1+ss+gg)/2); // monthly prediction for the actual cases
            Clima_mean(i)=mean(Clima(ss+1:ss+gg)); // monthly average of the climatology
            ss=ss+gg;
        end

        //    // substract the yearly annomaly
        //    EI_mean_ann=zeros(12,N_e)
        //    for k=1:N_e // Computation of the annomalies
        //        EI_mean_ann(:,k)=EI_mean(:,k)-mean(EI_mean(:,k)); // annomalies of the predicted cases
        //    end
        //    Cas_month_ann=Cas_month-mean(Cas_month); // Annomalies of the actual cases
        //    Clima_mean_ann=Clima_mean-mean(Clima_mean);

        // Compute the correction for the predicted cases
        for L=1:12 // Computation of the 
            cor(L)=cor(L)+Cas_month(L)-mean(EI_mean(L,:));
        end

        // Store the results
        I=[I,EI_mean];
        cas=[cas,Cas_month];
        Clima=[Clima_mean];
    end
    cor=cor ./ N_y // Divide by N_y
    
    //We put the months in the right order (Jan to Dec)
    cor=[cor(10:12),cor(1:9)];
endfunction