*! version 0.5.3 2014-10-24 | scott long | keep
* version 0.5.2 2014-10-24 | scott long | keep
* version 0.5.1 2014-10-24 | scott long | finish mata update

//  operations are on rows; if matcol, transpose first then tranpose at end

capture program drop _matrcmath
capture program drop _decodeoptions

* DO: see baum on how to retain row.col names

mata
    function mm_select(matsrc, matrows, matcols)
        {
            return(matsrc[matrows,matcols])
        }
end

program _matrcmath

    version 11.2

    syntax , ismatcol(string) matsource(string) type(string) ///
        [ indexout(string) list name(string) NOName ///
          operator(string) indexleft(string)    indexright(string) ///
          isselect(string) selectlist(string)   format(string) ///
          keeplist(string) iskeep(string) ]

//  get matrix and make everything row based

    tempname matcopy rowleft rowright rownew vectemp vecselect ///
        veckeep matsel matcols

    matrix `matcopy' = `matsource' // replace source at end

if "$debug"=="1" matlist `matcopy', title(matcopy is `matcopy')

    * all work is row operations, so if matcol, transpose matrix
    if `ismatcol' matrix `matcopy' = `matcopy'' // transpose
    local rowsN = rowsof(`matcopy')
    local colsN = colsof(`matcopy')
    local origcolnms : colfullnames `matcopy', quoted // reattached later
    local origrownms : rowfullnames `matcopy', quoted // reattached later

//  row/col selection for stat operations

    if `isselect' matrix `vecselect' = (`selectlist')
    else matrix `vecselect' = (.) // will select all

//  augment matrix if added row exceeds current dimensions

    if `indexout'>`rowsN' {
        local newrowsN = `indexout' - `rowsN'
        matrix `matcopy' = `matcopy' \ J(`newrowsN',`colsN',.z)
    }
    if `indexout'==0 {
        matrix `matcopy' = J(1,`colsN',.z) \ `matcopy'
        * since new row, increase all indices
        local indexleft = `indexleft' + 1
        local indexright = `indexright' + 1
        local indexout = `indexout' + 1
    }

if "$debug"=="1" matlist `matcopy', title(after augmentation)

//  binary

    if "`type'"=="bi1" {

        matrix `rowleft' = `matcopy'[`indexleft',1...]
        matrix `rowright' = `matcopy'[`indexright',1...]

        mata: st_matrix(st_local("rownew"), ///
             (st_matrix(st_local("rowleft")) ///
           :`operator' st_matrix(st_local("rowright")) ) )
    }

    if "`type'"=="bi2" {

        matrix `rowleft' = `matcopy'[`indexleft',1...]
        matrix `rowright' = `matcopy'[`indexright',1...]

        if "`operator'"=="%of" { // what % 1st is of 2nd
            mata: st_matrix(st_local("rownew"), ///
             100*(st_matrix(st_local("rowleft")) ///
               :/ st_matrix(st_local("rowright")) ) )
        }
        else if "`operator'"=="%dif" { // % change 1st to 2nd
            mata:  st_matrix(st_local("rownew"), ///
             100*( st_matrix(st_local("rowleft")) ///
                :- st_matrix(st_local("rowright")) ) ///
                :/ st_matrix(st_local("rowleft")) ///
               )
        }
    }

//  unary

    * one step on one row/col
    if "`type'"=="un1" {

        matrix `rowleft' = `matcopy'[`indexleft',1...]
        mata: st_matrix(st_local("rownew"), ///
         `operator'(st_matrix(st_local("rowleft")) ) )
    }

    * multiple steps on one row/col
    if "`type'"=="un2" {

        matrix `rowleft' = `matcopy'[`indexleft',1...]

        if "`operator'"=="sq" {
            mata: st_matrix(st_local("rownew"), ///
                  st_matrix(st_local("rowleft")) ///
               :* st_matrix(st_local("rowleft")) )
        }
        else if "`operator'"=="inv" {
            mata: st_matrix(st_local("rownew"), ///
             1  :/ st_matrix(st_local("rowleft")) )
        }
        else if "`operator'"=="logit" {
            mata:         st_matrix(st_local("rownew"), ///
               log (      st_matrix(st_local("rowleft")) ///
                 :/ ( 1:- st_matrix(st_local("rowleft")) ) ///
               )  )
        }
        else if "`operator'"=="odds" {
            mata:         st_matrix(st_local("rownew"), ///
                          st_matrix(st_local("rowleft")) ///
                 :/ ( 1:- st_matrix(st_local("rowleft")) ) ///
               )
        }

    }
//  utility

    if "`type'"=="util" {

        * add row or column with RC numbers
        if ("`operator'"=="num" | "`operator'"=="number") {
            local nomata = 1
            local matN = colsof(`matcopy')
            numlist "1/`matN'"
            local matnum "`r(numlist)'"
            local matnum : subinstr local matnum " " " , ", all
            matrix `rownew' = ( `matnum' )
            if ("`name'"=="") local name "Number"
        }
    }

//  stats operations on rows

    if "`type'"=="stat" {

        * create matsel with only rows to be summarized
        matrix `matcols' = J(1,1,.) // keep all cols
        mata: st_matrix(st_local("matsel"), ///
        mm_select(st_matrix(st_local("matcopy")), ///
            st_matrix(st_local("vecselect")), ///
            st_matrix(st_local("matcols")) ) )

*matlist `vecselect', title(vecselect)
*matlist `matsel', title(matsel)

        if "`operator'"=="sum" {
            mata: st_matrix(st_local("rownew"), ///
                colsum(st_matrix(st_local("matsel"))) )
        }
        else if "`operator'"=="mean" {
            mata: st_matrix(st_local("rownew"), ///
                mean(st_matrix(st_local("matsel"))) )
        }
        else if "`operator'"=="meanabs" {
            mata: st_matrix(st_local("rownew"), ///
                mean(abs(st_matrix(st_local("matsel")))) )
        }
        else if "`operator'"=="min" {
            mata: st_matrix(st_local("rownew"), ///
                colmin(st_matrix(st_local("matsel"))) )
        }
        else if "`operator'"=="max" {
            mata: st_matrix(st_local("rownew"), ///
                colmax(st_matrix(st_local("matsel"))) )
        }
        else if "`operator'"=="range" {
            mata: st_matrix(st_local("rownew"), ///
                colmax(st_matrix(st_local("matsel"))) ///
                :- colmin(st_matrix(st_local("matsel"))) ///
                )
        }
        else if "`operator'"=="var" | "`operator'"=="sd"{

            * get var first
            mata: st_matrix(st_local("rownew"), ///
                ( /// A
                    colsum( st_matrix(st_local("matsel")) ///
                         :* st_matrix(st_local("matsel")) ) /// colsum
                    :- ///
                    ( ( colsum( st_matrix(st_local("matsel")) ) ///
                        :* colsum( st_matrix(st_local("matsel")) ) ) ///
                      / rows(st_matrix(st_local("matsel"))) ///
                    ) /// B
                ) /// A
                / (rows(st_matrix(st_local("matsel"))) - 1) ///
            )
            if "`operator'"=="sd" {
                mata: st_matrix(st_local("rownew"), ///
                    sqrt(st_matrix(st_local("rownew")) ) )
            }
        } // sd and var

    } // stat

//  name

    if ("`type'"=="name") matrix `rownew' = `matcopy'[`indexout',1...]
    if ("`noname'"=="noname") local name "" // keep original names
    if "`name'"=="" & `indexout'<=`colsN' {
        local name : word `indexout' of `origrownms'
    }
    matrix rowname `rownew' = `"`name'"'

//  create output matrix

if "$debug"=="1" matlist `rownew', title(computed vector)

    capture matrix drop `matsource'
    local rowsN = rowsof(`matcopy')
    forvalues irow = 1/`rowsN' {
        if (`irow'==`indexout') ///
            matrix `matsource' = nullmat(`matsource') \ `rownew'
        else {
            matrix `vectemp' = `matcopy'[`irow',1...]
            matrix `matsource' = nullmat(`matsource') \ `vectemp'
        }
    }
    matrix colnames `matsource' = `origcolnms'

//  keep if specified keep()

    if `iskeep' {

        local rowsN = rowsof(`matsource')
        local colsN = colsof(`matsource')
        local keeplist : subinstr local keeplist "l" "`rowsN'", all
        numlist "`keeplist'"
        matrix `veckeep' = (`r(numlist)')
        matrix `matcols' = J(1,1,.) // keep all cols
*matlist `matsource'        , title(preselect)
        mata: st_matrix(st_local("matsource"), ///
            mm_select(st_matrix(st_local("matsource")), ///
            st_matrix(st_local("veckeep")), ///
            st_matrix(st_local("matcols")) ) )
*matlist `matsource'        , title(post)

    }
    else matrix `veckeep' = (.) // will keep all

//  output

    if `ismatcol' matrix `matsource' = `matsource'' // transpose

    if ("`list'"=="list") ///
        matlist `matsource', title(`matsource') format(`format')

end
exit

*  version 0.3.6 2014-08-08 | scott long | format
*  version 0.3.5 2014-08-08 | scott long | preserve returns
*  version 0.3.4 2014-08-08 | scott long | labels
*  version 0.3.3 2014-08-08 | scott long | abs
*  version 0.3.2 2014-08-04 | scott long | number
*  version 0.3.1 2014-08-03 | scott long |
*  version 0.3.0 2014-08-02 | scott long | tempnames, new features
*  version 0.2.0 2014-08-01 | scott long |
*! version 0.5.0 2014-08-21work | scott long | refactor
* version 0.3.7 2014-08-16 | scott long | missing /// fixed