*! enforce v1.0 Thomas Blanchet 20apr2020

// -------------------------------------------------------------------------- //
// Command to simultaneously enforce several accounting identities on
// several variables
// -------------------------------------------------------------------------- //

program enforce
    version 13

    // Check if gtools installed
    cap which gtools
    if (_rc == 1) {
        di as error "gtools required; run {stata ssc install gtools}"
        exit 198
    }

    syntax anything [if] [in], [FIXedvars(varlist) NOFILLmissing NOENForce ///
        NOCHeck REPlace SUFfix(string) PREfix(string) TOLerance(real 1) zero(real 1e-7) DIAGnostic force]
    marksample touse

    // ---------------------------------------------------------------------- //
    // Parse accounting identities
    // ---------------------------------------------------------------------- //

    if ("`anything'" == "") {
        di as err "no identity specified"
        exit 198
    }

    local i = 1 // Identity number

    // Loop over identities
    gettoken expr other: anything, match(parns) bind
    while ("`expr'" != "") {
        if ("`parns'" == "") {
            di as error "identities must be specified in parentheses"
            exit 198
        }

        local idenexpr`i' `expr'

        // Separate LHS from RHS
        tokenize "`expr'", parse("=")
        local lhs = "`1'"
        local rhs = "`3'"

        if ("`rhs'" == "") {
            di as err "expression (`expr') has no equal sign"
            exit 198
        }

        // Parse LHS
        gettoken left right: lhs, parse("+- ")
        while ("`left'" != "") {
            if ("`left'" == "+") {
                // Include with positive sign
                gettoken var right: right, parse("+- ")

                if ("`var'" == "0") {
                    gettoken left right: right, parse("+- ")
                    continue
                }

                local idensign`i' `idensign`i'' 1
                confirm numeric variable `var'
                local idenvars`i' `idenvars`i'' `var'
                local allvars `allvars' `var'

                gettoken left right: right, parse("+- ")
            }
            else if ("`left'" == "-") {
                // Include with negative sign
                gettoken var right: right, parse("+- ")

                if ("`var'" == "0") {
                    gettoken left right: right, parse("+- ")
                    continue
                }

                local idensign`i' `idensign`i'' -1
                confirm numeric variable `var'
                local idenvars`i' `idenvars`i'' `var'
                local allvars `allvars' `var'

                gettoken left right: right, parse("+- ")
            }
            else {
                // Include with (implicit) positive sign
                local var = "`left'"

                if ("`var'" == "0") {
                    gettoken left right: right, parse("+- ")
                    continue
                }

                local idensign`i' `idensign`i'' 1
                confirm numeric variable `var'
                local idenvars`i' `idenvars`i'' `var'
                local allvars `allvars' `var'

                gettoken left right: right, parse("+- ")
            }
        }

        // Parse RHS in the same way, except the signs are reversed
        gettoken left right: rhs, parse("+- ")
        while ("`left'" != "") {
            if ("`left'" == "+") {
                // Include with positive sign
                gettoken var right: right, parse("+- ")

                if ("`var'" == "0") {
                    gettoken left right: right, parse("+- ")
                    continue
                }

                local idensign`i' `idensign`i'' -1
                confirm numeric variable `var'
                local idenvars`i' `idenvars`i'' `var'
                local allvars `allvars' `var'

                gettoken left right: right, parse("+- ")
            }
            else if ("`left'" == "-") {
                // Include with negative sign
                gettoken var right: right, parse("+- ")

                if ("`var'" == "0") {
                    gettoken left right: right, parse("+- ")
                    continue
                }

                local idensign`i' `idensign`i'' 1
                confirm numeric variable `var'
                local idenvars`i' `idenvars`i'' `var'
                local allvars `allvars' `var'

                gettoken left right: right, parse("+- ")
            }
            else {
                // Include with (implicit) positive sign
                local var = "`left'"

                if ("`var'" == "0") {
                    gettoken left right: right, parse("+- ")
                    continue
                }

                local idensign`i' `idensign`i'' -1
                confirm numeric variable `var'
                local idenvars`i' `idenvars`i'' `var'
                local allvars `allvars' `var'

                gettoken left right: right, parse("+- ")
            }
        }

        gettoken expr other: other, match(parns) bind
        local i = `i' + 1
    }

    local niden = `i' - 1

    // List of (unique) variables involved in the problem
    local uniqvars: list uniq allvars
    local uniqvars: list sort uniqvars

    // Check if there fixed vars not involved in the problem (ignore with a warning)
    local ignoredvars: list fixedvars - allvars
    foreach v of local ignoredvars {
        di as text "note: variable {bf:`v'} not included in any identity; ignored"
    }

    local fixedvars: list uniq fixedvars
    local fixedvars: list fixedvars & allvars
    local fixedvars: list sort fixedvars

    // List of non-fixed variable
    local nonfixedvars: list uniqvars - fixedvars
    local nonfixedvars: list uniq nonfixedvars
    local nonfixedvars: list sort nonfixedvars

    local nuniqvars:     word count `uniqvars'
    local nfixedvars:    word count `fixedvars'
    local nnonfixedvars: word count `nonfixedvars'

    if (`nfixedvars' > 0) {
        confirm numeric variable `fixedvars'
    }

    if (`nnonfixedvars' == 0) {
        di as err "no nonfixed variables in the system"
        exit 412
    }

    // Create temporary version of variables
    foreach v of varlist `uniqvars' {
        tempvar tmp`v'
        qui generate `tmp`v'' = `v', after(`v')
        label variable `tmp`v'' "`v'"
        local tmpuniqvars `tmpuniqvars' `tmp`v''
    }

    // ---------------------------------------------------------------------- //
    // Check that variables to be generated do not already exist
    // ---------------------------------------------------------------------- //

    if ("`suffix'" != "") {
        foreach v of local uniqvars {
            cap confirm existence `v'
            if (_rc != 0) {
                di as err "{bf:`v'} already exists"
                exit 198
            }
        }

    }
    else if ("`prefix'" != "") {
        foreach v of local uniqvars {
            cap confirm existence `v'
            if (_rc != 0) {
                di as err "{bf:`v'} already exists"
                exit 198
            }
        }
    }
    else {
        if ("`replace'" == "") {
            di as error "{bf:replace}, {bf:suffix(...)}, or {bf:prefix(...)} required"
            exit 198
        }
    }

    if ("`diagnostic'" != "" & "nocheck" == "") {
        cap confirm existence _check
        if (_rc != 0) {
            di as err "{bf:`v'} already exists"
            exit 198
        }
    }

    // ---------------------------------------------------------------------- //
    // Construct matrix of the identities
    // ---------------------------------------------------------------------- //

    tempname matiden
    matrix define `matiden' = J(`niden', `nuniqvars', 0)

    forvalues i = 1/`niden' {
        local n: word count `idenvars`i''
        forvalues j = 1/`n' {
            local v: word `j' of `idenvars`i''
            local s: word `j' of `idensign`i''

            local pos : list posof "`v'" in uniqvars
            matrix `matiden'[`i', `pos'] = `s'
        }
    }

    matrix colnames `matiden' = `uniqvars'

    // Indexes of fixed variables
    tempname fixedvaridx
    matrix define `fixedvaridx' = J(1, `nuniqvars', 0)
    matrix colnames `fixedvaridx' = `uniqvars'
    foreach v of local fixedvars {
        matrix `fixedvaridx'[1, colnumb(`fixedvaridx', "`v'")] = 1
    }

    // ---------------------------------------------------------------------- //
    // Check the identities
    // ---------------------------------------------------------------------- //

    if ("`nocheck'" == "") {
        tempname check_sys check_zer
        qui generate `check_sys' = 0 if `touse'
        qui generate `check_zer' = 0 if `touse'
        label variable `check_sys' "check_sys"
        label variable `check_zer' "check_zer"

        mata: check_identities()

        di ""
        di as text "Data diagnostics"
        di "{hline 50}{c TT}{hline 20}"
        di "{lalign 49: type of issue} {c |}{ralign 20: # of obs.}"
        di "{hline 50}{c +}{hline 20}"

        di as text "{lalign 49: no issue detected} {c |}",, _continue
        qui count if `check_sys' == 0 & `check_zer' == 0 & `touse'
        di as text "{space 3}" as result %17.0g = r(N),, _continue
        if ("`diagnostic'" == "") {
            di as text ""
        }
        else {
            di as text "{space 2}(marked as {bf:_check = 0})"
        }

        di as text "{lalign 49: no solution to the system of identities} {c |}",, _continue
        qui count if `check_sys' == 1 & `touse'
        di as text "{space 3}" as result %17.0g = r(N),, _continue
        if ("`diagnostic'" == "") {
            di as text ""
        }
        else {
            di as text "{space 2}(marked as {bf:_check = 1})"
        }

        di as text "{lalign 49: no solution because of variables equal to zero} {c |}",, _continue
        qui count if `check_zer' == 1 & `touse'
        di as text "{space 3}" as result %17.0g = r(N),, _continue
        if ("`diagnostic'" == "") {
            di as text ""
        }
        else {
            di as text "{space 2}(marked as {bf:_check = 2})"
        }

        di "{hline 50}{c BT}{hline 20}"
        if ("`diagnostic'" == "") {
            di as text "Use option {bf:diagnostic} to identify problematic observations."
        }
        else {
            qui generate _check = 0 if `touse'
            qui replace _check = 1 if `touse' & `check_sys' == 1
            qui replace _check = 2 if `touse' & `check_zer' == 1

            qui label define _check 0 "ok" 1 "no solution" 2 "no solution due to zeros", replace
            qui label values _check _check
        }

        qui count if (`check_sys' == 1 | `check_zer' == 1) & `touse'
        if (r(N) > 0) {
            di as text ""
            di as text "The command detected issues with your data that make the system"
            di as text "of identities effectively unsolvable."
            if ("`force'" != "") {
                di as text ""
                di as text "You used the option {bf:force}, so the command will proceed"
                di as text "anyway and try to find the best solution in spite of"
                di as text "these issues."
            }
            else {
                di as text ""
                di as text "The issues may be false alerts due to rounding errors, in which"
                di as text "you can increase tolerance via the {bf:tolerance} option, or use"
                di as text "the options {bf:force} or {bf:nocheck} to proceed anyway. But be sure to"
                di as text "double-check your data or identities before you continue."
                exit 412
            }
        }
    }

    // ---------------------------------------------------------------------- //
    // Fill missing values
    // ---------------------------------------------------------------------- //

    if ("`nofillmissing'" == "") {
        // Group data by block with identical missing variable structure
        foreach v of varlist `tmpuniqvars' {
            tempvar miss`v'
            qui generate `miss`v'' = missing(`v') if `touse'
            local missvars `missvars' `miss`v''
        }
        tempvar missgroup
        qui gegen `missgroup' = group(`missvars') if `touse'

        qui glevelsof `missgroup' if `touse', local(missgroups)

        foreach group of local missgroups {
            tempname gind
            qui generate `gind' = (`missgroup' == `group') & `touse'

            mata: fill_missing()

            drop `gind'
        }

        drop `missgroup' `missvars'

        // Show report of missing values filled
        local colwidth   = 12
        local dbcolwidth = 2*`colwidth'
        local width2     = 4*`dbcolwidth' + 1
        local width3     = round(`colwidth'/2)
        local width4     = `dbcolwidth' - `width3'
        local width5     = `colwidth' - 3

        di ""
        di as text "Missing values recovered"
        di as text "{hline 13}{c TT}{hline `width2'}"

        di as text "{ralign 12:} {c |}", _continue
        di as text "{space `width3'}{rcenter `width4': nonmissing}{space `width3'}{rcenter `width4': mean}",, _continue
        di as text "{space `width3'}{rcenter `width4': min}{space `width3'}{rcenter `width4': max}"

        di as text "{ralign 12:Variable} {c |}", _continue
        di as text 4*"{ralign `colwidth': before}{ralign `colwidth': after}"

        di as text "{hline 13}{c +}{hline `width2'}"

        foreach v of varlist `uniqvars' {
            summarize `v', meanonly
            local nmiss_before = r(N)
            local mean_before  = r(mean)
            local min_before   = r(min)
            local max_before   = r(max)

            summarize `tmp`v'', meanonly
            local nmiss_after = r(N)
            local mean_after  = r(mean)
            local min_after   = r(min)
            local max_after   = r(max)

            local abbrvar = abbrev(`"`v'"', 12)

            di as text "{ralign 12:`abbrvar'} {c |}", _continue
            di as text "{space 3}" as result %`width5'.0g = `nmiss_before',, _continue
            di as text "{space 3}" as result %`width5'.0g = `nmiss_after',, _continue
            di as text "{space 3}" as result %`width5'.3g = `mean_before',, _continue
            di as text "{space 3}" as result %`width5'.3g = `mean_after',, _continue
            di as text "{space 3}" as result %`width5'.3g = `min_before',, _continue
            di as text "{space 3}" as result %`width5'.3g = `min_after',, _continue
            di as text "{space 3}" as result %`width5'.3g = `max_before',, _continue
            di as text "{space 3}" as result %`width5'.3g = `max_after'
        }
        di as text "{hline 13}{c BT}{hline `width2'}"
    }

    // ---------------------------------------------------------------------- //
    // Enforce identities
    // ---------------------------------------------------------------------- //

    if ("`noenforce'" == "") {
        // Group data by block with identical 1) missing variable structure
        // and 2) zero-value structure
        local missvars
        local zerovars
        foreach v of varlist `tmpuniqvars' {
            tempvar miss`v' zero`v'
            qui generate `miss`v'' = missing(`v') if `touse'
            qui generate `zero`v'' = (abs(`v') <= `zero') if `touse'
            local missvars `missvars' `miss`v''
            local zerovars `zerovars' `zero`v''
        }

        tempvar group
        qui gegen `group' = group(`missvars' `zerovars') if `touse'

        qui glevelsof `group' if `touse', local(groups)

        foreach g of local groups {
            tempname gind
            qui generate `gind' = (`g' == `group') & `touse'

            mata: force_identities()

            drop `gind'
        }
        drop `group' `missvars' `zerovars'

        // Show reports
        local colwidth = 12
        local width2   = 4*`colwidth' + 1
        local width3   = round(`colwidth'/2)
        local width4   = `dbcolwidth' - `width3'
        local width5   = `colwidth' - 3

        di ""
        di as text "Summary of adjustments"
        di as text "{hline 61}{c TT}{hline `width2'}"
        di as text "{ralign 61:}{c |}", _continue
        di as text "{ralign `colwidth': mean}{ralign `colwidth': sd}",, _continue
        di as text "{ralign `colwidth': min}{ralign `colwidth': max}"
        di as text "{hline 61}{c +}{hline `width2'}"

        local i = 1
        while ("`idenexpr`i''" != "") {
            local s = substr("`idenexpr`i''", 1, 59)
            di as text "{lalign 60: {bf:`s'}} {c |}"

            di as text "{ralign 60: pre-adjustment absolute discrepancies} {c |} ",, _continue

            tempvar discr
            qui generate `discr' = .
            mata: calc_discrepancy(st_local("uniqvars"))
            qui summarize `discr' if `touse'
            drop `discr'

            di as text "{space 3}" as result %9.3g = r(mean),, _continue
            di as text "{space 3}" as result %9.3g = r(sd),, _continue
            di as text "{space 3}" as result %9.3g = r(min),, _continue
            di as text "{space 3}" as result %9.3g = r(max)

            di as text "{ralign 60: post-adjustment absolute discrepancies} {c |} ",, _continue

            tempvar discr
            qui generate `discr' = .
            mata: calc_discrepancy(st_local("tmpuniqvars"))
            qui summarize `discr' if `touse'
            drop `discr'

            di as text "{space 3}" as result %9.3g = r(mean),, _continue
            di as text "{space 3}" as result %9.3g = r(sd),, _continue
            di as text "{space 3}" as result %9.3g = r(min),, _continue
            di as text "{space 3}" as result %9.3g = r(max)

            local i = `i' + 1
        }

        di as text "{hline 61}{c +}{hline `width2'}"

        foreach v of local uniqvars {
            di as text "{lalign 60: {bf:`v'}} {c |}"

            di as text "{ralign 60: before adjustment} {c |} ",, _continue
            qui summarize `v' if `touse'
            di as text "{space 3}" as result %9.3g = r(mean),, _continue
            di as text "{space 3}" as result %9.3g = r(sd),, _continue
            di as text "{space 3}" as result %9.3g = r(min),, _continue
            di as text "{space 3}" as result %9.3g = r(max)

            di as text "{ralign 60: after adjustment (excl. formerly missing)} {c |} ",, _continue
            qui summarize `tmp`v'' if `touse' & !missing(`v')
            di as text "{space 3}" as result %9.3g = r(mean),, _continue
            di as text "{space 3}" as result %9.3g = r(sd),, _continue
            di as text "{space 3}" as result %9.3g = r(min),, _continue
            di as text "{space 3}" as result %9.3g = r(max)

            di as text "{ralign 60: after adjustment (incl. formerly missing)} {c |} ",, _continue
            qui summarize `tmp`v'' if `touse'
            di as text "{space 3}" as result %9.3g = r(mean),, _continue
            di as text "{space 3}" as result %9.3g = r(sd),, _continue
            di as text "{space 3}" as result %9.3g = r(min),, _continue
            di as text "{space 3}" as result %9.3g = r(max)

            di as text "{ralign 60: difference (after minus before)} {c |} ",, _continue
            tempvar change
            qui generate `change' = `tmp`v'' - `v' if `touse'
            qui summarize `change' if `touse'
            di as text "{space 3}" as result %9.3g = r(mean),, _continue
            di as text "{space 3}" as result %9.3g = r(sd),, _continue
            di as text "{space 3}" as result %9.3g = r(min),, _continue
            di as text "{space 3}" as result %9.3g = r(max)

            drop `change'
        }
        di as text "{hline 61}{c BT}{hline `width2'}"
    }

    // ---------------------------------------------------------------------- //
    // Create output variables
    // ---------------------------------------------------------------------- //

    if ("`suffix'" != "") {
        foreach v of local uniqvars {
            qui generate `v'`suffix' = `tmp`v'', after(`v')
        }
    }
    else if ("`prefix'" != "") {
        foreach v of local uniqvars {
            qui generate `prefix'`v' = `tmp`v'', after(`v')
        }
    }
    else if ("`replace'" != "") {
        foreach v of local uniqvars {
            qui replace `v' = `tmp`v''
        }
    }

end

mata:

void calc_discrepancy(string scalar varlist) {
    i = strtoreal(st_local("i"))

    st_view(vars, ., varlist, st_local("touse"))

    iden = st_matrix(st_local("matiden"))[i, .]

    // Remove variables with a coefficient equal to zero, so that it doesn't
    // create undue missing values
    zerovars = (iden :== 0)

    st_store(., st_local("discr"), st_local("touse"), abs(vars[., selectindex(!zerovars)]*iden[selectindex(!zerovars)]'))
}

void check_identities() {
    tolerance = strtoreal(st_local("tolerance"))

    // Identities in matrix form
    matiden = st_matrix(st_local("matiden"))

    // Variable names
    varnames = st_matrixcolstripe(st_local("matiden"))[., 2]

    // First, look at the system ignoring wether variables are fixed or not

    // Full SVD decomposition to analyse the system and get the nullspace
    fullsvd(matiden, U, s, V)
    _transpose(V)
    rank = rank_from_singular_values(s, tol = tolerance)
    if (rank < cols(matiden)) {
        nullspace = V[., (rank + 1)::cols(matiden)]
    } else {
        nullspace = J(cols(matiden), 0, 0)
    }

    // Store nullspace for fixed variables for later
    fixedvaridx = st_matrix(st_local("fixedvaridx"))

    if (rank == cols(matiden)) {
        // System has full rank
        display("{txt}warning: the identities imply that all the variables are equal to zero")
    } else {
        // System is rank-deficient, so there is some moving parts, but some
        // variables might still be fixed. Look at the nullspace of the
        // matrix to see if this is the case
        tol = 1000*cols(V)*epsilon(1)*tolerance
        for (i = 1; i <= cols(matiden); i++) {
            if (all(abs(nullspace[i, .]) :<= tol)) {
                display("{txt}warning: the identities imply {bf:" + varnames[i, 1] + " = 0}")
            }
            for (j = i + 1; j <= cols(matiden); j++) {
                if (all(abs(nullspace[i, .] :- nullspace[j, .]) :<= tol)) {
                    display("{txt}warning: the identities imply {bf:" + varnames[i, 1] + " = " + varnames[j, 1] + "}")
                }
                if (all(abs(nullspace[i, .] :+ nullspace[j, .]) :<= tol)) {
                    display("{txt}warning: the identities imply {bf:" + varnames[i, 1] + " = -" + varnames[j, 1] + "}")
                }
            }
        }
    }

    // Additional checks needed if there are fixed variables
    if (any(fixedvaridx :== 1)) {
        // Check that the rank of the system is less than the number of variables
        // (so that the system is not fully determined)
        matidennofix = matiden[., selectindex(!fixedvaridx)]
        nofixvarnames = varnames[selectindex(!fixedvaridx), 1]

        fullsvd(matidennofix, U, s, V)
        _transpose(V)
        rank = rank_from_singular_values(s, tol = tolerance)
        // Store the range of the system for later
        nofixrange = U[., 1::rank]

        if (rank == cols(matidennofix)) {
            display("{txt}warning: fixed variables perfectly determine nonfixed variables; their initial values will be irrelevant")
        } else {
            // System is rank-deficient, so there is some moving parts even
            // after fixing fixed variables, but some variables may still
            // be perfectly determined
            nullspace = V[., (rank + 1)::cols(matidennofix)]
            tol = 1000*cols(V)*epsilon(1)*tolerance
            for (i = 1; i <= cols(matidennofix); i++) {
                if (all(abs(nullspace[i, .]) :<= tol)) {
                    display("{txt}warning: the identities imply that {bf:" + nofixvarnames[i, 1] + ///
                        "} is perfectly determined by the fixed variables; its initial value will be irrelevant")
                }
            }
        }

        // Check that the system is solvable for each observations

        // Import variables from the dataset
        st_view(vars, ., st_local("tmpuniqvars"), st_local("touse"))
        // Variables for storing potential issue with observations
        st_view(check_sys, ., st_local("check_sys"), st_local("touse"))
        st_view(check_zer, ., st_local("check_zer"), st_local("touse"))

        // Identities with only fixed variables selected
        matidenfix = matiden[., selectindex(fixedvaridx)]
        rhs = -vars[., selectindex(fixedvaridx)]*matidenfix'

        for (i = 1; i <= rows(vars); i++) {
            // Remove constraints with a missing RHS
            missrhs = (rhs[i, .] :>= .)

            // First: rank of the coefficient matrix of the system
            fullsvd(matidennofix[selectindex(!missrhs), .], U, s, V)
            rkcoef = rank_from_singular_values(s, tol = tolerance)

            // Second: rank of the augmented matrix
            fullsvd((matidennofix[selectindex(!missrhs), .], rhs[i, selectindex(!missrhs)]'), U, s, V)
            rkaug = rank_from_singular_values(s, tol = tolerance)

            // Apply the Rouche–Capelli theorem
            if (rkaug > rkcoef) {
                // Signal error
                check_sys[i, 1] = 1
            } else {
                // Zero-value variables are implicitely fixed, so they can prevent
                // the system from being solvable: we check if this is the case
                zerovars = (abs(vars[i, selectindex(!fixedvaridx)]) :<= strtoreal(st_local("zero")))

                // Check that the RHS of the system belongs to its range even
                // after removing zero-valued variables
                matidennozero = matidennofix[selectindex(!missrhs), selectindex(!zerovars)]

                fullsvd(matidennozero, U, s, V)
                rknozero = rank_from_singular_values(s, tol = tolerance)
                fullsvd((matidennozero, rhs[i, selectindex(!missrhs)]'), U, s, V)
                rknozeroaug = rank_from_singular_values(s, tol = tolerance)

                // Use the Rouche–Capelli theorem
                if (rknozeroaug > rknozero) {
                    // Signal error
                    check_zer[i, 1] = 1
                }
            }
        }
    }
}

void fill_missing() {
    tolerance = strtoreal(st_local("tolerance"))

    // Identities in matrix form
    matiden = st_matrix(st_local("matiden"))

    // View to all the variables of the problem
    st_view(vars, ., st_local("tmpuniqvars"), st_local("gind"))

    // View to the variable with the missing strcuture of the data
    st_view(missvars, ., st_local("missvars"), st_local("gind"))
    missstruct = missvars[1, .] // They are all the same, just need the first line

    // Ignore trivial cases
    if (all(missstruct :== 0)) {
        // Nothing missing, nothing to be done
        return
    } else if (all(missstruct :== 1)) {
        // Everything missing, nothing can be done
        return
    }

    // Create the proper system of equations
    A = matiden[., selectindex(missstruct)]
    B = -matiden[., selectindex(!missstruct)]*vars[., selectindex(!missstruct)]'

    // Find a generalized solution of the system: the use of generalized
    // solution allows us to get a reasonable solution even if the system
    // isn't yet consistent
    X = svsolve(A, B, rank, tol = tolerance)

    // Look at the nullspace of the matrix to see which variables may be
    // fully determined
    fullsvd(A, U, s, V)
    _transpose(V)
    rank = rank_from_singular_values(s, tol = tolerance)
    missvaridx = selectindex(missstruct)

    if (rank == cols(A)) {
        // System is perfectly determined, we fill all variables
        vars[., missvaridx] = X'
    } else {
        // System is undetermined, but some variables might be determined,
        // and for that we must check the nullspace
        nullspace = V[., (rank + 1)::cols(V)]
        tol = 1000*cols(V)*epsilon(1)*tolerance
        for (i = 1; i <= cols(A); i++) {
            // For fully determined variables, all coefficients of the basis of
            // the nullspace are zero
            if (all(abs(nullspace[i, .]) :<= tol)) {
                vars[., missvaridx[i]] = X[i, .]'
            }
        }
    }
}

void force_identities() {
    tolerance = strtoreal(st_local("tolerance"))

    // Identities in matrix form
    matiden = st_matrix(st_local("matiden"))

    // Variable names
    varnames = st_matrixcolstripe(st_local("matiden"))[., 2]

    // Fixed variables of the system
    fixedvaridx = st_matrix(st_local("fixedvaridx"))

    // View to all the variables of the problem
    st_view(vars, ., st_local("tmpuniqvars"), st_local("gind"))

    // View to the variable with the missing strcuture of the data
    st_view(missvars, ., st_local("missvars"), st_local("gind"))
    missstruct = missvars[1, .] // They are all the same, just need the first line

    // View to the variable with the zero strcuture of the data
    st_view(zerovars, ., st_local("zerovars"), st_local("gind"))
    zerostruct = zerovars[1, .] // They are all the same, just need the first line

    // Construct a system of equality that only involves nonmissing variables
    fullsvd(matiden, U, s, V)
    _transpose(V)
    rank = rank_from_singular_values(s, tol = tolerance)
    nullspace = V[selectindex(!missstruct), (rank + 1)::cols(V)]
    fullsvd(nullspace, U, s, V)
    rank = rank_from_singular_values(s, tol = tolerance)
    if (rank == cols(U)) {
        // No enforcable constraints (all variables linearly independent)
        return
    }
    // Get group-specific constraints (left null space = orthogonal complement of the range)
    matidengrp = U[., (rank + 1)::cols(U)]'

    // Fixed and moving within all vars
    varfix = (fixedvaridx :| zerostruct) :& !missstruct
    varnfix = (!fixedvaridx :& !zerostruct) :& !missstruct

    // Fixed and moving variables within nonmissings
    grpfix = fixedvaridx[1, selectindex(!missstruct)] :| zerostruct[1, selectindex(!missstruct)]
    grpnfix = !fixedvaridx[1, selectindex(!missstruct)] :& !zerostruct[1, selectindex(!missstruct)]

    // Estimate the LHS and the RHS of the equality constraints
    lhs = matidengrp[., selectindex(grpnfix)]
    rhs = -matidengrp[., selectindex(grpfix)]*vars[, selectindex(varfix)]'

    // Enforce the constraints
    for (i = 1; i <= rows(vars); i++) {
        // Matrix of quadratic coefficients
        quadcoefs = 1 :/ abs(vars[i, selectindex(varnfix)])
        // Ensure reasonable conditioning
        meancoef = mean(quadcoefs')
        Q = diag(quadcoefs/meancoef)

        // Vector of linear coefficients
        c = sign(vars[i, selectindex(varnfix)])'/meancoef

        // Build the system to be solved
        A = (Q, lhs' \ lhs, J(rows(lhs), rows(lhs), 0))
        b = (c \ rhs[, i])

        // Solve
        rank = _qrsolve(A, b, tol = tolerance)

        vars[i, selectindex(varnfix)] = b[1::cols(Q), 1]'
    }
}

end