%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                                     %
%%           comp_poiss_casc2.m        %
%%                                     %
%%   P. Chainais & R. Riedi & P. Abry  %
%%                                     %
%%              mars 2002              %
%%                                     %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Synthetizes three processes Q_r, A_r and V_r 
% associated to a scale invariant compound Poisson cascade 
% as described in the paper entitled 
%         "Compound Poisson cascades"
% submitted to 
%         "Colloque Auto-Similarite et Applications"
% in Clermont-Ferrand in may 2002
%
% Input :
%			T    : length (duration) of processes,
%			rmin : scale resolution of the cascade,
%			law  : corresponds to distribution F (two examples are implemented here) 
%                    -> 1 : F = Log-Normal distribution,
%                    -> 2 : F = exponential distribution,
%					 -> 3 : F = Dirac (Wi=Cte=W0) = modele log-Poisson / She-Leveque
%			param : parameter of law F.
% Output :
%           Qr   : Poisson cascading noise,
%           Ar   : Poisson cascading motion,
%           Vr   : Poisson cascading Brownian motion,
%           time : time vector for Qr(t), Ar(t) and Vr(t),
%           q    : order for which theoretical exponents are given, (may be changed in the program)
%           TtheoQ : (loosely) "scaling" exponents for Qr,
%           TtheoA : scaling exponents for Ar,
%           TtheoV : scaling exponents for Vr.
%
% Example :
% 
%  T=50 ; rmin=0.0002 ; law=1 ; param=0.06 ;
% [Qr,Ar,Vr,time,q,TtheoQ,TtheoA,TtheoV] = comp_poiss_casc2(T,rmin,law,param) ;
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [Qr,Ar,Vr,time,q,TtheoQ,TtheoA,TtheoV] = comp_poiss_casc2(T,rmin,law,param) ;

if law ~= 3
	C = 1;
else
	C = 2 ;
end

umax = 1/rmin-1 ;          % u will be a uniform variable in [0 1/rmin-1]
lambda = C*(1/rmin-1)*T ; 

%%%%% Time sampling rate = rmin
Dt=rmin/100;

%%%%% Nombre de multiplicateurs et Dt %%%
N = round(lambda) ;         % in theory, should be replaced by a random Poisson variable of mean lambda
% N = poissrnd(lambda,1,1)  % Poisson r.v. of mean lambda
fprintf('Initial number of random multipliers : %d \n',N);

%%%%% (ti,ri) in [0 T]x[rmin 1]
ti = rand(1,N) * T ;
ui = rand(1,N) * umax ;
ri = (1+ui).^(-1) ;

%%%%% Instants for which Qr is computed (1<t<T-1) so that boundary effects are avoided
time = 0:Dt:T-1 ;
ti = ti-1/2 ;   
n=length(time) ;
fprintf('Length of vectors Qr,Ar,Vr : %d \n',n)
di = max(1,ceil((ti-ri/2)/Dt)) ;  %%%%% diviser ri / 2 ????????????
fi = min(n,floor((ti+ri/2)/Dt)) ;
influence = find(di<=fi) ;        % eliminates boundary effects
di = di(influence) ;
fi = fi(influence) ;
Nutil=length(influence) ;
fprintf('Number of multipliers actually used : %d \n',Nutil);

clear influence

%%%%% multipliers Wi
if law == 1 %% F = Log-Normal law
	sigma2 = param ; 
	m = -sigma2/2 ; 
	Wi = exp(randn(Nutil,1)*sqrt(sigma2) + m) ; % such that E[W]=1
end
if law == 2 %% F = Log-exponential law
	Temp = param ;
	Wi = ((1+Temp)^(1/Temp)*rand(Nutil,1)).^Temp ; % power law, E[W]=1
end
if law == 3 %% F = dirac W0 (cf Log Poisson)
	W0 = param ;
	Wi = W0*ones(Nutil,1) ;
end

%%%%%  computation of Qr, Ar and Vr
%% Qr %%
fprintf('Please, be patient...\n')
Qr = ones(1,n) ;
for l=1:Nutil
	Qr(di(l):fi(l)) = Qr(di(l):fi(l)).*ones(1,fi(l)-di(l)+1)*Wi(l) ;
end

%%%%%%%% Normalisation correcte pour cas SI seulement %%%%%%%%%%
if law == 3
	Qr=rmin^(1-W0)*Qr;    % modele log Poisson, W=W0 => normalisation importante
end

%% Ar %%
Ar = cumsum(Qr*Dt) ;

%% Vr %%
Gn=randn(1,length(Qr)) ; 
BmMt= Gn .* sqrt(Qr)*Dt ; 
Vr = cumsum(BmMt) ;

%%%%% Theoretical expressions for scaling exponents
q=1:6;
if law == 1
	sigma2 = param ;
	TtheoQ = -C*(exp(q.*(q-1)*sigma2/2)-1) ;        % Qr
	TtheoA = -C*(exp(q.*(q-1)*sigma2/2)-1) + q ;    % Ar
	H = 1/2 ; qh = H  * q ;
	TtheoV = -C*(exp(qh.*(qh-1)*sigma2/2)-1) + qh  ; % Vr
end

if law == 2 
	Temp = param ;
	TtheoQ = -C*((1+Temp).^q./(1+q*Temp)-1) ;     % Qr 
	TtheoA = -C*((1+Temp).^q./(1+q*Temp)-1) + q ; % Ar
	H = 1/2 ; qh = H  * q ;
	TtheoV = -C*((1+Temp).^qh./(1+qh*Temp)-1) + qh ; % Vr
end

if law == 3
	W0 = param ;
	TtheoQ = C*(-q*(1-W0)+1-W0.^q) ;     % Qr 
	TtheoA = (1-C+C*W0)*q+C*(1-W0.^q); % Ar
	H = 1/2 ; qh = H  * q ;
	TtheoV = (1-C+C*W0)*qh+C*(1-W0.^qh);; % Vr
end