function chi2 = GompertzRandomLogrank(cumsurv1,cumsurv2,maxduration,treatment,RR,startportion)
%main use:
%100*mean(GompertzRandomLogrank(1000,50,1,1,0.9354)>3.84) 
%proba that two groups of 50 mice aged 500 days to be different

%returns the chi2(dll=1) value for thte Logrank test.
%if it is more than 3.84, then P<0.05: les courbes sont differentes

%GompertzRandomLogrank([91 91 90 90 87 83 83 78],[142 142 142 141 141 141 140 134])
%n=5;GompertzRandomLogrank(GompertzRandomSurvival([0.1 0.005 5],n),GompertzRandomSurvival([0.1 0.005 5],n))

maxTime=max(length(cumsurv1),length(cumsurv2));
for i=length(cumsurv1)+1:maxTime;cumsurv1(i)=0;end
for i=length(cumsurv2)+1:maxTime;cumsurv2(i)=0;end



if maxTime==1
    p=[0.0003 0.000006 35]; % C elegans
    p=[0 0.000012 150]; % Wistar
    p=[0 0.000007 120]; % C57BL6

    chi2=zeros(1,cumsurv1);
    for i=1:cumsurv1
        chi2(i)=GompertzRandomLogrank(GompertzRandomSurvival(p,cumsurv2,1,1,startportion,maxduration),GompertzRandomSurvival(p,cumsurv2,treatment,RR,startportion,maxduration),maxduration);
    end
    %disp( sum(chi2>3.84)/cumsurv1 );
    return;
    %[0.1 0.005 5]. pourcentage de fois où c'est jugé différent (chi2 > 3.84, 10, 15 :)


    %******** SI CHAQUE JOUR COMPTE. Wistar N=25 546+400 jours *********
'Faux positifs';
chisq=GompertzRandomLogrank(100,25,1.00);disp([sum(chisq>1),sum(chisq>2),sum(chisq>3),sum(chisq>3.84),sum(chisq>4),sum(chisq>5)]);%24 13 10  6  5  3
chisq=GompertzRandomLogrank(100,25,1.00);disp([sum(chisq>1),sum(chisq>2),sum(chisq>3),sum(chisq>3.84),sum(chisq>4),sum(chisq>5)]);%19 11  6  5  5  2

chisq=GompertzRandomLogrank(100,25,1.00);disp([sum(chisq>1),sum(chisq>2),sum(chisq>3),sum(chisq>3.84),sum(chisq>4),sum(chisq>5)]);%24 10  3  2  2  2
'si étude jusqu au bout:';

'déceler mal 15% de modif';
chisq=GompertzRandomLogrank(100,25,0.85);disp([sum(chisq>1),sum(chisq>2),sum(chisq>3),sum(chisq>3.84),sum(chisq>4),sum(chisq>5)]);%53 18 10  5  3  0 coef*(t+546)
chisq=GompertzRandomLogrank(100,25,1.15);disp([sum(chisq>1),sum(chisq>2),sum(chisq>3),sum(chisq>3.84),sum(chisq>4),sum(chisq>5)]);%93 87 78 75 74 66 coef*(t+546)
chisq=GompertzRandomLogrank(100,25,1.15);disp([sum(chisq>1),sum(chisq>2),sum(chisq>3),sum(chisq>3.84),sum(chisq>4),sum(chisq>5)]);%52 26 20 17 17 14 (coef*t+546)
chisq=GompertzRandomLogrank(100,25,0.85);disp([sum(chisq>1),sum(chisq>2),sum(chisq>3),sum(chisq>3.84),sum(chisq>4),sum(chisq>5)]);%30 17 8   5  5  3 coef*(t+546)

chisq=GompertzRandomLogrank(100,25,0.85);disp([sum(chisq>1),sum(chisq>2),sum(chisq>3),sum(chisq>3.84),sum(chisq>4),sum(chisq>5)]);%27 15  7  2  1  0 coef*t+546
chisq=GompertzRandomLogrank(100,25,1.15);disp([sum(chisq>1),sum(chisq>2),sum(chisq>3),sum(chisq>3.84),sum(chisq>4),sum(chisq>5)]);%49 33 18 12 12  9 coef*t+546
'si étude jusqu au bout :';
chisq=GompertzRandomLogrank(100,25,0.85);disp([sum(chisq>1),sum(chisq>2),sum(chisq>3),sum(chisq>3.84),sum(chisq>4),sum(chisq>5)]);%92 88 85 81 80 77 coef*t+546
chisq=GompertzRandomLogrank(100,25,1.15);disp([sum(chisq>1),sum(chisq>2),sum(chisq>3),sum(chisq>3.84),sum(chisq>4),sum(chisq>5)]);%89 82 79 71 70 66 coef*t+546
    
    
    
    %******** SI CHAQUE JOUR COMPTE. [0.0003 0.000006 35] *********
'20 animaux décèlent mal 20% de modif';
chisq=GompertzRandomLogrank(100,20,1.2);disp([sum(chisq>10),sum(chisq>15),sum(chisq>25)]);%1 96 94 81
chisq=GompertzRandomLogrank(100,20,1.1);disp([sum(chisq>10),sum(chisq>15),sum(chisq>25)]);%0.92 71% 54% 34%
chisq=GompertzRandomLogrank(100,20,1.0);disp([sum(chisq>10),sum(chisq>15),sum(chisq>25)]);%0.28 9% 6% 1%
chisq=GompertzRandomLogrank(100,20,1.2);disp([sum(chisq>10),sum(chisq>15),sum(chisq>25)]);%0.99 97 94 85
chisq=GompertzRandomLogrank(100,20,1.1);disp([sum(chisq>10),sum(chisq>15),sum(chisq>25)]);%0.90 77% 66% 40%
chisq=GompertzRandomLogrank(100,20,1.0);disp([sum(chisq>10),sum(chisq>15),sum(chisq>25)]);%0.23 5% 2% 1%

    %******** SI CHAQUE JOUR COMPTE. [0.1 0.005 5] et pas [0.0003 0.000006 35] *********
'Faux positifs';
chisq=GompertzRandomLogrank(100,5,1);disp([sum(chisq>10),sum(chisq>15),sum(chisq>25)]);%0.29 11% 3%
chisq=GompertzRandomLogrank(100,10,1);disp([sum(chisq>10),sum(chisq>15),sum(chisq>25)]);%0.32 9% 2% 0%
chisq=GompertzRandomLogrank(100,15,1);disp([sum(chisq>10),sum(chisq>15),sum(chisq>25)]);%0.24 10% 2% 0%
chisq=GompertzRandomLogrank(100,20,1);disp([sum(chisq>10),sum(chisq>15),sum(chisq>25)]);%0.22 6% 1% 0%

chisq=GompertzRandomLogrank(100,5,0.5);disp([sum(chisq>10),sum(chisq>15),sum(chisq>25)]);%0.61 29% 11% 1%

'20 animaux décèlent mal 20% de modif';
chisq=GompertzRandomLogrank(100,20,0.8);disp([sum(chisq>10),sum(chisq>15),sum(chisq>25)]);%0.49 17% 11% 2%
chisq=GompertzRandomLogrank(100,20,0.9);disp([sum(chisq>10),sum(chisq>15),sum(chisq>25)]);%0.34 7% 1% 0%
chisq=GompertzRandomLogrank(100,20,1.1);disp([sum(chisq>10),sum(chisq>15),sum(chisq>25)]);%0.27 11% 6% 0%
chisq=GompertzRandomLogrank(100,20,1.2);disp([sum(chisq>10),sum(chisq>15),sum(chisq>25)]);%0.43 15% 8% 2%

    
    %******** SI CHAQUE JOUR COMPTE. [0.1 0.005 5] *********
'Quand les paramètres sont les memes il faut plutot prendre chi2>15';
chisq5=GompertzRandomLogrank(100,5,1); %0.31/0.29/0.27, 12%/7%, 3%/3%
chisq10=GompertzRandomLogrank(100,10,1); %0.29/0.31, 11%, 5%
chisq15=GompertzRandomLogrank(100,15,1); %0.21/0.28, 9%
chisq20=GompertzRandomLogrank(100,20,1); %0.25,0.19 7%/7%, 5%/2%

'20 animaux décèlent mal 20% de modif';
chisq=GompertzRandomLogrank(100,20,0.8);disp([sum(chisq>10),sum(chisq>15)]);%0.94 79% 68%
chisq=GompertzRandomLogrank(100,20,0.9);disp([sum(chisq>10),sum(chisq>15)]);%0.53 23% 14%
chisq=GompertzRandomLogrank(100,20,1.1);disp([sum(chisq>10),sum(chisq>15)]);%0.59 33% 22%
chisq=GompertzRandomLogrank(100,20,1.2);disp([sum(chisq>10),sum(chisq>15)]);%0.90 70% 54%

'5 animaux décèlent très très mal 20% de modif et peu 50% de modif';
chisq5=GompertzRandomLogrank(100,5,0.5); %0.98, 83%, 68%
chisq5=GompertzRandomLogrank(100,5,0.8); %0.58, 25%, 16%
chisq5=GompertzRandomLogrank(100,5,0.9); %0.36,0.40, 16%13%, 6%/8%
chisq5=GompertzRandomLogrank(100,5,1.1); %0.26/0.35, 6%/15%, 3%/4%
chisq5=GompertzRandomLogrank(100,5,1.2); %0.5, 22%, 10%

'10 animaux décèlent très mal 20% de modif';
chisq=GompertzRandomLogrank(100,10,0.8);disp([sum(chisq>10),sum(chisq>15)]);%0.70 46% 30%
chisq=GompertzRandomLogrank(100,10,0.9);disp([sum(chisq>10),sum(chisq>15)]);%0.48 26% 14%
chisq=GompertzRandomLogrank(100,10,1.1);disp([sum(chisq>10),sum(chisq>15)]);%0.39 16% 7%
chisq=GompertzRandomLogrank(100,10,1.2);disp([sum(chisq>10),sum(chisq>15)]);%0.69 40% 23%

'15 animaux décèlent mal 20% de modif';
chisq=GompertzRandomLogrank(100,15,0.8);disp([sum(chisq>10),sum(chisq>15)]);%0.85 70% 48%
chisq=GompertzRandomLogrank(100,15,0.9);disp([sum(chisq>10),sum(chisq>15)]);%0.60 30% 16%
chisq=GompertzRandomLogrank(100,15,1.1);disp([sum(chisq>10),sum(chisq>15)]);%0.55 23% 8%
chisq=GompertzRandomLogrank(100,15,1.2);disp([sum(chisq>10),sum(chisq>15)]);%0.83 55% 32%


    %******** SI CHAQUE EVENEMENT COMPTE. [0.1 0.005 5] *********
'Faux positifs';
chisq=GompertzRandomLogrank(100,5,1);disp([sum(chisq>10),sum(chisq>15),sum(chisq5>25)]);%0.46 25% 15% 6% (4%>30)
chisq=GompertzRandomLogrank(100,10,1);disp([sum(chisq>10),sum(chisq>15),sum(chisq5>25)]);%0.43 25% 19% 6%
chisq=GompertzRandomLogrank(100,15,1);disp([sum(chisq>10),sum(chisq>15),sum(chisq5>25)]);%0.44 25% 17% 6%
chisq=GompertzRandomLogrank(100,20,1);disp([sum(chisq>10),sum(chisq>15),sum(chisq5>25)]);%0.47 19% 10% 6%

'20 animaux décèlent mal 20% de modif';
chisq=GompertzRandomLogrank(100,20,0.8);disp([sum(chisq>10),sum(chisq>15),sum(chisq>25)]);%0.97 91% 89%
chisq=GompertzRandomLogrank(100,20,0.9);disp([sum(chisq>10),sum(chisq>15),sum(chisq>25)]);%0.81 68% 54%
chisq=GompertzRandomLogrank(100,20,1.1);disp([sum(chisq>10),sum(chisq>15),sum(chisq>25)]);%0.50 19% 12%
chisq=GompertzRandomLogrank(100,20,1.2);disp([sum(chisq>10),sum(chisq>15),sum(chisq>25)]);%0.73 53% 38% 13%

Ea=0;
Eb=0;
Exposesa=cumsurv1(1);
Exposesb=cumsurv2(1);
Exposes=Exposesa+Exposesb;
for i=2:maxTime
    if Exposes>0
        morts= Exposes-cumsurv1(i)+cumsurv2(i);
        if morts>0
            Ea = Ea + Exposesa * morts/Exposes; %AttendusGpe1
            Eb = Eb + Exposesb * morts/Exposes; %AttendusGpe2
            %disp([AttendusGpe1 AttendusGpe2]);
            Exposesa=cumsurv1(i);
            Exposesb=cumsurv2(i);
            Exposes=Exposesa+Exposesb;
        end
    end
end
%disp([Ea,Eb]);
chi2= (Ea-Eb)^2/(Ea*Eb/(Ea+Eb));


Ea=0;
Eb=0;
for i=2:maxTime
    exposes = cumsurv1(i-1)+cumsurv2(i-1);
    if exposes>0
        mortalite = 1-(cumsurv1(i)+cumsurv2(i))/exposes;
        Ea = Ea + cumsurv1(i-1) * mortalite;
        Eb = Eb + cumsurv2(i-1) * mortalite;
        %disp([AttendusGpe1 AttendusGpe2]);
    end
end
%disp([Ea,Eb]);
chi2= (Ea-Eb)^2/(Ea*Eb/(Ea+Eb));


gain_de_vie=0;
for i=2:maxTime
    gain_de_vie = gain_de_vie + (cumsurv2(i)-cumsurv1(i));
end
%disp([Ea,Eb]);
chi2= gain_de_vie;

end


%figure;
%plot(1:maxTime, cumsurv1,'g',1:maxTime, cumsurv2,'r')

%maxTime=min(maxTime,400); %pour le cas où je finis mi-juillet 2006.

Ea=0;
Eb=0;
Exposesa=cumsurv1(1);
Exposesb=cumsurv2(1);
Exposes=Exposesa+Exposesb;
for i=2:maxTime
    if Exposes>0
        morts= Exposes-(cumsurv1(i)+cumsurv2(i));
        if morts>0
            Ea = Ea + Exposesa * morts/Exposes; %AttendusGpe1
            Eb = Eb + Exposesb * morts/Exposes; %AttendusGpe2
            %disp([AttendusGpe1 AttendusGpe2]);
            Exposesa=cumsurv1(i);
            Exposesb=cumsurv2(i);
            Exposes=Exposesa+Exposesb;
        end
    end
end
%disp([Ea,Eb]);
Oa = cumsurv1(1) - Exposesa;
chi2= (Oa-Ea)^2/(Ea*Eb/(Ea+Eb));


%if chi2>3.84;resultat='different';else;resultat='pas different';end 
%title(['n=' num2str(mean(cumsurv1(1),cumsurv2(1))) '  chi2=' num2str(chi2) '  : ' resultat]);

return