%
% FIXSEIS.M;1    2/4        12-JUN-2006 18:44:44.06
%
% complsim: complete simulation of 3 yrs of seismicity
% with all tests.
%
% parameters
%
clear all
fs1 = 15;
nx = 500;
x = linspace(0,1,nx);
y1 = linspace(0,1,nx);
y = fliplr(y1);
%
rate=15;
al=1;
bet=1/rate;
%
gmc=1;
% gmc=2;
% gmc=3;
% gmc=5;
% gmc=10;
% gmc=20;
%
al2=1/gmc;
bet2=gmc/rate;
% sc=2500
% sc=100
% sc=50
% sc=45
% sc=44
% sc=44.5
% sc=44.4
% sc=42.5
% sc=25
% sc=16
% sc=10
% sc=5
sc=4
% sc=2
scn=2
% scn=5
%
tmax=25000;
% tmax=15000;
% tmax=5000;
% tmax=2500;
tmax=1000;
% tmax=500;
% tmax=250;
% tmax=100;
% tmax=10;
% tmax=3;
% tmax=1.0;
% tmax=0.5;
%
% w=2.5/52.1786;
% w=5/52.1786;
% w=25/52.1786;
% w=52.1786/52.1786;
% w=tmax;
%
% nsim=10000;
% nsim=1000;
% nsim=10;
nsim=1;
%
len=20*tmax;
len=30*tmax;
%
% p=1/2;
p=1;
ws=0.0000000001;
wd=11;
wd=12;
ws=0.000001;
wd=7;
wd=8;
% ws=0.0001;
% wd=5;
% ws=0.00001;
% wd=5;
%
% t=gamseis1(len,al,bet);
% numit=1;
% numit=10;
numit=20;
% numit=33;
numit=50;
numit=100;
% numit=250;
for j=1:numit
%     w1(j)=j*tmax/numit;
%     w1(j)=2*j*bet/numit;
    w1(j)=ws*bet*10^(wd*(j-1)/numit);
    w=w1(j);
%     w=0.1
%     t=gamseis1(len,al2,bet2);
%     t=gamseis2(len,al2,bet2,sc,scn);
    t=gamseis3(len,al2,bet2,sc,scn);  %Random permutation
    [s,l, np, nu, tau]=pred7(t,w,tmax,p);
    nu1(j)=1-nu;
    tau1(j)=tau;
end
nu2s=[1, nu1, 0];
tau2s=[0, tau1, 1];
%
for j=1:numit+1
    dnu2(j)=nu2s(j)-nu2s(j+1);
    dtau2(j)=tau2s(j+1)-tau2s(j);
%     t(t >= tmax)=[];
    if dnu2(j) < 0
        dnu2(j)=0
    end
        if dtau2(j) < 0
        dtau2(j)=0
    end
    inf(j)=(log(dnu2(j)./dtau2(j))/log(2)).*dnu2(j);
end
inf(isnan(inf))=[0];
inf(isinf(inf))=[0];
infs=sum(inf);
kap = 1/gmc;
gj=(1-kap-gammaln(kap)+log(kap)+(kap-1)*psi(kap))/log(2);
% numit=100
wd1=wd;
wd1=wd+1;
for j=1:numit
    w1(j)=ws*bet*10^(wd1*(j-1)/numit);
    w=w1(j);
    nu3(j)=gammainc(kap*w,kap);
    nu3j=nu3(j);
%     tau32(j)=(nu3j/(1-nu3j))*gammainc(kap*w,1+kap);
    tau31(j)=w*(1-nu3j);
    tau32(j)=gammainc(kap*w,1+kap);
    tau3(j)=tau31(j)+tau32(j);
%     tau3(j)=tau31(j);
%     tau3(j)=gammainc(kap*w,1+kap);
end
nu3=1-nu3;
nu4=[1, nu3, 0];
tau4=[0, tau3, 1];
% tau4=[0, tau31, 1];
% tau4=[0, tau32, 1];
%
nu0=0.5;
%
numit=13
for j=1:numit
    d=2*10^((j-1)/30);
[nu,fval] = fsolve(@(nu) twop(nu,d),nu0);
% nu
d1(j)=d;
nu1o(j)=1-nu;
tau1o(j)=nu/d;
%
end
%
numit=51
for j=1:numit
    d=5*10^((j-1)/10);
[nu,fval] = fsolve(@(nu) twop(nu,d),nu0);
% nu
d2(j)=d;
nu2o(j)=1-nu;
tau2o(j)=nu/d;
%
end
%

figure
% plot(tau2,nu2,'ro',x,y,'k-',tau4,nu4,'b-',...
%     tau31,nu3,'g',tau32,nu3,'m', 'LineWidth', 2.00)
% grid
plot(tau2s,nu2s,'ro',x,y,'k-',tau1o,nu1o,'g-',tau2o,nu2o,'g-', 'LineWidth', 2.00)
% grid
% title(['Error diagram, gamma renewal process, \kappa = ',...
%     num2str(1/gmc),...
% ', tmax = ',num2str(tmax)], 'FontSize', fs1)
% % xlabel('\tau', 'FontSize', fs1)
% xlabel('Fraction alarm time, \tau', 'FontSize', fs1)
% % ylabel('\nu', 'FontSize', fs1)
% ylabel('Fraction failures to predict, \nu', 'FontSize', fs1)
% %
text(0.15,.95,['infs = ' num2str(infs),...
    ';  infg = ', num2str(gj), '; sc =',num2str(sc)], 'FontSize', fs1)
text(0.20,.85,['start = ' num2str(ws),...
    ';  range = ', num2str(wd),';  numit = ', num2str(numit)], 'FontSize', fs1)
text(0.3,.75,['scn =',num2str(scn)], 'FontSize', fs1)
axis square
%
axis([0 1 0 1])
%
%
hold on
% [nu,fval] = fsolve(@(nu) twop(nu,d),nu0,options)
%
d=2;
[nu,fval] = fsolve(@(nu) twop(nu,d),nu0);
nu1=[1, 1-nu, 0];
tau1=[0, nu/d, 1];
d=2.5;
[nu,fval] = fsolve(@(nu) twop(nu,d),nu0);
nu2=[1, 1-nu, 0];
tau2=[0, nu/d, 1];
d=3;
[nu,fval] = fsolve(@(nu) twop(nu,d),nu0);
nu3=[1, 1-nu, 0];
tau3=[0, nu/d, 1];
d=4
[nu,fval] = fsolve(@(nu) twop(nu,d),nu0);
nu4=[1, 1-nu, 0];
tau4=[0, nu/d, 1];
d=6
[nu,fval] = fsolve(@(nu) twop(nu,d),nu0);
nu5=[1, 1-nu, 0];
tau5=[0, nu/d, 1];
d=10;
[nu,fval] = fsolve(@(nu) twop(nu,d),nu0);
nu6=[1, 1-nu, 0];
tau6=[0, nu/d, 1];
d=20;
[nu,fval] = fsolve(@(nu) twop(nu,d),nu0);
nu7=[1, 1-nu, 0];
tau7=[0, nu/d, 1];
d=50;
[nu,fval] = fsolve(@(nu) twop(nu,d),nu0);
nu8=[1, 1-nu, 0];
tau8=[0, nu/d, 1];
d=100;
[nu,fval] = fsolve(@(nu) twop(nu,d),nu0);
nu9=[1, 1-nu, 0];
tau9=[0, nu/d, 1];
d=250;
[nu,fval] = fsolve(@(nu) twop(nu,d),nu0);
nu10=[1, 1-nu, 0];
tau10=[0, nu/d, 1];
d=1000;
[nu,fval] = fsolve(@(nu) twop(nu,d),nu0);
nu11=[1, 1-nu, 0];
tau11=[0, nu/d, 1];
d=10000;
[nu,fval] = fsolve(@(nu) twop(nu,d),nu0);
nu12=[1, 1-nu, 0];
tau12=[0, nu/d, 1];
%
%
% figure
% plot(tau2, nu2)
plot(tau1,nu1,'k-',tau2,nu2,'k-',tau3,nu3,'k-',tau4,nu4,'k-',tau5,nu5,'k-',....
tau6,nu6,'k-',tau7,nu7,'k-',tau8,nu8,'k-',tau9,nu9,'k-',tau10,nu10,'k-',tau11,nu11,'k-',...
tau12,nu12,'k-', 'LineWidth', 0.5)
%
% grid
axis square
title('Error diagram, renewal process, information content -- 1 bit', 'FontSize', fs1)
%
xlabel('Fraction alarm time, \tau', 'FontSize', fs1)
% ylabel('\nu', 'FontSize', fs1)
ylabel('Fraction failures to predict, \nu', 'FontSize', fs1)
grid
%
hold off
%
% %
% figure
% %
% xs=2*10^-8
% semilogx(tau2s,nu2s,'ro',x,y,'k-',tau1o,nu1o,'g-',tau2o,nu2o,'g-', 'LineWidth', 2.00)
% text(xs,.25,['infs = ' num2str(infs),...
%     ';  infg = ', num2str(gj), '; sc =',num2str(sc)], 'FontSize', fs1)
% text(xs,.15,['start = ' num2str(ws),...
%     ';  range = ', num2str(wd),';  numit = ', num2str(numit), '; scn =',num2str(scn)],...
%     'FontSize', fs1)
% axis square
% % axis([0 1 0 1])
% %
% % %
% % hold on
% % %
% % semilogx(tau1,nu1,'k-',tau2,nu2,'k-',tau3,nu3,'k-',tau4,nu4,'k-',tau5,nu5,'k-',....
% % tau6,nu6,'k-',tau7,nu7,'k-',tau8,nu8,'k-',tau9,nu9,'k-',tau10,nu10,'k-',tau11,nu11,'k-',...
% % tau12,nu12,'k-', 'LineWidth', 0.5)
% % %
% % % grid
% % axis square
% title('Error diagram, renewal process, information content -- 1 bit', 'FontSize', fs1)
% %
% xlabel('Fraction alarm time, \tau', 'FontSize', fs1)
% % ylabel('\nu', 'FontSize', fs1)
% ylabel('Fraction failures to predict, \nu', 'FontSize', fs1)
% grid
% % %
% % hold off