%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                                     %
%%           warpIDCnoise.m            %
%%                                     %
%%   P. Chainais & R. Riedi & P. Abry  %
%%                                     %
%%              mars 2004              %
%%                                     %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Synthesizes infinitely divisible cascading noise/density Q_r(t), 
% associated to a multifractal (power law scaling) Compound Poisson Cascade or Log Normal cascade,
%
%    >>>>>>>>>>>       WITH WARPED CONTROL MEASURE                <<<<<<<<<<<<
% 
% See references : ( www.isima.fr/~chainais/PUB/publications.html ) 
%
% On non scale invariant infinitely divisible cascades 
% LIMOS/RR-04-06 - Auteurs : P. Chainais1, R. Riedi, P. Abry . 
% at    http://www.isima.fr/limos/publications.htm 
%
% On non scale invariant infinitely divisible cascades ,
% P. Chainais, R. Riedi, et P. Abry,Ê in Actes du GRETSI, Paris, 2003. 
% 
% Scale invariant infinitely divisible  cascades ,
% P. Chainais, R. Riedi, et P. Abry,Ê in Proceedings of  PSIP'03, Grenoble, 2003. 
% 
% Compound  Poisson Cascades ,
% P. Chainais, R. Riedi, et P. Abry,Ê Annales de l'UniversitŽ Blaise Pascal 
% Workshop Self-Similarity and Applications, Clermont-Ferrand 2002. 
%
%%%%%%%%%%%%%%%%%
% Input :
%			T    : length (time duration) of processes,
%			rmin : resolution of the cascade,
%           Dt   : time sampling rate,
%			Law  : corresponds to the chosen infinitely divible distribution.
%                  four examples are implemented here : 0               is a Log Normal cascade,
%                                                       others (1 to 3) are compound Poisson cascades 
%                                                                           with distribution F.
%           ParamLaw : parameters of the inf. div. law
%                    -> Law : (ParamLaw)
%                    ->   0 : Log-Normal cascade [sigma2],
%                    ->   1 : F = Log-Normal distribution [mu,sigma2],
%                    ->   2 : F = exponential distribution [Temp], 
%                    ->   3 : F = Dirac (Wi=Cte=W0) (log-Poisson ou She-Leveque) [W0], (W0<1)
%			C : multiplicative constant [default is C=1],
%           beta : parameters of the control measure dm(t,r)=dt.dr/r^(2+beta).
%
%%%%%%%%%%%%%%%%%
% Output :
%           Qr   : IDC noise/density,
%           time : time vector for Qr(t),
%           warpfun : warping term, of the form m(cone_tau) =
%           (1-tau^(-beta))/(-beta)
%           q    : order for which theoretical exponents are given, (may be changed in the program)
%           TtheoQ : theoretical scaling exponents for Qr, q + C.phi(q)
%
%%%%%%%%%%%
% Example 1 (CPC) :
%-------------------
%  T=50 ; rmin=0.002 ; Dt=2*rmin ; Law=3 ; ParamLaw=2/3 ; beta=-0.3; C=1 ; 
%  [Qr,time,TtheoQ,q,warpfun,scales] = warpIDCnoise(T,rmin,Dt,Law,ParamLaw,beta,C,q) ;
%
% Example 2 (Log Normal) :
%-------------------------
%  T=10 ; rmin=0.02 ; Dt=rmin/2 ; Law=0 ; ParamLaw=0.2 ; beta=-0.4; C=1 ; q=(0:8)/2 ; 
%  [Qr,time,TtheoQ,q,warpfun,scales] = warpIDCnoise(T,rmin,Dt,Law,ParamLaw,beta,C,q) ;
%
% >>>>   ! !  N.B.: if Law = 0, (rmin/Dt) and (1/Dt) must be integers  ! !  <<<<<<<<
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [Qr,time,TtheoQ,q,warpfun,scales] = warpIDCnoise(T,rmin,Dt,Law,ParamLaw,beta,C,q) ;

%%%% input tests
if nargin<6 error('Some input arguments are missing...') ;
elseif nargin<7
    C = 1 ;
    if nargin<8
        q = 1:6 ;
    end
end

% if Law>3 error('Law must be 0, 1, 2 or 3 !');
        
%%%% Warping term m(C_tau)
%scales=logspace(log10(rmin/5),1,100);
scales=logspace(log10(Dt),0,100);
warpfun = -log(scales) ;

if beta ~= 0
	warpfun = (1-(scales).^(-beta))/(-beta) ;
end

%%%%%%%%%%%%%%%%%%%%%%% Initialisation   %%%%%%%%%%%%%%%%%%%%%%%
time = 0:Dt:T ;
n=length(time) ;
fprintf('Length of vector Qr: %d \n',n)

            %%%%%%%% Compound Poisson Cascades %%%%%%%%%%

if Law  % Law > 0 =>  Compound Poisson Cascade based on distribution F 
    Tmax = T + 1 ;             % because of border effects
 	umax = rmin^(-beta-1)-1 ;        % u will be a uniform variable in [0 1/rmin^(1+beta)-1]

    %%%%% Number of multipliers %%%
	lambda = C/(1+beta)*(1/rmin^(1+beta)-1)*Tmax ; % warped cascade with dm(t,r)=Cdrdt/r^(2+beta)
    N = round(lambda) ;        % Number of multipliers, 
                               % 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 Tmax]x[rmin 1]
	ti = -1/2 + rand(1,N) * Tmax ;
    ui = rand(1,N) * umax ;
    ri = (1+ui).^(-1/(beta+1)) ;
    
    di = max(1,ceil((ti-ri/2)/Dt)) ;  % definition of determining Wi for Qr(t)
    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

    %%%%% Throwing multipliers Wi
    switch Law
        case 1, % F = Log-Normal Law
            mu = ParamLaw(1) ;
            sigma2 = ParamLaw(2) ;
            H1 = C * (1 - exp(mu+sigma2/2)) ;  % normalization \E Qr = 1
	        TtheoQ = -q*H1 + C * (1 - exp(mu*q + q.^2*sigma2/2)) ;
            Wi = exp(randn(Nutil,1)*sqrt(sigma2) + mu) ; 
            nondegenere = max(0,1-H1-C*(mu+sigma2)*exp(mu+sigma2/2)) ; 
       case 2, % F = Log-exponential Law
    	    Temp = ParamLaw(1) ;
            H1 = C*0 ;
        	TtheoQ = -C*((1+Temp).^q./(1+q*Temp)-1) ;     % Qr 
	        Wi = ((1+Temp)^(1/Temp)*rand(Nutil,1)).^Temp ; % power law, E[W]=1
            nondegenere = max(0,1+C/(1+T)*(T-(1+T)*log(1+T))) ; 
        case 3, % F = dirac W0 (cf Log Poisson)
	        W0 = ParamLaw(1) ; 
            H1 = C*(1-W0) ;
        	TtheoQ = C*(-q*(1-W0)+1-W0.^q) ;     % Qr 
	        Wi = W0*ones(Nutil,1) ;
            nondegenere = max(0,1+C*((W0-1)-W0*log(W0))) ; 
        otherwise,
            error('Error : Law must be 0, 1, 2 or 3 !\n')
    end
    
    %%%%% Test of degeneracy %%%%%%
    if ((~beta) & (~nondegenere)) fprintf('> > > > >>> Be careful : the measure degenerates ! [ zeta''(1) < -1 ] \n')
    end

    %%%%%  Computation of Qr(t) %%%%
    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  %%%%%%%% 
    if (~beta)
        % normalisation pour cas SI 
        Qr = Qr / rmin^H1 ;  
    else
        % normalisation pour cas non SI, 1/r^(2+beta)
        Qr = Qr * exp( (rmin^(-beta)-1)/beta * H1 ) ;
    end

    
else               % Law = 0 => Log Normal cascade

    %%%%%%%%%%%%%%%%%%%%  LOG NORMAL  CASCADE  %%%%%%%%%%%%%%%%%%
    
    if (((rmin/Dt)-fix(rmin/Dt)) | ((1/Dt)-fix(1/Dt)))
        error('rmin/Dt and 1/Dt must be integers in that case !') ;
    end
    rho0 = round(rmin/Dt) ; % rho0 : rmin/Dt is an integer, scale resolution of the cascade,
    rho1 = round(1/Dt) ;    % rho1 : 1/Dt is an integer, time resolution,

    mu = -ParamLaw(1)/2 ;
    sigma2 = ParamLaw(1) ;
    
    %%%%% Theoretical expressions for scaling exponents
    TtheoQ = -mu*q-sigma2/2*q.^2 ;        % Qr

    %%%%% Test of degeneracy %%%%%%
    nondegenere = max(0,1-sigma2/2) ;
    if (~beta & ~nondegenere) fprintf('> > > > >>> Be careful : the measure degenerates ! [ zeta''(1) < -1 ] \n')
    end

    %%%%% Control measures for alpha, beta and gamma "multipliers" %%%
    if ~beta
        %% alpha : triangles "tete en bas", sous r=1
        mes1 = C*(log(rho1/(rho1-1))-1/rho1) ;
    
        %% beta : losanges "interieurs"
        j = (rho1-1):-1:(rho0+1) ; 
        mes2 = C*(log(j.^2./((j-1).*(j+1)))) ;
    
        %% gamma : triangles "tete en haut", sur r=rmin
        mes3 = C*(1/rho0-log((rho0+1)/rho0)) ;
    else
        %% alpha : triangles "tete en bas", sous r=1
        rbot=(rho1-1)*Dt;
        rtop=rho1*Dt;
        mes1 = C*(rtop^(-beta)*(rbot/rtop/(1+beta)-1/beta)-rbot^(-beta)*(1/(1+beta)-1/beta)) ;
    
        %% beta : losanges "interieurs"
        j = (rho1-1):-1:(rho0+1) ;
        rtop=(j+1)*Dt;
        rmid=j*Dt;
        rbot=(j-1)*Dt;
        mes2 = C*(1/(1+beta)-1/beta)*(2*rmid.^(-beta)-rtop.^(-beta)-rbot.^(-beta)) ;
        clear rbot rtop
    
        %% gamma : triangles "tete en haut", sur r=rmin
        rbot=rho0*Dt;
        rtop=(rho0+1)*Dt;
        mes3 = C*(rbot^(-beta)*(rtop/rbot/(1+beta)-1/beta)-rtop^(-beta)*(1/(1+beta)-1/beta)) ;
    
        clear rbot rtop rmid
    end


    %%%%% Initialization of memory matrix MatMem
    MatMem = zeros(rho1,rho1) ;

    %% alpha = first line
    MatMem(1,:) = dilatrandn(mu*mes1,sigma2*mes1,1,rho1) ;

    %% beta = "inside of MatMem"
    for jj=2:(rho1-rho0)
	    MatMem(jj,jj:rho1) = dilatrandn(mu*mes2(jj-1),sigma2*mes2(jj-1),1,rho1-jj+1) ;
    end

    %% gamma = last line
    MatMem(rho1-rho0+1,rho1-rho0+1:rho1) = dilatrandn(mu*mes3,sigma2*mes3,1,rho0) ;


    %%%%% Initialization of logQr(0)   <<<<<<<<<<<<<<<<<<<<<
    logQr(1) = sum(sum(MatMem)) ;

    %%%%% Random set of values of last column LastColumn for the whole
    %%%%% realization
    LastColumn=zeros(rho1,1) ;  

    muvec = [mu*mes1, mu*mes2, mu*mes3];
    sigvec = [sigma2*mes1,sigma2*mes2,sigma2*mes3];

    %%%%% Indices 
    ii = 1:rho1 ;
    jj = 1:rho1 ;
    [J,I] = meshgrid(ii,jj) ;


    for k=2:(rho1*T+1)
        %    pause
        % 	if ~mod(k,1000) 
        % 		fprintf(' %d',k) 
        %     end
	
        %% computes logQr(k)
        LastColumn(1:rho1-rho0+1) = dilatrandn(muvec,sigvec);
        tempmat = MatMem(find(J>I)) ;
        logQr(k) = sum(tempmat) + sum(LastColumn) ;
    
	    %% update MatMem :
        MatMem(find((J>=I) & (J<rho1))) = tempmat ;
	    MatMem(:,end) = LastColumn ;
    end

    %%%% Qr = exp(M(cone))
    Qr = exp(logQr) ;

    %%%% Normalisation
    H1 = -mu-sigma2/2 ;
    if ~beta
        Qr = (rho1/rho0)^H1*exp(logQr) ;
    else
        Qr = exp(H1*(1-(rho0/rho1)^(-beta))/(-beta))*exp(logQr) ;
    end

end
    
    
