// Related Software: Scilab
// License     : CeCILL-B : http://www.cecill.info/licences
// Functions   : cummax(), cummin()
// Author (C)  : Samuel GOUGEON - France
// Release     : 1.0
// Date        : 2015-05-25
// Distribution: http://fileexchange.scilab.org/toolboxes/348000
// History     :
//   1.0 - 2015-05-25: first publication on FileExchange
//

function R = cummax(A, orien, fromlast)
//
// Cumulative maxima of values in a vector, matrix or hypermatrix
//
// CALLING SEQUENCES
// cummax           // displays this help
// cummax demo ;    // runs a randomized demo
// R = cummax(A)
// R = cummax(A, orien)
// R = cummax(A, orien, fromlast)
//
// PARAMETERS
// A, R : vector, matrix or hypermatrix of decimals, complex numbers, encoded
//        integers, or booleans. R gets the sizes and the type of A.
//
// orien: "r" | 1 | "c" | 2 | n>2: the index of the dimension across which A
//        must be processed. "r" or 1 means across rows ; "c" or 2 means 
//        across columns ; if A is an hypermatrix  with N dimensions, any 
//        integer N>= n > 2 means across the dimension #n. Default value = 1.
//
// fromlast: optional scalar boolean (default %F). If it is true, the processing 
//        starts from the last element across the chosen orientation, instead 
//        of starting from the first one.
//
// DESCRIPTION
// With R = cummax(A) or R = cummax(A,1) or R = cummax(A,"r"): 
//  For all columns, pages etc, the line #i of R takes the maximal value met
//   across lines #1 to #i of A.
//  If the option fromlast is %T, the line #i of R takes the maximal value met 
//   across lines #i to #$ of A (where $ stands for the index of the last line).
//  If A is boolean, %T is considered as 1 and %F is considered as 0.
//  If A contains complex numbers, numbers with maximal moduli are retained ; 
//   in case of equal moduli, numbers with maximal phases are retained.
//
// With R = cummax(A, "c") or R = cummax(A, 2): A is processed accros columns:
//  For all lines, pages, etc, the column #i of R takes the maximal value met
//   across columns #1 to #i of A.
//
// EXAMPLES
// --------
//// DECIMAL NUMBERS
//-->a = grand(3,5,"uin",0,4)
// a  = [
//    0.    0.    1.    2.    3.  
//    4.    1.    0.    3.    2.  
//    0.    4.    3.    2.    0.  
//    ]
//-->cummax(a)
// ans  = 
//    0.    0.    1.    2.    3.  
//    4.    1.    1.    3.    3.  
//    4.    4.    3.    3.    3.  
//   
//-->cummax(a,2, fromlast=%t)
// ans  =
//    3.    3.    3.    3.    3.  
//    4.    3.    3.    3.    2.  
//    4.    4.    3.    2.    0.  
//
//// ENCODED INTEGERS
//-->a = int8(a)
// a  =
//  0  0  1  2  3  
//  4  1  0  3  2  
//  0  4  3  2  0  
// 
//-->cummax(a,1, fromlast=%t)
// ans  =
//  4  4  3  3  3  
//  4  4  3  3  2  
//  0  4  3  2  0  
//  
//// BOOLEANS
//-->a = a>2
// a  = 
//  F F F F T  
//  T F F T F  
//  F T T F F  
// 
//-->cummax(a)
// ans  =
//  F F F F T  
//  T F F T T  
//  T T T T T  
//
//// COMPLEX NUMBERS
//-->c = grand(3,4,"uin",-2,2) + %i*grand(3,4,"uin",-2,2)
// c  =
//    1. + i      1. - i      1. + i      i         
//    1. + i      0           1. + 2.i    1. + 2.i  
//    1. - i    - 2.        - 1. + i    - 2. - i    
// 
//-->phases = atan(imag(c),real(c))
// phases  =
//    0.7853982  - 0.7853982    0.7853982    1.5707963  
//    0.7853982    0.           1.1071487    1.1071487  
//  - 0.7853982    3.1415927    2.3561945  - 2.677945   
// 
//-->cummax(c)
// ans  =
//    1. + i      1. - i      1. + i      i         
//    1. + i      1. - i      1. + 2.i    1. + 2.i  
//    1. + i    - 2.          1. + 2.i    1. + 2.i  
//
// // HYPERMATRIX:
//-->a = int8(grand(2,4,3,"uin",0,4))
// a  =
//(:,:,1)
//  3  1  2  3  
//  1  2  2  2  
//(:,:,2)
//  4  1  0  3  
//  1  3  1  1  
//(:,:,3)
//  0  4  2  0  
//  4  4  0  3  
// 
//-->cummax(a,3)
// ans  =
//(:,:,1)
//  3  1  2  3  
//  1  2  2  2  
//(:,:,2)
//  4  1  2  3  
//  1  3  2  2  
//(:,:,3)
//  4  4  2  3  
//  4  4  2  3  
//

    // DEMO
    // ----
    if argn(2)==0 then
        head_comments cummax
        R = []
        return
    end
    if argn(2)==1 & typeof(A)=="string" & A=="demo" then
        old = mode()
        mode(6)
        clc, mprintf("a, cummax(a)")
        a = grand(4,3,"uin",0,5)
        cummax(a)
        mprintf("cummax(a, 2)")
        cummax(a, 2)
        mprintf("cummax(a, ""r"", fromlast=%%t)")
        cummax(a, "r", fromlast=%t)
        clc, mprintf("HYPERMATRIX: a, cummax(a,3)")
        a = int8(grand(2,4,3,"uin",0,5))
        cummax(a,3)
        
        clc, mprintf("COMPLEX NUMBERS:")
        c = grand(4,3,"uin",-2,2) + %i* grand(4,3,"uin",-2,2)
        phases = atan(imag(c), real(c))
        mprintf("cummax(c)")
        cummax(c)
        clc, mprintf("c, phases, cummax(c, 2, fromlast=%%t)")
        c=c, phases=phases, cummax(c,2,%t)
        mode(old)
        R = []
        return
    end
    
    // DEFAULT PARAMETERS
    // ------------------
    if argn(2)<2 | orien=="r" then
        orien = 1
    elseif orien=="c"
        orien = 2
    end
    if argn(2)<3 then
        fromlast = %f
    end
    
    imax = size(A,orien)

    if fromlast then 
        i = imax
        irange = imax-1:-1:1
        prev = "i+1"
    else
        i = 1
        irange = 2:imax
        prev = "i-1"
    end

    // BOOLEANS: pre-processing
    bool = type(A)==4
    if bool then
        A = bool2s(A)
    end

    // ENCODED INTEGERS & REALS
    // ------------------------
    if type(A(1))==8 | isreal(A) then
        tmp = repmat(":",1,ndims(A))
        tmp(orien) = "i"
        cn = strcat(tmp, ",")
        imax = size(A,orien)
        R = A
        tmp(orien) = prev
        c = strcat(tmp, ",")
        for i = irange
            execstr("R("+cn+") = max(R("+c+"),A("+cn+"))")
        end 
        
    // COMPLEXES: comparing magnitudes. For equal magnitudes, comparing phases
    // ---------
    else
        sA = size(A)
        if orien>1
            dims = 1:length(sA)
            dims([1 orien]) = [orien 1]
            A = permute(A,dims)
            nsA = size(A)
        else
            nsA = sA
        end
        A = matrix(A, nsA(1), -1)
        R = A
        Xmoduli = abs(A(i,:))
        Xphases = atan(imag(A(i,:)), real(A(i,:)))
        for i = irange
            if fromlast
                p = i+1
            else
                p = i-1
            end
            moduli = abs(A(i,:))
            phases = atan(imag(A(i,:)), real(A(i,:)))
            klt = find( moduli < Xmoduli )
            R(i,klt) = R(p,klt) 
            kgt = find( moduli > Xmoduli )
            keq = find( moduli == Xmoduli )  
            Xmoduli(kgt) = moduli(kgt) 
            Xphases(kgt) = phases(kgt) 
            if keq~=[]
                klt = find( phases(keq) < Xphases(keq) ) 
                R(i,keq(klt)) = R(p,keq(klt)) 
                kgt = find( phases(keq) > Xphases(keq) )
                Xphases(keq(kgt)) = phases(keq(kgt)) 
            end 
        end 
        R = matrix(R, nsA)
        if orien>1
            R = permute(R, dims)
        end
    end
    // BOOLEANS: post-processing
    if bool then
        R = R==1
    end
endfunction

// ===========================================================================

function R = cummin(A, orien, fromlast)
//
// Cumulative minima of values in a vector, matrix or hypermatrix
//
// CALLING SEQUENCES
// cummin           // displays this help
// cummin demo ;    // runs a randomized demo
// R = cummin(A)
// R = cummin(A, orien)
// R = cummin(A, orien, fromlast)
//
// PARAMETERS
// A, R : vector, matrix or hypermatrix of decimals, complex numbers, encoded
//        integers, or booleans. R gets the sizes and the type of A.
//
// orien: "r" | 1 | "c" | 2 | n>2: the index of the dimension across which A
//        must be processed. "r" or 1 means across rows ; "c" or 2 means 
//        across columns ; if A is an hypermatrix  with N dimensions, any 
//        integer N>= n > 2 means across the dimension #n. Default value = 1.
//
// fromlast: optional scalar boolean (default %F). If it is true, the processing 
//        starts from the last element across the chosen orientation, instead 
//        of starting from the first one.
//
// DESCRIPTION
// With R = cummin(A) or R = cummin(A,1) or R = cummin(A,"r"): 
//  For all columns, pages etc, the line #i of R takes the minimal value met
//   across lines #1 to #i of A.
//  If the option fromlast is %T, the line #i of R takes the minimal value met 
//   across lines #i to #$ of A (where $ stands for the index of the last line).
//  If A is boolean, %T is considered as 1 and %F is considered as 0.
//  If A contains complex numbers, numbers with minimal moduli are retained ; 
//   in case of equal moduli, numbers with minimal phases are retained.
//
// With R = cummin(A, "c") or R = cummin(A, 2): A is processed accros columns:
//  For all lines, pages, etc, the column #i of R takes the minimal value met
//   across columns #1 to #i of A.
//
// SEE ALSO: cummax
//
// EXAMPLES
// --------
//-->a = int8(grand(3,4,"uin",0,5))
// a  = [
//  4  3  5  2  
//  2  1  4  4  
//  1  0  5  3  
//  ]
//-->cummin(a)
// ans  =
//  4  3  5  2  
//  2  1  4  2  
//  1  0  4  2  
//  
//-->cummin(a, 2)
// ans  =
//  4  3  3  2  
//  2  1  1  1  
//  1  0  0  0  
// 
//-->cummin(a, 2, fromlast=%t)
// ans  =
//  2  2  2  2  
//  1  1  4  4  
//  0  0  3  3  
// 

    if argn(2)==0 then
        head_comments cummin
        R = []
        return
    end

    // DEMO
    // ----
    if argn(2)==1 & typeof(A)=="string" & A=="demo" then
        old = mode()
        mode(6)
        clc, mprintf("a, cummin(a)")
        a = grand(4,3,"uin",0,5)
        cummin(a)
        mprintf("cummin(a, 2)")
        cummin(a, 2)
        mprintf("cummin(a, ""r"", fromlast=%%t)")
        cummin(a, "r", fromlast=%t)
        clc, mprintf("HYPERMATRIX: a, cummin(a,3)")
        a = int8(grand(2,4,3,"uin",0,5))
        cummin(a,3)
        
        clc, mprintf("COMPLEX NUMBERS:")
        c = grand(4,3,"uin",-2,2) + %i* grand(4,3,"uin",-2,2)
        phases = atan(imag(c), real(c))
        mprintf("cummin(c)")
        cummin(c)
        clc, mprintf("c, phases, cummin(c, 2, fromlast=%%t)")
        c=c, phases=phases, cummin(c,2,%t)
        mode(old)
        R = []
        return
    end
    
    // DEFAULT PARAMETERS
    // ------------------
    if argn(2)<2 | orien=="r" then
        orien = 1
    elseif orien=="c"
        orien = 2
    end
    if argn(2)<3 then
        fromlast = %f
    end
    
    imax = size(A,orien)

    if fromlast then 
        i = imax
        irange = imax-1:-1:1
        prev = "i+1"
    else
        i = 1
        irange = 2:imax
        prev = "i-1"
    end

    // BOOLEANS: pre-processing
    bool = type(A)==4
    if bool then
        A = bool2s(A)
    end

    // ENCODED INTEGERS & REALS
    // ------------------------
    if type(A(1))==8 | isreal(A) then
        tmp = repmat(":",1,ndims(A))
        tmp(orien) = "i"
        cn = strcat(tmp, ",")
        imax = size(A,orien)
        R = A
        tmp(orien) = prev
        c = strcat(tmp, ",")
        for i = irange
            execstr("R("+cn+") = min(R("+c+"),A("+cn+"))")
        end 
        
    // COMPLEXES: comparing magnitudes. For equal magnitudes, comparing phases
    // ---------
    else
        sA = size(A)
        if orien>1
            dims = 1:length(sA)
            dims([1 orien]) = [orien 1]
            A = permute(A,dims)
            nsA = size(A)
        else
            nsA = sA
        end
        A = matrix(A, nsA(1), -1)
        R = A
        Xmoduli = abs(A(i,:))
        Xphases = atan(imag(A(i,:)), real(A(i,:)))
        for i = irange
            if fromlast
                p = i+1
            else
                p = i-1
            end
            moduli = abs(A(i,:))
            phases = atan(imag(A(i,:)), real(A(i,:)))
            kgt = find( moduli > Xmoduli )
            R(i,kgt) = R(p,kgt) 
            klt = find( moduli < Xmoduli )
            keq = find( moduli == Xmoduli )  
            Xmoduli(klt) = moduli(klt) 
            Xphases(klt) = phases(klt) 
            if keq~=[]
                kgt = find( phases(keq) > Xphases(keq) ) 
                R(i,keq(kgt)) = R(p,keq(kgt)) 
                klt = find( phases(keq) < Xphases(keq) )
                Xphases(keq(klt)) = phases(keq(klt)) 
            end 
        end 
        R = matrix(R, nsA)
        if orien>1
            R = permute(R, dims)
        end
    end
    // BOOLEANS: post-processing
    if bool then
        R = R==1
    end
endfunction

// ===========================================================================

// Utility for cummin() and cummax() (for Scilab <6. useless with Scilab 6)
function a = %hm_atan(y,x)
    s = size(y)
    if argn(2)==1
        a = matrix(atan(y(:)), s)
    else
        a = matrix(atan(y(:), x(:)), s)
    end
endfunction
