/*****************************************************************************
*                                                                            *
*   CUMVAR.SPL   Copyright (C) 2022 DSP Development Corporation              *
*                               All Rights Reserved                          *
*                                                                            *
*   Author:      Randy Race                                                  *
*                                                                            *
*   Synopsis:    Cumulative sum                                              *
*                                                                            *
*   Revisions:   19 Jul 2022  RRR  Creation                                  *
*                                                                            *
*****************************************************************************/


#if @HELP_CUMVAR

    CUMVAR

    Purpose: Calculates the cumulative variance of a series.

    Syntax:  CUMVAR(series)

             series        - Any series, multi-series table, or expression
                          resulting in a series or table.

    Returns: A series or table.

    Example:
             cumvar({20, 15, 30, 10, 25})

             returns a new series containing the cumulative
             variance {nan, 12.5, 58.333333, 72.916667, 62.5}.

    Example:
             W1: integ(gnorm(1000, 1))
             W2: cumvar(W1)

             var(W1) == W2[end] 
       
             The last point of CUMVAR is the total var of the input data.

    Remarks:
             CUMVAR calculates the cumulative var of a series.

             The nth value of the output series is equal to the variance
             of the first n points of the input series.

             See CUMAVG for the cumulative average.

    See Also:
             CUMAVG
             CUMMAX
             CUMMIN
             CUMPROD
             VAR
#endif


/* cumulative variance */
cumvar(s, dof, dim, nanmode, method)
{
        local sumx, sumxx, v, nanidx, vidx, N, L;

        /* parse args */
        if (argc < 1 || not(isarray(s)))
        {
                error(sprintf("%s - input series required", __FUNC__));
        }
        else
        {
                (dof, dim, nanmode, method) = stdvar_parse_args(TRUE, dof, dim, nanmode, method);
        }

        if (iscomplex(s))
        {
                v = cumvar(real(s), dof, dim, nanmode, method) + cumvar(imag(s), dof, dim, nanmode, method);

                return(v);
        }

        if (dim > 1) 
        {
                if (dim > 2) 
                {
                        return(s);
                }

                s = s';
        }

        if (isempty(s))
        {
                return(s);
        }

        if (nanmode == "omitnan")
        {
                /* indices of nan */
                nanidx = find(isnan(s));

                if (length(nanidx) > 0)
                {
                        /* variance with nan removed */
                        v = cumvar(removena(s), dof, dim, nanmode, method);

                        /* insert previous values at nan indices */
                        vidx = nanidx - (1..length(nanidx));
                        v = insert(v, v[vidx], vidx, 0);

                        return(v);
                }
        }
        
        s = s - mean(s);

        sumx  = cumsum(s, dim, nanmode, method);
        sumxx = cumsum(s * s, dim, nanmode, method);

        L = length(sumx);
        N = 1..L;

        v = sumxx - (sumx * sumx) / N;

        N = (dof > 0) ? 1..L : 0..(L-1);

        v /= N;

        v[1, ..] = (dof == 0) ? nan : 0.0;

        if (dim > 1)
        {
                v = v';
        }
        
        return(v);
}