*! version 0.6.0 2015-04-21 refactoring
* version 0.5.4 2015-04-19 | scott long | remove drop keep

//  matmath: math operations on matrices

// mata function row/col selecting, sorting, or re-ordering
mata
    function mm_sortmat(matsrc, sortcol)
        {
            return(sort(matsrc,sortcol))
        }

    function mm_sortorder(matsrc, sortcol)
        {
            return(order(matsrc,sortcol))
        }
end

capture program drop _matmath
program _matmath, sclass

    version 11.2
    syntax, matout(string) matleft(string) optype(string) ///
    [   ///
        matright(string) operator(string)  ///
        cols(string) rows(string) ///
        value(string) name(string) transpose ///
        rowdec(string) coldec(string) fmttype(string) ///
     ]

// list format(string) droprows(string) dropcols(string) keeprows(string) keepcols(string)

    tempname matrows matcols origreturns
    tempname newvec matinsert sortcol ordervec mat1 mat2
    _return hold `origreturns'
    _return restore `origreturns', hold

if "$debug"=="1" di "matout `matout'"

//  characteristics of the left matrix

    local matleftrownms : rowfullnames `matleft', quoted
    local matleftcolnms : colfullnames `matleft', quoted
    local matleftrowsN = rowsof(`matleft')
    local matleftcolsN = colsof(`matleft')

//  binary: A = B operator C

    if "`optype'"=="OPbinmath" {

        mata: st_matrix(st_local("matout"), ///
              st_matrix(st_local("matleft")) ///
        :`operator' st_matrix(st_local("matright")) )

        matrix colnames `matout' = `matleftcolnms'
        matrix rownames `matout' = `matleftrownms'
    }

//  binary: percent change operations

    else if "`optype'"=="OPbinpct" {

        if "`operator'"=="%of" { // what % is 1st of 2nd matrix
            mata: st_matrix(st_local("matout"), ///
             100*(st_matrix(st_local("matleft")) ///
               :/ st_matrix(st_local("matright")) ) )
        }

        else if "`operator'"=="%dif" { // % change from 1st to 2nd
            mata:  st_matrix(st_local("matout"), ///
             100*( st_matrix(st_local("matleft")) ///
                :- st_matrix(st_local("matright")) ) ///
                :/ st_matrix(st_local("matleft")) ///
               )
        }

        matrix colnames `matout' = `matleftcolnms'
        matrix rownames `matout' = `matleftrownms'
    }

//  binary: appending right or left using mat_c/rapp

    else if "`optype'"=="OPbinconcat" {

        if "`operator'"=="before" ///
            mat_capp `matout' : `matleft' `matright', miss(.)
        else if "`operator'"=="above" ///
            mat_rapp `matout' : `matleft' `matright', miss(.)
    }

//  unary: simple elment by element transformation like sq, log etc.

   else if "`optype'"=="OPunisimple" {

        mata: st_matrix(st_local("matout"), ///
         `operator'(st_matrix(st_local("matleft"))))

        matrix rownames `matout' = `matleftrownms'
        matrix colnames `matout' = `matleftcolnms'
   }

//  unary: these require binary operators (e.g., sq(A))

   else if "`optype'"=="OPunicomplex" {

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

//  unary operators: without math operations (fill and diag);

    else if "`optype'"=="OPbookuni" {

        if "`operator'"=="diag" {
            matrix `matout' = vecdiag(`matleft')
            if ("`name'"=="") local name "diag_`matleft'"
            matrix rowname `matout' = `name' // vec_matname
        }

        if "`operator'"=="fill" {
            matrix `matout' = J(`matleftrowsN',`matleftcolsN',`value')
            matrix colnames `matout' = `matleftcolnms'
            matrix rownames `matout' = `matleftrownms'
        }
    }

//  bookkeeping operations

   else if "`optype'"=="OPbookkeep" {

        mata: `matleft' = st_matrix(st_local("matleft"))

//  select rows/cols to keep

    if "`operator'"=="noomit" { // remove rows and columns with o. omitted
        _noomitted `matleft' row
        _noomitted `matleft' col
        matrix `matout' = `matleft'
    }

    else if "`operator'"=="round" {

        if ("`fmttype'"=="") local fmttype "f"
        matrix `matout' = `matleft'
        local dec `rowdec'
        local isrow = 1
        local iscol = 0
        if "`dec'"=="" {
            local dec `coldec'
            local isrow = 0
            local iscol = 1
        }
        local matoutrowsN = rowsof(`matout')
        local matoutcolsN = colsof(`matout')
        local i = 0
        foreach d in `dec' {
            local ++i
            if `d'<. { // no change if .
                if `isrow' & `i'<=`matoutrowsN' {
                    forvalues c = 1/`matoutcolsN' {
                        local cell = `matout'[`i',`c']
                        local cell = string(`cell',"%16.`d'`fmttype'")
                        matrix `matout'[`i',`c'] = `cell'
                    }
                }
                if `iscol' & `i'<=`matoutcolsN' {
                    forvalues r = 1/`matoutrowsN' {
                        local cell = `matout'[`r',`i']
                        local cell = string(`cell',"%16.`d'`fmttype'")
                        matrix `matout'[`r',`i'] = `cell'
                    }
                }
            } // not missing d
        }
    }

//  insert rows or columns

    else if "`operator'"=="insertafter" {

        if "`rows'"!="" { // insert row

            matrix `newvec' = J(1,`matleftcolsN',`value')
            if ("`name'"!="") matrix rowname `newvec' = `name'

            if (`rows'== 0) matrix `matout' = `newvec' \ `matleft'
            else if `rows'==`matleftrowsN' ///
                matrix `matout' = `matleft' \ `newvec'
            else {
                local st1 = `rows' + 1
                matrix `mat1' = `matleft'[1..`rows',1...]
                matrix `mat2' = `matleft'[`st1'...,1...]
                matrix `matout' = `mat1' \ `newvec' \ `mat2'
            }

        } // insert row

        if "`cols'"!="" { // insert column

            if "`rows'"!="" {
                matrix `matinsert' = `matout' // in case rows were inserted
                local matleftrowsN = rowsof(`matinsert')
            }
            else matrix `matinsert' = `matleft'
            matrix `newvec' = J(`matleftrowsN',1,`value')
            if ("`name'"!="") matrix colname `newvec' = `name'

            if (`cols'==0) matrix `matout' = `newvec' , `matinsert'
            else if `cols'==`matleftcolsN' ///
                matrix `matout' = `matinsert' , `newvec'
            else {
                local st1 = `cols' + 1
                matrix `mat1' = `matinsert'[1...,1..`cols']
                matrix `mat2' = `matinsert'[1...,`st1'...]
                matrix `matout' = `mat1' , `newvec' , `mat2'
            }

        } // insert col

    } // insert

//  sort

    else if "`operator'"=="sort" {

        scalar `sortcol' = `rows'`cols'

        if ("`rows'"!="") local sortbyrow = 1
        else local sortbyrow = 0

        matrix `matout' = `matleft'
        if `sortbyrow' matrix `matout' = `matout''

        * get order of rearranged rows/cols for reassinging labels
        mata: st_matrix(st_local("ordervec"), ///
           mm_sortorder(st_matrix(st_local("matout")), ///
                        st_numscalar(st_local("sortcol"))))

        mata: st_matrix(st_local("matout"), ///
             mm_sortmat(st_matrix(st_local("matout")), ///
                        st_numscalar(st_local("sortcol"))))

        * rearrange labels
        local orderrowsN = rowsof(`ordervec')
        forvalues i = 1/`orderrowsN' {
            local irow = `ordervec'[`i',1]
            if `sortbyrow' local rnm : word `irow' of `matleftcolnms'
            else local rnm : word `irow' of `matleftrownms'
            local newnms `" `newnms' `"`rnm'"' "'
        }
        if `sortbyrow' {
            matrix `matout' = `matout''
            matrix rownames `matout' = `matleftrownms'
            matrix colnames `matout' = `newnms'
        }
        else {
            matrix colnames `matout' = `matleftcolnms'
            matrix rownames `matout' = `newnms '
        }

    } // sort

} // book

end

* version 0.1.0 2015-01-20 from _rm_matrix_noomitted
*! version 0.2.0 2015-04-22 use matsel

//  _noomitted `inputmatrix' row|col

program _noomitted, sclass

    version 11.2

    args inputmat rowcol

    tempname omitmat rerror
    matrix `rerror' = r(error) // from margins if not estimable

    * operations below are column based, so transpose for row selectin
    if ("`rowcol'"=="row") matrix `inputmat' = `inputmat''

    local names : colfullnames `inputmat'
    local rnames : rowfullnames `inputmat'

    mata: inputmat = st_matrix(st_local("inputmat"))
    local colsN = colsof(`inputmat')

    * drop o. in stripe unless not estimable
    _ms_omit_info `inputmat'
    matrix `omitmat' = J(1,`colsN',1) - r(omit)
    forvalues i = 1(1)`colsN' {
        if (`omitmat'[1,`i']==1) local selset "`selset' `i'"
    }
    matsel `inputmat', colkeep(`selset')

    if ("`rowcol'"=="row") matrix `inputmat' = `inputmat''

end

exit

* version 0.4.0 2014-08-18 | scott long | sort
* version 0.5.0 2014-08-21 | scott long | refactor with mata
* version 0.5.1 2014-10-24 | scott long | document code
* version 0.5.4 2015-04-19 | scott long | remove drop keep
* version 0.5.3 2015-04-18 | scott long | adjust to droprows etc change
* version 0.5.2 2014-10-27 | scott long | round; noomit
/*
    * initially select if not o. in stripe
    mata: selectmat = st_matrix(st_local("omitmat"))
    mata: newmat = select(inputmat,selectmat)
    mata: st_matrix(st_local("inputmat"),newmat)
    * names for selected columns
    local names_sel ""
    local i = 0
    foreach nm of local names {
        local ++i
        local keep = `omitmat'[1,`i']
        if `keep' local names_sel `names_sel' `nm'
    }
    matrix colnames `inputmat' = `names_sel'
    matrix rownames `inputmat' = `rnames'
*/

/*
matrix `matcols' = J(1,1,.)
matrix `matrows' = J(1,1,.)
local matleftrowsN = rowsof(`matleft')
local matleftcolsN = colsof(`matleft')
* process keeprows and keepcols
if "`keeprows'"!="" {
    foreach r in `keeprows' {
        if (`r'>`matleftrowsN') | (`r'<1) {
            di as error "row `r' in keeprows() is invalid"
            local iserror = 1
        }
    }
    local tmp = subinstr("`keeprows'"," ","\",.)
    matrix `matrows' = (`tmp')
}
if "`keepcols'"!="" {
    foreach c in `keepcols' {
        if (`c'>`matleftcolsN') | (`c'<1) {
            di as error "column `c' in keepcolumns() is invalid"
            local iserror = 1
        }
    }
    local tmp = subinstr("`cols'"," ","\",.)
    matrix `matcols' = (`tmp')
}
if "`iserror'"=="1" {
    _return restore `origreturns'
    exit
}
mata: st_matrix(st_local("matout"), ///
    mm_select(st_matrix(st_local("matleft")), ///
           st_matrix(st_local("matrows")), ///
           st_matrix(st_local("matcols")) ) )
    if ("`transpose'"=="transpose") matrix `matout' = `matout''

    if ("`list'"=="list") ///
        matlist `matout', title(matrix `matout') format(`format')
*/


//  reattach labels for drop and keep

        if "`operator'"=="keep" | "`operator'"=="drop" {

            if "`rows'"!="" {
                foreach row in `rows' {
                    local rnm : word `row' of `matleftrownms'
                    * `"` '"' needed since names have ""s in them
                    local newrownms `" `newrownms' `"`rnm'"' "'
                }
                matrix rownames `matout' = `newrownms'
            }
            else matrix rownames `matout' = `matleftrownms'

            if "`cols'"!="" { // selected columns
                foreach col in `cols' {
                    local cnm : word `col' of `matleftcolnms'
                    local newcolnms `" `newcolnms' `"`cnm'"' "'
                }
                matrix colnames `matout' = `newcolnms'
            }
            else matrix colnames `matout' = `matleftcolnms'
        }
 * use matsel now...

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