/****************************************************************************
*                                                                           *
*   STD.SPL     Copyright 2007, 2023 (C) DSP Development Corporation        *
*                                                                           *
*   Author:     Randy Race                                                  *
*                                                                           *
*   Synopsis:   Column based standard deviation                             *
*                                                                           *
*   Revisions:  13 Aug 2007     RRR     Creation                            *
*                8 Feb 2023     RRR     optional args                       *
*                                                                           *
****************************************************************************/

#if @HELP_STD

    STD

    Purpose: Calculates the standard deviation of each column of a series

    Syntax:  STD(s, mode, dim)

                 s - input series

              mode - Optional integer. The normalization mode.
                     For N = length(s)

                       0: normalize by 1/(N-1) (default)
                       1: normalize by 1/N

               dim - Optional integer. Summation dimension.

                      1: calculate row-wise (default)
                      2: calculate column-wise

    Returns: A real scalar for a one column series or a 1xM real table for
             an M column series.

    Example:
             W1: 1..10
             std(w1)

             Returns 3.0277, the sample standard deviation.

    Example:
             W1: 1..10
             std(w1+1e7)

             Returns 3.0277, same as above since the true standard deviation
             is independent of the mean value.

    Example:
             W1: 1..10
             std(w1, 1)

             Returns 2.8723, the population standard deviation.

    Example:
             W1: ravel(1..9, 3)^2
             W2: std(w1)

             Returns the series {{4.0415, 10.0167, 16.0104}}

             the standard deviation of each column.

    Remarks:
             If s is the input series and N the length of the series,
             the basic form of the sample standard deviation (mode == 0)
             is:

                std(s) = sqrt(sum((s - mean(s))^2/(N-1)))

             The population standard deviation (mode == 1) is defined as:

                std(s) = sqrt(sum((s - mean(s))^2/N))


             STD uses a fast, highly accurate running sums algorithm that
             exhibits insensitivity to round-off errors.

             As shown in the second example, the standard deviation is
             independent of an additive mean value.

             STD returns the standard deviation of each column for a
             multi-column series.  See STDEV to return the standard
             deviation of a multi-column series as a whole.

    See Also:
             Mean
             Stdev
#endif


/* standard deviation of each column */
std(argv)
{
        local N, sd, mn;

        (s, mode, dim, flags) = std_parse_args(argv);

        /* check if stdev on full array */
        if ((isstring(dim)   && (tolower(dim)   == "all")) || 
            (isstring(flags) && (tolower(flags) == "all")) || 
            (isscalar(dim)   && dim < 0))
        {
                dim = 0;
        }
        else if (dim > 1)
        {
                s = s';
        }

        N = length(s);

        if (dim == 0)
        {
                /* full array */
                (sd, mn) = stdev(s, -1, -1, mode, flags);
        }
        else
        {
                /* stdev of each column or row */
                (sd, mn) = std_calc(s, mode, flags);
        }

        if (dim > 1)
        {
                sd = sd';
        }
        
        return(sd, mn);
}


/* calculate stdev of each column */
ITERATE std_calc(s, mode, flags)
{
        local sd, mn;

        (sd, mn) = stdev(s, -1, -1, mode, flags);
        
        return(sd, mn);
}


/* get args */
std_parse_args(argv)
{
        local s, mode, dim, flags;

        s = 0;
        mode = dim = flags = {};

        loop (j = 1..argc)
        {
                if (isarray(getargv(j)))
                {
                        s = refseries(getargv(j));
                }
                else if (isscalar(getargv(j)))
                {
                        if (isempty(mode))
                        {
                                mode = getargv(j);
                        }
                        else if (isempty(dim))
                        {
                                dim = getargv(j);
                        }
                }
                else if (isstring(getargv(j)))
                {
                        if (isempty(flags))
                        {
                                flags = getargv(j);
                        }
                        else if (tolower(getargv(j)) == "all")
                        {
                                dim = 0;
                        }
                }
        }

        s     = isscalar(s)    ? refseries(w0) : refseries(s);
        mode  = isempty(mode)  ? 0  : mode;
        dim   = isempty(dim)   ? 1  : dim;
        flags = isempty(flags) ? "" : flags;

        if (not(isscalar(mode)) || not(mode == 0 || mode == 1))
        {
                error(sprintf("%s - mode must be 0 or 1", __CALLER__));
        }

        return(s, mode, dim, flags);
}