*! version 1.1.1 04June2017 Mauricio Caceres Bravo, caceres@nber.org
*! Mata files for alogit.ado

version 13.1

///////////////////////////////////////////////////////////////////////
//                       Main alogit function                        //
///////////////////////////////////////////////////////////////////////

// capture mata: mata drop alogit()
// capture mata: mata drop _alogit_exact()
// capture mata: mata drop _alogit_exact_fast()
// capture mata: mata drop _alogit_exact_loop()
//
// capture mata: mata drop alogit_logisticden()
// capture mata: mata drop alogit_combos()
// capture mata: mata drop _alogit_binadd_reverse()

// dlabels = `dlabels'
// aloglik = &`mata_loglik'()
// predict = &`mata_predict'()

mata

// Main alogit function to perform the optimization
transmorphic function alogit(string rowvector dlabels,
                             pointer(function) aloglik,
                             pointer(function) predict)
{
    transmorphic M

    // Stata options
    string scalar model, method, algorithm, defgood, consider
    real scalar asexp, getc, getdef, reps, pll, fast

    // Data from Stata
    real matrix x, xz, info
    real colvector y, defcons
    real rowvector b

    // Auxiliary variables for optimization
    real colvector phi0, Pc0, nj_i, nc_i
    real matrix C, U0, C0, PC0

    // moptimize options
    real scalar difficult,
                technique,
                iterate,
                tolerance,
                ltolerance,
                nrtolerance

    string scalar exp,
                  nolog,
                  gradient,
                  hessian,
                  showstep,
                  showtolerance,
                  coefdiffs,
                  trace,
                  nonrtolerance,
                  norescale

    // Auxiliary variables for this function
    real scalar loglik, i, nj
    string scalar nstr, jstr
    real rowvector panel

    // Parse options from Stata
    // ------------------------

    exp       = st_local("exp")             // Exponential version
    model     = st_local("model")           // Model (dsc or alogit)
    method    = st_local("method")          // Method (exact, importance)
    algorithm = st_local("algorithm")       // Algorithm (matrix, loop, fast, faster)
    reps      = strtoreal(st_local("reps")) // Repetitions for importance
    defgood   = st_local("default")         // Default good
    consider  = st_local("consider")        // Consider goods

    // whether to compute the likelihood for exp(b)
    asexp = ( exp == "exp" )

    // Whether to process a default good or consideration goods
    // getcons = (consider != "") & (defgood == "")
    getdef = (defgood != "") & (consider == "")

    // Print the likelihood; don't optimize
    pll  = (st_local("loglik") == "loglik")

    // Will only generate C for model(alogit) and alg(matrix)
    getc = ((model == "alogit") & (algorithm == "matrix"))

    // If the computations are done in C, skip (*predict) since that
    // will be called directly from alogit.ado
    fast = ((algorithm == "fast") | (algorithm == "faster"))

    // Initialize optimization programs
    // --------------------------------

    M = moptimize_init()

    // Turn off intercepts (will add constant in .ado files as applicable)
    moptimize_init_eq_cons(M, 1, "off")
    moptimize_init_eq_cons(M, 2, "off")

    // Mark sample to use; grab grouping variable
    moptimize_init_touse(M, st_local("touse"))
    moptimize_init_by(M, st_local("group"))

    // Set up choice sets by grouping variable
    info  = panelsetup(*moptimize_util_by(M), 1)
    panel = panelstats(info)
    if (panel[3] == panel[4]) {
        nstr = strtrim(sprintf("%21.0gc", panel[1]))
        jstr = strtrim(sprintf("%21.0gc", panel[3]))
        printf("Balanced panel. n = " + nstr + ", j = " + jstr + "\n")
    }
    else {
        nstr = strtrim(sprintf("%21.0gc", panel[1]))
        jstr = strtrim(sprintf("%21.0gc", panel[3]))
        jstr = jstr + " to " + strtrim(sprintf("%21.0gc", panel[4]))
        printf("Unbalanced panel. n = " + nstr + ", j = " + jstr + "\n")
    }

    // Set up optimization variables
    // -----------------------------

    if (fast) {
        moptimize_init_userinfo(M, 1, algorithm)
        moptimize_init_userinfo(M, 8, asexp)
    }
    else {
        // Grab the data
        x  = st_data(., st_local("xvars"),  st_local("touse"))
        xz = st_data(., st_local("xzvars"), st_local("touse"))
        // xb  = x  * st_matrix("b0")'
        // xzd = xz * st_matrix("d0")'

        // Default good
        defcons = st_data(., defgood + consider,  st_local("touse"))
        moptimize_init_userinfo(M, 9, defcons)

        if (method == "exact") {
            // Aux for computing choice sets
            nj   = panel[4]
            nj_i = info[, 2] :- info[, 1] :+ 1
            nc_i = getc? (2:^nj_i :- 1): .
            C    = getc? alogit_combos(nj): .

            // Pass user variables to compute the exact likelihood
            moptimize_init_userinfo(M, 1, C)
            moptimize_init_userinfo(M, 2, nc_i)
            moptimize_init_userinfo(M, 3, nj_i)
            moptimize_init_userinfo(M, 4, getdef)
            moptimize_init_userinfo(M, 5, x)
            moptimize_init_userinfo(M, 6, xz)
            moptimize_init_userinfo(M, 7, info)
            moptimize_init_userinfo(M, 8, asexp)
        }
        else if (method == "importance") {
            // Get a random sample of choice sets
            phi0 = invlogit(xz * st_matrix("d0")')
            U0   = runiform(rows(xz), reps)
            C0   = phi0 :> U0
            Pc0  = C0 :* log(phi0) :+ (1 :- C0) :* log(1 :- phi0)

            // Get the probabilities of the choice sets sample
            PC0  = J(rows(info), reps, .)
            for (i = 1; i <= rows(info); i++) {
                PC0[i, ] = exp(colsum(panelsubmatrix(Pc0, i, info)))
            }

            // Pass user variables to compute the importance likelihood
            moptimize_init_userinfo(M, 1, C0)
            moptimize_init_userinfo(M, 2, PC0)
            moptimize_init_userinfo(M, 3, reps)
            moptimize_init_userinfo(M, 4, U0)
            moptimize_init_userinfo(M, 5, x)
            moptimize_init_userinfo(M, 6, xz)
            moptimize_init_userinfo(M, 7, info)
            moptimize_init_userinfo(M, 8, 0)
        }
    }

    // Set up optimization
    // -------------------

    // Set the evaluator to compute the importance likelihood
    moptimize_init_evaluator(M, aloglik)
    moptimize_init_evaluatortype(M, st_local("eval"))

    // Set dependent variable, independent vars
    moptimize_init_depvar(M, 1, st_local("depvar"))
    moptimize_init_eq_indepvars(M, 1, st_local("xvars"))
    moptimize_init_eq_indepvars(M, 2, st_local("xzvars"))

    // Label the things
    moptimize_init_eq_name(M, 1, "PY")
    moptimize_init_eq_name(M, 2, "PA")
    moptimize_init_eq_colnames(M, 2, dlabels)

    // Set up the starting parameters
    moptimize_init_eq_coefs(M, 1, asexp? exp(st_matrix("b0")): st_matrix("b0"))
    moptimize_init_eq_coefs(M, 2, asexp? exp(st_matrix("d0")): st_matrix("d0"))

    // Make the display pretty
    moptimize_init_valueid(M, "log likelihood")

    // Pass user options to moptimize()
    // --------------------------------

    difficult   = st_numscalar("difficult")
    technique   = st_local("technique")
    iterate     = st_numscalar("iterate")
    tolerance   = strtoreal(st_local("tolerance"))
    ltolerance  = strtoreal(st_local("ltolerance"))
    nrtolerance = strtoreal(st_local("nrtolerance"))

    if (difficult)        moptimize_init_singularHmethod(M, "hybrid")
    if (technique != "")  moptimize_init_technique(M, technique)
    if (iterate >= 0)     moptimize_init_conv_maxiter(M, iterate)
    if (tolerance   > 0)  moptimize_init_conv_ptol(M , tolerance)
    if (ltolerance  > 0)  moptimize_init_conv_vtol(M , ltolerance)
    if (nrtolerance > 0)  moptimize_init_conv_nrtol(M, nrtolerance)

    nolog         = (st_local("nolog")         == "")? "on":  "off"
    gradient      = (st_local("gradient")      == "")? "off": "on"
    hessian       = (st_local("hessian")       == "")? "off": "on"
    showstep      = (st_local("showstep")      == "")? "off": "on"
    showtolerance = (st_local("showtolerance") == "")? "off": "on"
    coefdiffs     = (st_local("coefdiffs")     == "")? "off": "on"
    trace         = (st_local("trace")         == "")? "off": "on"
    nonrtolerance = (st_local("nonrtolerance") == "")? "off": "on"
    norescale     = (st_local("norescale")     == "")? "on":  "off"

    moptimize_init_trace_value(M,      nolog)
    moptimize_init_trace_gradient(M,   gradient)
    moptimize_init_trace_Hessian(M,    hessian)
    moptimize_init_trace_step(M,       showstep)
    moptimize_init_trace_tol(M,        showtolerance)
    moptimize_init_trace_coefdiffs(M,  coefdiffs)
    moptimize_init_trace_coefs(M,      trace)
    moptimize_init_conv_ignorenrtol(M, nonrtolerance)
    moptimize_init_search_rescale(M ,  norescale)
    moptimize_init_nmsimplexdeltas(M, 1)

    // Run the routine and put the results in Stata
    // --------------------------------------------

    if (pll) {
        // Print likelihood at default/user-provided parameters
        b = asexp? exp((st_matrix("b0"), st_matrix("d0"))): st_matrix("b0"), st_matrix("d0")
        y = moptimize_init_depvar(M, 1)
        loglik = (*predict)(M, y, x, xz, defcons, method, "ll", b)
        printf("log-likelihood = " + strtrim(sprintf("%21.4f\n", loglik)))
        st_numscalar("ll", loglik)
    }
    else {
        // Run the optimization
        printf("\n")
        moptimize(M)
        if (!fast) {

            // Get average P(empty set)
            y = moptimize_init_depvar(M, 1)
            st_numscalar("mp0", (*predict)(M, y, x, xz, defcons, method, "mp0"))

            // Get average P(attention)
            if (model == "dsc") {
                st_numscalar("mp1", (*predict)(M, y, x, xz, defcons, method, "mp1"))
            }
            else {
                st_numscalar("mp1", .)
            }
        }
        // For fast, these are done in alogit.ado
    }
    return(M)
}

///////////////////////////////////////////////////////////////////////
//                         Exact Likelihood                          //
///////////////////////////////////////////////////////////////////////

// Loglik for exact alogit, matrix version
void function _alogit_exact(transmorphic M,
                            real scalar todo,
                            real rowvector b,
                            fv, g, H)
{

    // Function variables
    real matrix x, xz, info, C, xbmat, xzdmat
    real colvector nj_i, nc_i
    real colvector y, defcons, xb, xzd
    real rowvector bb, gd
    real scalar asexp, kx, kxz, ktot, getdef, nc_0

    // Function helpers
    real matrix xdphi_ij, xphi_ij
    real colvector exb, phi_ij, lphi_c, lphi_nc

    // Loop variables to use for the likelihood
    real scalar i
    real matrix Ci, Ciall, PY
    real colvector y_ij, exb_ij, dc_ij
    real colvector Cisel, lPA
    real rowvector PA, PY_ij

    // Loop variables to use for the gradient
    real colvector dLdP
    real rowvector dLdP_aux, dLdb_k, dLdg_k
    real matrix x_ij, xz_ij, xphi_i
    real matrix dPdb_k_aux, dPdg_k_aux

    // Loop variables to use for the Hessian
    real scalar l, k
    real matrix hb, hg, hbg, haux
    real matrix dPdb_l0, dLdb_kl1, dLdb_kl2, dLdb_kl3, dLdb_kl0, dPdb_kl
    real matrix dPdb_kl_aux, dPdg_lk_aux, dPdbg_lk
    real matrix dPdg_kl0, dLdg_kl

    // Variables to use in the computation
    // -----------------------------------

    C       = moptimize_init_userinfo(M, 1) // Choice sets
    nc_i    = moptimize_init_userinfo(M, 2) // Goods per group
    nj_i    = moptimize_init_userinfo(M, 3) // 2^nc_i - 1
    getdef  = moptimize_init_userinfo(M, 4) // Get default or consider
    x       = moptimize_init_userinfo(M, 5) // utility vars
    xz      = moptimize_init_userinfo(M, 6) // attention vars
    info    = moptimize_init_userinfo(M, 7) // pos of each group in data
    asexp   = moptimize_init_userinfo(M, 8) // Estimate exp(b) isntead of b
    defcons = moptimize_init_userinfo(M, 9) // default good/consider goods
    y       = moptimize_util_depvar(M, 1)   // outcome (choice)
    kx      = cols(x)
    kxz     = cols(xz)
    ktot    = kx + kxz

    if ( asexp ) {
        bb      = b[|1 \ kx|]       // beta
        gd      = b[|kx + 1 \ .|]   // gamma, delta
        xbmat   = (bb :^ x)         // exp(b)^x * b
        xzdmat  = (gd :^ xz)        // exp(g, d)^(x, z)
        xb      = xbmat[, 1]        // prod(exp(b)^x over kx) = exp(xb)
        xzd     = xzdmat[, 1]       // prod(exp(g, d)^(x, z) over kx + kz) = exp(xg + zd)
        for (k = 2; k <= kx; k++)  xb  = xb  :* xbmat[, k]
        for (k = 2; k <= kxz; k++) xzd = xzd :* xzdmat[, k]
    }
    else {
        xb  = moptimize_util_xb(M, b, 1)  // x * b
        xzd = moptimize_util_xb(M, b, 2)  // x * g + z * d
    }

    phi_ij = asexp? (xzd :/ (1 :+ xzd)): invlogit(xzd) // Get phi for P(attention)
    if (!getdef) phi_ij[selectindex(defcons), ] = J(sum(defcons), 1, 1)

    _editvalue(phi_ij, 1, 1 - epsilon(1)) // Numerical precision foo
    _editvalue(phi_ij, 0, epsilon(1))     // Numerical precision foo

    xphi_ij  = phi_ij :* xz       // Helper for the gradient
    exb      = asexp? xb: exp(xb) // Helper for P(Y)
    lphi_c   = log(phi_ij)        // Helpers for P(A) since Mata
    lphi_nc  = log(1 :- phi_ij)   // does not have a prod function.

    // Helper for the Hessian
    xdphi_ij = ( asexp? (xzd :/ (1 :+ xzd):^2): alogit_logisticden(xzd) ) :* xz
    if (!getdef) xdphi_ij[selectindex(defcons), ] = J(sum(defcons), kxz, 0)

    H  = J(ktot, ktot,  0) // Hessian
    g  = J(1, ktot,  0)    // Gradient
    fv = 0                 // Log likelihood

    // Loop over individuals
    // ---------------------

    for (i = 1; i <= rows(info); i++) {
        // Drop individuals with > 1 or 0 goods chosen
        y_ij = panelsubmatrix(y, i, info)
        if(sum(y_ij) != 1) continue
        dc_ij = panelsubmatrix(defcons, i, info)

        // Get choice sets by individual
        nc_0  = getdef? 1: nc_i[i] - 2^(nj_i[i] - sum(dc_ij)) + 2
        Ciall = C[|nc_0, 1 \ nc_i[i] + 1, nj_i[i]|]
        Cisel = selectindex(Ciall[, selectindex(y_ij)])
        Ci    = J(1, nj_i[i], 0) \ Ciall[Cisel, ]

        // Compute the likelihood
        // ----------------------

        // Probability of each choice set. We want to compute P(Y |
        // c, theta) P(c | theta); we compute P(c | theta) / (total
        // for P(Y)) and then sum over all the choice sets
        lPA    = Ci * panelsubmatrix(lphi_c, i, info) +
                 (1 :- Ci) * panelsubmatrix(lphi_nc, i, info)
        PA     = rowshape(exp(lPA), 1)
        exb_ij = panelsubmatrix(exb, i, info)
        PY     = rowshape(exb_ij, 1) :* (Ci :/ (Ci * exb_ij))

        _editmissing(PY, 0)
        if (getdef) PY[1, selectindex(dc_ij)] = 1
        PY_ij = PA * PY

        // Compute the likelihood
        fv = fv + log(PY_ij) * y_ij

        // Compute the gradient and the Hessian
        // ------------------------------------

        if (todo >= 1) {
            dLdP     = y_ij :/ colshape(PY_ij, 1)
            dLdP_aux = PY * dLdP
            x_ij     = panelsubmatrix(x,  i, info)
            xz_ij    = panelsubmatrix(xz, i, info)
            xphi_i   = panelsubmatrix(xphi_ij, i, info)'

            // Compute the gradient
            dPdb_k_aux = (PY * x_ij)'
            dPdg_k_aux = ((Ci * xz_ij) :- rowshape(rowsum(xphi_i), 1))'

            dLdb_k = PY_ij * (x_ij  :* dLdP) - ((dPdb_k_aux :* PA) * dLdP_aux)'
            dLdg_k = ((PA :* dPdg_k_aux) * dLdP_aux)'
            g = g + (dLdb_k, dLdg_k)

            // Compute the Hessian
            if (todo == 2) {
                hb  = dLdb_k' * dLdb_k
                hg  = dLdg_k' * dLdg_k
                hbg = dLdb_k' * dLdg_k

                // For betas
                for (k = 1; k <= kx; k++) {
                    for (l = k; l <= kx; l++) {
                        dPdb_l0  = PY :* x_ij[, l]' - PY :* (PY * x_ij[, l])
                        dLdb_kl1 = dPdb_l0 :* x_ij[, k]'
                        dLdb_kl2 = dPdb_l0 :* colshape(dPdb_k_aux[k, ], 1)
                        dLdb_kl3 = PY :* (dPdb_l0 * x_ij[, k])
                        dLdb_kl0 = dLdb_kl1 -  dLdb_kl2 - dLdb_kl3
                        dPdb_kl  = (PA * dLdb_kl0) * dLdP
                        H[k, l]  = H[k, l] + dPdb_kl - hb[k, l]
                    }

                    // Cross partials
                    dPdb_kl_aux = PY :* x_ij[, k]' - (PY * x_ij[, k]) :* PY
                    for (l = 1; l <= kxz; l++) {
                        dPdg_lk_aux = rowshape((Ci * xz_ij[, l]) :- sum(xphi_i[l, ]), 1) :* PA
                        dPdbg_lk = (dPdg_lk_aux * dPdb_kl_aux) * dLdP
                        H[k, kx + l] = H[k, kx + l] + dPdbg_lk - hbg[k, l]
                    }
                }

                // For gamma, delta
                haux = uppertriangle(panelsubmatrix(xdphi_ij, i, info)' * xz_ij)
                for (k = 1; k <= kxz; k++) {
                    for (l = k; l <= kxz; l++) {
                        dPdg_kl0   = (dPdg_k_aux[l, ] :* dPdg_k_aux[k, ] :- haux[k, l]) :* PA
                        dLdg_kl    = dPdg_kl0 * dLdP_aux - hg[k, l]
                        H[kx + k, kx + l] = H[kx + k, kx + l] + dLdg_kl
                    }
                }
            }
        }
    }
    if ( todo == 2 ) H = makesymmetric(uppertriangle(H)')
    if ( asexp ) {
        if (todo >= 1) g = g :/ b
        if (todo == 2) H = H :/ (b' * b) :- diag(g :/ b)
    }
}

// Compute the likelihood, gradient, and Hessian in C
void function _alogit_exact_fast(transmorphic M,
                                 real scalar todo,
                                 real rowvector b,
                                 fv, g, H)
{
    real scalar ktot
    ktot = st_numscalar("__alogit_ktot")
    H    = J(ktot, ktot,  0);
    g    = J(1, ktot,  0);
    st_numscalar("__alogit_ll", .)
    st_matrix("__alogit_b", b)
    st_numscalar("__alogit_todo", todo)
    if (moptimize_init_userinfo(M, 1) == "fast") {
        stata("plugin call aloglik_fast " + st_local("fastcall") + " alogit")
    }
    else {
        stata("plugin call aloglik_faster " + st_local("fastcall") + " alogit")
    }
    if (todo >= 1) g = st_matrix("__alogit_g")
    if (todo == 2) H = makesymmetric(uppertriangle(st_matrix("__alogit_H"))')
    fv = st_numscalar("__alogit_ll")
    if ( moptimize_init_userinfo(M, 8) ) {
        if (todo >= 1) g = g :/ b
        if (todo == 2) H = H :/ (b' * b) :- diag(g :/ b)
    }
}

// Loglik for exact alogit, loop version
void function _alogit_exact_loop(transmorphic M,
                                 real scalar todo,
                                 real rowvector b,
                                 fv, g, H)
{

    // Function variables
    real matrix x, xz, info, xbmat, xzdmat
    real colvector y, defcons, xb, xzd, nj_i
    real rowvector bb, gd
    real scalar asexp, kx, kxz, ktot, getdef

    // Function helpers
    real matrix xdphi_ij, xphi_ij
    real colvector exb, phi_ij, lphi_c, lphi_nc

    // Loop variables to use for the likelihood
    real scalar i, k, sel
    real scalar PY_ijc, PY_ij, PA
    real colvector y_ij, exb_ij
    real colvector lphi_c_ij, lphi_nc_ij, PY
    real rowvector c

    // Loop variables to use for the gradient
    real matrix x_ij, xz_ij, xphi_i, PYx_ij_mat
    real rowvector x_ij_sel, PYx_ij, cxz_ij, xphi_isum
    real rowvector dPdb_aux, dPdg_aux
    real rowvector dLdb_k, dLdg_k

    // Loop variables to use for the Hessian
    real matrix hb, hg, hbg, haux, xdphi_i
    real rowvector x_ijPY
    real matrix dPdb_kl_1, dPdb_kl_2, dLdb_aux
    real matrix dPdg_l_aux, dLdbg_aux
    real matrix dLdg_aux

    // Variables to use in the computation
    // -----------------------------------

    // Get choice sets, xb, and xg + zd
    nj_i    = moptimize_init_userinfo(M, 3) // Goods per group
    getdef  = moptimize_init_userinfo(M, 4) // Get default or consider
    x       = moptimize_init_userinfo(M, 5) // utility vars
    xz      = moptimize_init_userinfo(M, 6) // attention vars
    info    = moptimize_init_userinfo(M, 7) // pos of each group in data
    asexp   = moptimize_init_userinfo(M, 8) // Estimate exp(b) isntead of b
    defcons = moptimize_init_userinfo(M, 9) // default good/consider goods
    y       = moptimize_util_depvar(M, 1)   // outcome (choice)
    kx      = cols(x)
    kxz     = cols(xz)
    ktot    = kx + kxz

    if ( asexp ) {
        bb      = b[|1 \ kx|]       // beta
        gd      = b[|kx + 1 \ .|]   // gamma, delta
        xbmat   = (bb :^ x)         // exp(b)^x * b
        xzdmat  = (gd :^ xz)        // exp(g, d)^(x, z)
        xb      = xbmat[, 1]        // prod(exp(b)^x over kx) = exp(xb)
        xzd     = xzdmat[, 1]       // prod(exp(g, d)^(x, z) over kx + kz) = exp(xg + zd)
        for (k = 2; k <= kx; k++)  xb  = xb  :* xbmat[, k]
        for (k = 2; k <= kxz; k++) xzd = xzd :* xzdmat[, k]
    }
    else {
        xb  = moptimize_util_xb(M, b, 1)  // x * b
        xzd = moptimize_util_xb(M, b, 2)  // x * g + z * d
    }

    phi_ij = asexp? (xzd :/ (1 :+ xzd)): invlogit(xzd) // Get phi for P(attention)
    if (!getdef) phi_ij[selectindex(defcons), ] = J(sum(defcons), 1, 1)

    _editvalue(phi_ij, 1, 1 - epsilon(1))  // Numerical precision foo
    _editvalue(phi_ij, 0, epsilon(1))      // Numerical precision foo

    xphi_ij  = phi_ij :* xz        // Helper for the gradient
    exb      = asexp? xb: exp(xb)  // Helper for P(Y)
    lphi_c   = log(phi_ij)         // Helper for P(A)
    lphi_nc  = log(1 :- phi_ij)    // Helper for P(A)

    // Helper for the Hessian
    xdphi_ij = ( asexp? (xzd :/ (1 :+ xzd):^2): alogit_logisticden(xzd) ) :* xz
    if (!getdef) xdphi_ij[selectindex(defcons), ] = J(sum(defcons), kxz, 0)

    // Initialize values
    H  = J(ktot, ktot,  0) // Hessian
    g  = J(1, ktot,  0)    // Gradient
    fv = 0                 // Log likelihood

    // Loop over individuals
    // ---------------------

    for (i = 1; i <= rows(info); i++) {

        // Drop individuals with > 1 or 0 goods chosen
        // -------------------------------------------

        y_ij = panelsubmatrix(y, i, info)
        if(sum(y_ij) != 1) continue

        lphi_c_ij  = panelsubmatrix(lphi_c, i, info)
        lphi_nc_ij = panelsubmatrix(lphi_nc, i, info)
        exb_ij     = panelsubmatrix(exb, i, info)
        x_ij       = (todo >= 1)? panelsubmatrix(x,  i, info): .
        xz_ij      = (todo >= 1)? panelsubmatrix(xz, i, info): .
        xphi_i     = (todo >= 1)? panelsubmatrix(xphi_ij, i, info)': .
        xdphi_i    = (todo >= 1)? panelsubmatrix(xdphi_ij, i, info): .

        // Initialize P(Y)
        // ---------------

        sel    = selectindex(y_ij)
        if (getdef) {
            // If you specify a default good, start with the empty set
            c  = J(1, nj_i[i], 0)
            PY = panelsubmatrix(defcons, i, info)
        }
        else {
            // Otherwise, start with the first set that includes the
            // choice and the goods that must be considered.
            c  = panelsubmatrix(defcons, i, info)'
            while (c[sel] == 0) _alogit_binadd_reverse(c)
            PY = c' :* exb_ij :/ (c * exb_ij)
        }
        PA     = exp(c * lphi_c_ij + (1 :- c) * lphi_nc_ij)
        PY_ijc = PA * PY[sel]
        PY_ij  = PY_ijc

        //  Gradient and Hessian
        // ---------------------

        if (todo >= 1) {
            x_ij_sel = x_ij[sel, ]

            // These are defined this way to use in the Hessian
            PYx_ij_mat = PY :* x_ij
            PYx_ij     = colsum(PYx_ij_mat)
            cxz_ij     = c * xz_ij

            // Gradient helpers
            dPdb_aux = PYx_ij * PY_ijc
            dPdg_aux = cxz_ij * PY_ijc

            // Hessian helpers
            xphi_isum = rowshape(rowsum(xphi_i), 1)
            if (todo >= 2) {
                // For beta
                x_ijPY     = x_ij_sel :- PYx_ij
                dPdb_kl_1  = x_ijPY' * x_ijPY
                dPdb_kl_2  = x_ij' * PYx_ij_mat - PYx_ij' * PYx_ij
                dLdb_aux   = (dPdb_kl_1 - dPdb_kl_2) :* PY_ijc

                // Cross partials
                dPdg_l_aux = cxz_ij :- xphi_isum
                dLdbg_aux  = (x_ijPY' * dPdg_l_aux) :* PY_ijc

                // For gamma, delta
                haux       = uppertriangle(xdphi_i' * xz_ij)
                dLdg_aux   = (dPdg_l_aux' * dPdg_l_aux) :* PY_ijc
            }
        }

        // Loop through all consideration sets
        // -----------------------------------

        while (!min(c)) {
             _alogit_binadd_reverse(c)
            if (c[sel] == 0) continue

            // Contribute to P(Y)
            // ------------------

            PY     = c' :* exb_ij :/ (c * exb_ij)
            PA     = exp(c * lphi_c_ij + (1 :- c) * lphi_nc_ij)
            PY_ijc = PA * PY[sel]
            PY_ij  = PY_ij + PY_ijc

            // Gradient and Hessian
            // --------------------

            if (todo >= 1) {

                // Defined this way to use in the Hessian
                PYx_ij_mat = PY :* x_ij
                PYx_ij     = colsum(PYx_ij_mat)
                cxz_ij     = c * xz_ij

                // Gradient helpers
                dPdb_aux = dPdb_aux + PYx_ij * PY_ijc
                dPdg_aux = dPdg_aux + cxz_ij * PY_ijc

                // Hessian
                if (todo == 2) {
                    // For beta
                    x_ijPY     = x_ij_sel :- PYx_ij
                    dPdb_kl_1  = x_ijPY' * x_ijPY
                    dPdb_kl_2  = x_ij' * PYx_ij_mat - PYx_ij' * PYx_ij
                    dLdb_aux   = dLdb_aux + (dPdb_kl_1 - dPdb_kl_2) :* PY_ijc

                    // Cross partials
                    dPdg_l_aux = cxz_ij :- xphi_isum
                    dLdbg_aux  = dLdbg_aux + (x_ijPY' * dPdg_l_aux) :* PY_ijc

                    // For gamma
                    dLdg_aux   = dLdg_aux + (dPdg_l_aux' * dPdg_l_aux) :* PY_ijc
                }
            }
        }

        // Compute the likelihood
        // ----------------------

        fv = fv + log(PY_ij)

        //  Gradient and Hessian
        // ---------------------

        if (todo >= 1) {
            dLdb_k = x_ij_sel :- dPdb_aux :/ PY_ij
            dLdg_k = dPdg_aux :/ PY_ij :- xphi_isum

            g = g + (dLdb_k, dLdg_k)

            // Hessian
            if (todo == 2) {
                hb  = (dLdb_aux  / PY_ij) :- dLdb_k' * dLdb_k
                hbg = (dLdbg_aux / PY_ij) :- dLdb_k' * dLdg_k
                hg  = (dLdg_aux  / PY_ij) :- dLdg_k' * dLdg_k :- haux
                H   = H :+ ((hb, hbg) \ (hbg', hg))
            }
        }

    }
    if (todo == 2) H = makesymmetric(uppertriangle(H)')
    if ( asexp ) {
        if (todo >= 1) g = g :/ b
        if (todo == 2) H = H :/ (b' * b) :- diag(g :/ b)
    }
}

///////////////////////////////////////////////////////////////////////
//                         Helper Functions                          //
///////////////////////////////////////////////////////////////////////

// Logistic density (prior to Stata 14, this was not provided by Mata)
real function alogit_logisticden(real x)
{
    return(invlogit(x) :/ (1 :+ exp(x)))
}

// All the choice sets of nj products in one matrix
real matrix function alogit_combos(real scalar nj)
{
    real scalar nc, i, j
    real matrix C
    string scalar z
    nc = 2^nj
    C  = J(nc, nj, 0)
    for (j = 1; j < nc; j++) {
        z = strreverse(inbase(2, j))
        for (i = 1; i <= strlen(z); i++) {
            C[j + 1, i] = strtoreal(substr(z, i, 1))
        }
    }
    return(C)
}

// Add one to a vector representing a binary number in reverse
// (i.e. first entry is 0s, last is 2^J-1). In testing, using this
// function was faster than using a function that adds one to a vector
// representing a binary number in order.
void function _alogit_binadd_reverse(real rowvector c)
{
    real scalar i
    i = 1
    while (c[i]) {
        c[i] = 0
        i++
    }
    c[i] = 1
}
end
///////////////////////////////////////////////////////////////////////
//                        Alternative Models                         //
///////////////////////////////////////////////////////////////////////

// capture mata: mata drop _alogit_importance()

mata

// Loglik for importance sampler, only one loop
void function _alogit_importance(transmorphic M,
                                 real scalar todo,
                                 real rowvector b,
                                 fv, g, H)
{
    // Function variables
    real matrix x, xz, info, C, PC0
    real colvector y, def, xb, xzd
    real scalar R, kx, kxz, ktot

    // Function helpers
    real matrix xdphi_ij, xphi_ij
    real colvector exb, phi_ij, lphi_c, lphi_nc

    // Loop variables to use for the likelihood
    real scalar i
    string scalar err_noprob
    real matrix Ci, PY
    real colvector y_ij, exb_ij, d_ij
    real colvector lPA
    real rowvector PAPC0, PY_ij, lPY_ij, ll

    // Loop variables to use for the gradient
    real colvector dLdP
    real rowvector dPdb_k, dLdb_k, dPdg_k
    real matrix x_ij, xz_ij, xphi_i
    real matrix dPdb_k_aux, dPdg_k_aux

    // Loop variables to use for the Hessian
    real scalar l, k
    real matrix haux
    real matrix dPdb_l0, dLdb_kl1, dLdb_kl2, dLdb_kl3, dLdb_kl0, dPdb_kl, dLdbg_kl
    real matrix dPdb_kl_aux, dPdg_lk_aux, dPdbg_lk, dPdg_l
    real matrix dPdg_l_aux, dPdg_kl0, dLdg_kl, dPdg_l0

    // Variables to use in the computation
    // -----------------------------------

    // Get choice sets, xb, and xg + zd
    C      = moptimize_init_userinfo(M, 1) // Choice sets
    PC0    = moptimize_init_userinfo(M, 2) // Importance weights
    R      = moptimize_init_userinfo(M, 3) // Repetitions
    x      = moptimize_init_userinfo(M, 5) // utility vars
    xz     = moptimize_init_userinfo(M, 6) // attention vars
    info   = moptimize_init_userinfo(M, 7) // pos of each group in data
    def    = moptimize_init_userinfo(M, 9) // default good
    xb     = moptimize_util_xb(M, b, 1)    // x * b
    xzd    = moptimize_util_xb(M, b, 2)    // x *g + z * d
    y      = moptimize_util_depvar(M, 1)   // outcome (choice)
    kx     = cols(x)
    kxz    = cols(xz)
    ktot   = kx + kxz

    phi_ij   = invlogit(xzd)              // Get phi for P(attention)
    _editvalue(phi_ij, 1, 1 - epsilon(1)) // Numerical precision foo
    _editvalue(phi_ij, 0, epsilon(1))     // Numerical precision foo

    xdphi_ij = alogit_logisticden(xzd) :* xz // Helper for the Hessian
    xphi_ij  = phi_ij :* xz                  // Helper for the gradient
    exb      = exp(xb)                       // Helper for P(Y)
    lphi_c   = log(phi_ij)                   // Helpers for P(A) since Mata
    lphi_nc  = log(1 :- phi_ij)              // does not have a prod function.

    H  = J(ktot, ktot,  0) // Hessian
    g  = J(1, ktot,  0)    // Gradient
    fv = 0                 // Log likelihood

    // Loop over individuals
    // ---------------------

    for (i = 1; i <= rows(info); i++) {
        y_ij = panelsubmatrix(y, i, info)
        if(sum(y_ij) != 1) continue
        d_ij = selectindex(panelsubmatrix(def, i, info))

        // Grab the sampled choice sets
        Ci = panelsubmatrix(C, i, info)'

        // Compute the likelihood
        // ----------------------

        // Probability of each choice set. We want to compute P(Y |
        // c, theta) P(c | theta); we compute P(c | theta) / (total
        // for P(Y)) and then sum over all the choice sets
        lPA    = Ci * panelsubmatrix(lphi_c, i, info) +
                 (1 :- Ci) * panelsubmatrix(lphi_nc, i, info)
        PAPC0  = rowshape(exp(lPA), 1) :/ PC0[i, ]
        exb_ij = panelsubmatrix(exb, i, info)
        PY     = rowshape(exb_ij, 1) :* (Ci :/ (Ci * exb_ij))
        _editmissing(PY, 0)
        _editmissing(PAPC0, 0)
        PY_ij  = PAPC0 * PY

        if (st_local("noprob") == "error") {
            if (sum(!(PY_ij :> 0)) > 0) {
                err_noprob = "The sampler returned probabilities = 0\n" +
                "To fix this, try\n" +
                "- Increasing 'reps'\n" +
                "- Run with -noprob(drop)- to force 0s in the likelihood.\n"
                errprintf(err_noprob)
                exit()
            }
        }

        // Compute the likelihood
        lPY_ij       = PY_ij / R
        lPY_ij[d_ij] = lPY_ij[d_ij]
        ll = log(lPY_ij)
        _editmissing(ll, 0)
        fv = fv + ll * y_ij

        // Figure out if none of the goods were chosen, and the prob of that
        if(sum(y_ij) == 0) {
            fv = fv + sum(panelsubmatrix(lphi_nc, i, info))
        }

        // Compute the gradient and the Hessian
        // ------------------------------------

        if (todo >= 1) {
            dLdP   = y_ij :/ colshape(PY_ij, 1)
            x_ij   = panelsubmatrix(x,  i, info)
            xz_ij  = panelsubmatrix(xz, i, info)
            xphi_i = panelsubmatrix(xphi_ij, i, info)'
            _editmissing(dLdP, 0)

            // For betas
            for (k = 1; k <= kx; k++) {
                dPdb_k_aux = PY * x_ij[, k]
                dPdb_k = PY_ij :* x_ij[, k]' - ((dPdb_k_aux' :* PAPC0) * PY)
                dLdb_k = dPdb_k * dLdP
                g[, k] = g[, k] + dLdb_k

                if (todo == 2) {
                    for (l = k; l <= kx; l++) {
                        dPdb_l0  = PY :* x_ij[, l]' - PY :* (PY * x_ij[, l])
                        dLdb_kl1 = dPdb_l0 :* x_ij[, k]'
                        dLdb_kl2 = dPdb_l0 :* dPdb_k_aux
                        dLdb_kl3 = PY :* (dPdb_l0 * x_ij[, k])
                        dLdb_kl0 = dLdb_kl1 -  dLdb_kl2 - dLdb_kl3
                        dPdb_kl  = ((PAPC0 * dLdb_kl0) * dLdP)
                        dLdbg_kl = (dPdb_k :* (PAPC0 * dPdb_l0)) * (dLdP:^2)
                        H[k, l]  = H[k, l] + dPdb_kl - dLdbg_kl
                    }

                    // Cross partials
                    dPdb_kl_aux = PY :* x_ij[, k]' - (PY * x_ij[, k]) :* PY
                    for (l = 1; l <= kxz; l++) {
                        dPdg_lk_aux = rowshape((Ci * xz_ij[, l]) :- sum(xphi_i[l, ]), 1) :* PAPC0
                        dPdbg_lk = (dPdg_lk_aux * dPdb_kl_aux) * dLdP
                        dPdg_l   =  dPdg_lk_aux * PY
                        H[k, kx + l] = H[k, kx + l] + dPdbg_lk - (dPdb_k :* dPdg_l) * (dLdP:^2)
                    }
                }
            }

            // For gamma, delta
            haux = uppertriangle(panelsubmatrix(xdphi_ij, i, info)' * xz_ij)
            for (k = 1; k <= kxz; k++) {
                dPdg_k_aux = rowshape((Ci * xz_ij[, k]) :- sum(xphi_i[k, ]), 1)
                dPdg_k = (dPdg_k_aux :* PAPC0) * PY
                g[, kx + k] = g[, kx + k] + dPdg_k * dLdP

                if (todo == 2) {
                    for (l = k; l <= kxz; l++) {
                        dPdg_l_aux = rowshape((Ci * xz_ij[, l]) :- sum(xphi_i[l, ]), 1)
                        dPdg_kl0 = (dPdg_l_aux :* dPdg_k_aux :- haux[k, l]) :* PAPC0
                        dPdg_l0  = (dPdg_l_aux :* PAPC0) * PY
                        dLdg_kl  = (dPdg_kl0 * PY) * dLdP - (dPdg_k :* dPdg_l0) * (dLdP:^2)
                        H[kx + k, kx + l] = H[kx + k, kx + l] + dLdg_kl
                    }
                }
            }
            if(sum(y_ij) == 0) {
                g = g - (J(1, kx, 0), rowsum(xphi_i)')
                if (todo == 2) {
                    H[|kx + 1, kx + 1 \ ., .|] = H[|kx + 1, kx + 1 \ ., .|] - haux
                }
            }
        }
    }
    if (todo == 2) H = makesymmetric(uppertriangle(H)')
}
end
///////////////////////////////////////////////////////////////////////
//                     alogit prediction helpers                     //
///////////////////////////////////////////////////////////////////////

// capture mata: mata drop alogit_predict()
// capture mata: mata drop alogit_counterfactual()
// capture mata: mata drop alogit_predict_loop()
// capture mata: mata drop alogit_counterfactual_loop()

mata

// Prediction (various model fits, elasticity, etc.)
real function alogit_predict(transmorphic M,
                             real matrix y,
                             real matrix x,
                             real matrix xz,
                             real matrix defcons,
                             string scalar method,
                             string scalar predict,
                             | real rowvector b)
{

    // Function variables
    real matrix info, C
    real colvector select, nj_i, nc_i
    real colvector xb, xzd
    real scalar kx, getdef, nc_0, ipos, apos, db, dg, R, dydx, nj
    string scalar afit
    real rowvector b0

    // Function helpers
    real colvector exb, phi_ij, lphi_c, lphi_nc, phi0, Pc0
    real matrix PC0

    // Loop variables to use for the likelihood
    real scalar i, p0
    real matrix Ci, PY
    real colvector y_ij, exb_ij, dc_ij, phi_i
    real colvector lPA, dPYPC
    real rowvector PA, PY_ij, PAPC0, ll, dPY_ij, lPY_ij

    // Variables to return
    real scalar mp0, loglik
    real colvector P_ij, PYC_ij, P0_i, dP_ij

    // Variables to use in the computation
    // -----------------------------------

    moptimize_init_touse(M, st_local("touse"))
    moptimize_init_by(M, st_local("group"))

    info   = panelsetup(*moptimize_util_by(M), 1)
    select = selectindex(st_data(., st_local("touse"), st_local("e(sample)")))

    // Parse vector to use
    if (args() == 7) {
        afit = st_local("afit")
        b0   = moptimize_init_eq_coefs(M, 1), moptimize_init_eq_coefs(M, 2)
        b    = (afit == "loglik")? b0: moptimize_result_coefs(M)
    }
    kx = cols(x)
    if ( moptimize_init_userinfo(M, 8) ) b = log(b)

    // For derivatives, elasticities
    dydx = (predict == "dpycdpc") | (predict == "dpyc") | (predict == "dpc")
    if (dydx) {
        ipos = strtoreal(st_local("ipos"))
        apos = strtoreal(st_local("apos"))
        db   = (ipos == 0)? 1 : b[ipos]
        dg   = (apos == 0)? 1 : b[kx + apos]
    }

    // Linear predictions
    xb  = x  * b[|1 \ kx|]'
    if (predict == "u") return(xb)
    xzd = xz * b[|kx + 1 \ .|]'
    if (predict == "a") return(xzd)

    // P(Attention)
    getdef = (st_local("e(default)") != "") & (st_local("e(consider)") == "")
    getdef = getdef | ((st_local("default") != "") & (st_local("consider") == ""))
    phi_ij = invlogit(xzd)
    if (!getdef) phi_ij[selectindex(defcons), ] = J(sum(defcons), 1, 1)
    if (predict == "pa") return(phi_ij)

    // Helpers for loop
    _editvalue(phi_ij, 1, 1 - epsilon(1))
    _editvalue(phi_ij, 0, epsilon(1))
    exb      = exp(xb)
    lphi_c   = log(phi_ij)
    lphi_nc  = log(1 :- phi_ij)

    // Various parsing for different methods
    if (method == "importance") {
        R    = moptimize_init_userinfo(M, 3)
        phi0 = invlogit(xz * moptimize_init_eq_coefs(M, 2)')
        C    = phi0 :> moptimize_init_userinfo(M, 4)[select, ]
        Pc0  = C :* log(phi0) :+ (1 :- C) :* log(1 :- phi0)
        PC0  = J(rows(info), R, .)
        for (i = 1; i <= rows(info); i++) {
            PC0[i, ] = exp(colsum(panelsubmatrix(Pc0, i, info)))
        }
    }
    else {
        nj   = panelstats(info)[4]
        nj_i = info[, 2] :- info[, 1] :+ 1
        nc_i = 2:^nj_i :- 1
        C    = alogit_combos(nj)
    }

    loglik = mp0 = 0
    P_ij   = PYC_ij = P0_i = dP_ij = J(rows(y), 1, .)
    for (i = 1; i <= rows(info); i++) {
        y_ij = panelsubmatrix(y, i, info)
        if(sum(y_ij) != 1) {
            // if (sum(y_ij) > 1) {
            //     errprintf("dropped group with multiple positive outcomes\n")
            // }
            // if (sum(y_ij) == 0) {
            //     errprintf("dropped group with no positive outcomes\n")
            // }
            continue
        }

        exb_ij = panelsubmatrix(exb, i, info)
        PYC_ij[|info[i, ]'|] = exb_ij / sum(exb_ij)
        if (predict == "pyc") continue

        p0  = exp(sum(panelsubmatrix(lphi_nc, i, info)))
        mp0 = mp0 + p0
        if (predict == "mp0") continue
        P0_i[|info[i, ]'|] = J(length(y_ij), 1, p0)
        if (predict == "p0") continue

        dc_ij  = panelsubmatrix(defcons, i, info)
        if (method == "exact") {
            if (getdef) {
                Ci = C[|1, 1 \ nc_i[i] + 1, nj_i[i]|]
            }
            else {
                nc_0 = nc_i[i] - 2^(nj_i[i] - sum(dc_ij)) + 2
                Ci = J(1, nj_i[i], 0) \ C[|nc_0, 1 \ nc_i[i] + 1, nj_i[i]|]
            }
        }
        else if (method == "importance") {
            Ci = panelsubmatrix(C, i, info)'
        }
        lPA    = Ci * panelsubmatrix(lphi_c, i, info) +
                 (1 :- Ci) * panelsubmatrix(lphi_nc, i, info)
        PA     = rowshape(exp(lPA), 1)
        PY     = rowshape(exb_ij, 1) :* (Ci :/ (Ci * exb_ij))
        _editmissing(PY, 0)
        if (getdef & (method == "exact")) PY[1, selectindex(dc_ij)] = 1

        ///////////////
        //  Margins  //
        ///////////////
        if (predict == "dpycdpc") {
            phi_i = rowshape(panelsubmatrix(phi_ij, i, info), 1)
            dPYPC = (1 :- PY) :* db :+ (Ci :- phi_i) :* dg
        }
        else if (predict == "dpyc") {
            dPYPC = (1 :- PY) :* db
        }
        else if (predict == "dpc") {
            phi_i = rowshape(panelsubmatrix(phi_ij, i, info), 1)
            dPYPC = (Ci :- phi_i) :* dg
        }
        else {
            dPYPC = 1
        }
        ///////////////
        //  Margins  //
        ///////////////

        if (method == "exact") {
            PY_ij  = PA * PY
            dPY_ij = PA * (PY :* dPYPC)
            ll     = log(PY_ij)
        }
        else if (method == "importance")  {
            PAPC0  = PA :/ PC0[i, ]
            _editmissing(PAPC0, 0)
            PY_ij  = (PAPC0 * PY) / R
            dPY_ij = (PAPC0 * (PY :* dPYPC)) / R
            ll     = log(PY_ij)
            _editmissing(ll, 0)
        }

        loglik = loglik + ll * y_ij
        P_ij[|info[i, ]'|]  = colshape(PY_ij, 1)
        dP_ij[|info[i, ]'|] = colshape(dPY_ij, 1)
    }
    mp0 = mp0 / (i - 1)

         if (predict == "py")  return(P_ij)
    else if (predict == "pyc") return(PYC_ij)
    else if (predict == "ll")  return(loglik)
    else if (predict == "p0")  return(P0_i)
    else if (predict == "mp0") return(mp0)
    else if (dydx)             return(dP_ij)
}

// Counterfactual (user-provider likelihood and/or P(A))
real function alogit_counterfactual(real colvector u,
                                    real colvector phi,
                                    real colvector defcons,
                                    real matrix info)
{
    real scalar i, nj, getdef, nc_0
    real matrix C, Ci, PY
    real colvector P_ij, select, nj_i, nc_i
    real colvector exb, lphi_c, lphi_nc
    real colvector u_ij, y_ij, exb_ij, dc_ij, lPA
    real rowvector PA, PY_ij

    getdef  = (st_local("e(default)") != "") & (st_local("e(consider)") == "")
    getdef  = getdef | ((st_local("default") != "") & (st_local("consider") == ""))
    if (!getdef) phi[selectindex(defcons), ] = J(sum(defcons), 1, 1)
    _editvalue(phi, 1, 1 - epsilon(1))
    _editvalue(phi, 0, epsilon(1))

    select  = selectindex(st_data(., st_local("touse")))
    nj      = panelstats(info)[4]
    nj_i    = info[, 2] :- info[, 1] :+ 1
    nc_i    = 2:^nj_i :- 1
    C       = alogit_combos(nj)
    exb     = exp(u)
    lphi_c  = log(phi)
    lphi_nc = log(1 :- phi)
    P_ij    = J(rows(u), 1, .)

    for (i = 1; i <= rows(info); i++) {
        if(!anyof(select, info[i, 1]) & !anyof(select, info[i, 2])) continue
        u_ij  = panelsubmatrix(u, i, info)
        y_ij  = max(u_ij) :== u_ij
        dc_ij = panelsubmatrix(defcons, i, info)
        if (getdef) {
            Ci = C[|1, 1 \ nc_i[i] + 1, nj_i[i]|]
        }
        else {
            nc_0 = nc_i[i] - 2^(nj_i[i] - sum(dc_ij)) + 2
            Ci = C[|nc_0, 1 \ nc_i[i] + 1, nj_i[i]|]
        }

        if(sum(y_ij) != 1) {
            // if(sum(y_ij) > 1) {
            //     errprintf("dropped group with multiple positive outcomes\n")
            // }
            // else if(sum(y_ij) == 0) {
            //     errprintf("dropped group with no positive outcomes\n")
            // }
            continue
        }

        lPA    = Ci * panelsubmatrix(lphi_c, i, info) +
                 (1 :- Ci) * panelsubmatrix(lphi_nc, i, info)
        PA     = rowshape(exp(lPA), 1)
        exb_ij = panelsubmatrix(exb, i, info)
        PY     = rowshape(exb_ij, 1) :* (Ci :/ (Ci * exb_ij))
        _editmissing(PY, 0)
        if (getdef) PY[1, selectindex(dc_ij)] = 1
        PY_ij  = PA * PY
        P_ij[|info[i, ]'|] = colshape(PY_ij, 1)
    }
    return(P_ij)
}

// Predict using loop algorithm
real function alogit_predict_loop(transmorphic M,
                                  real matrix y,
                                  real matrix x,
                                  real matrix xz,
                                  real matrix defcons,
                                  string scalar method,
                                  string scalar predict,
                                  | real rowvector b)
{
    // Function variables
    real matrix info
    real colvector nj_i, xb, xzd
    real scalar kx, getdef, ipos, apos, db, dg, dydx
    string scalar afit
    real rowvector b0

    // Function helpers
    real colvector exb, phi_ij, lphi_c, lphi_nc

    // Loop variables to use for the likelihood
    real scalar i, sel, PA, p0
    real colvector y_ij, exb_ij, phi_i, lphi_c_ij, lphi_nc_ij
    real colvector dPYPC, PY, PY_ijc, PY_ij
    real rowvector c, dPY_ij

    // Variables to return
    real scalar mp0, loglik
    real colvector P_ij, PYC_ij, P0_i, dP_ij

    if (method != "exact") {
        errprintf("-algorithm(loop)- only available for -method(exact)-")
        exit()
    }

    // Variables to use in the computation
    // -----------------------------------

    // Re-initialize info, group if user changed them
    moptimize_init_touse(M, st_local("touse"))
    moptimize_init_by(M, st_local("group"))

    info = panelsetup(*moptimize_util_by(M), 1)
    nj_i = info[, 2] :- info[, 1] :+ 1

    // Use user coefficients if passed
    if (args() == 7) {
        afit = st_local("afit")
        b0   = moptimize_init_eq_coefs(M, 1), moptimize_init_eq_coefs(M, 2)
        b    = (afit == "loglik")? b0: moptimize_result_coefs(M)
    }
    kx  = cols(x)
    if ( moptimize_init_userinfo(M, 8) ) b = log(b)

    // Helper for derivatives/elasticities
    dydx = (predict == "dpycdpc") | (predict == "dpyc") | (predict == "dpc")
    if (dydx) {
        ipos = strtoreal(st_local("ipos"))
        apos = strtoreal(st_local("apos"))
        db   = (ipos == 0)? 1 : b[ipos]
        dg   = (apos == 0)? 1 : b[kx + apos]
    }

    // New fit!
    xb  = x  * b[|1 \ kx|]'
    if (predict == "u") return(xb)
    xzd = xz * b[|kx + 1 \ .|]'
    if (predict == "a") return(xzd)
    phi_ij = invlogit(xzd)
    getdef = (st_local("e(default)") != "") & (st_local("e(consider)") == "")
    getdef = getdef | ((st_local("default") != "") & (st_local("consider") == ""))
    if (!getdef) phi_ij[selectindex(defcons), ] = J(sum(defcons), 1, 1)
    if (predict == "pa") return(phi_ij)

    // For the likelihood, P(A), P(Y)
    _editvalue(phi_ij, 1, 1 - epsilon(1))  // Numerical precision foo
    _editvalue(phi_ij, 0, epsilon(1))      // Numerical precision foo
    exb      = exp(xb)                     // Helper for P(Y)
    lphi_c   = log(phi_ij)                 // Helper for P(A)
    lphi_nc  = log(1 :- phi_ij)            // Helper for P(A)

    // Loop through individuals
    loglik = mp0 = 0
    P_ij   = PYC_ij = P0_i = dP_ij = J(rows(y), 1, .)
    for (i = 1; i <= rows(info); i++) {

        y_ij = panelsubmatrix(y, i, info)
        if(sum(y_ij) != 1) {
            // if(sum(y_ij) > 1) {
            //     errprintf("dropped group with multiple positive outcomes\n")
            // }
            // else if(sum(y_ij) == 0) {
            //     errprintf("dropped group with no positive outcomes\n")
            // }
            continue
        }

        exb_ij = panelsubmatrix(exb, i, info)
        PYC_ij[|info[i, ]'|] = exb_ij / sum(exb_ij)
        if (predict == "pyc") continue

        lphi_c_ij  = panelsubmatrix(lphi_c, i, info)
        lphi_nc_ij = panelsubmatrix(lphi_nc, i, info)

        p0  = exp(sum(lphi_nc_ij))
        mp0 = mp0 + p0
        if (predict == "mp0") continue
        P0_i[|info[i, ]'|] = J(length(y_ij), 1, p0)
        if (predict == "p0") continue

        // Initialize P(Y)
        // ---------------

        sel    = selectindex(y_ij)
        if (getdef) {
            // If you specify a default good, start with the empty set
            c  = J(1, nj_i[i], 0)
            PY = panelsubmatrix(defcons, i, info)
        }
        else {
            // Otherwise, start with the first set that includes the
            // the goods that must be considered.
            c  = panelsubmatrix(defcons, i, info)'
            PY = c' :* exb_ij :/ (c * exb_ij)
        }
        PA     = exp(c * lphi_c_ij + (1 :- c) * lphi_nc_ij)
        PY_ijc = PA :* PY
        PY_ij  = PA * PY[sel]

        ///////////////
        //  Margins  //
        ///////////////
        if (predict == "dpycdpc") {
            phi_i = panelsubmatrix(phi_ij, i, info)
            dPYPC = (1 :- PY) :* db :+ (c' :- phi_i) :* dg
        }
        else if (predict == "dpyc") {
            dPYPC = (1 :- PY) :* db
        }
        else if (predict == "dpc") {
            phi_i = panelsubmatrix(phi_ij, i, info)
            dPYPC = (c' :- phi_i) :* dg
        }
        else {
            dPYPC = 1
        }
        dPY_ij = PA :* PY :* dPYPC
        ///////////////
        //  Margins  //
        ///////////////

        while (!min(c)) {
             _alogit_binadd_reverse(c)
            PA     = exp(c * lphi_c_ij + (1 :- c) * lphi_nc_ij)
            PY     = c' :* exb_ij :/ (c * exb_ij)
            PY_ijc = PY_ijc :+ PA :* PY
            PY_ij  = PY_ij + PA * PY[sel]
            ///////////////
            //  Margins  //
            ///////////////
            if (predict == "dpycdpc") {
                phi_i = panelsubmatrix(phi_ij, i, info)
                dPYPC = (1 :- PY) :* db :+ (c' :- phi_i) :* dg
            }
            else if (predict == "dpyc") {
                dPYPC = (1 :- PY) :* db
            }
            else if (predict == "dpc") {
                phi_i = panelsubmatrix(phi_ij, i, info)
                dPYPC = (c' :- phi_i) :* dg
            }
            else {
                dPYPC = 1
            }
            dPY_ij = dPY_ij :+ PA :* (PY :* dPYPC)
            ///////////////
            //  Margins  //
            ///////////////
        }

        loglik = loglik + log(PY_ij)
        P_ij[|info[i, ]'|]  = PY_ijc
        dP_ij[|info[i, ]'|] = dPY_ij
    }
    mp0 = mp0 / (i - 1)

         if (predict == "py")  return(P_ij)
    else if (predict == "pyc") return(PYC_ij)
    else if (predict == "ll")  return(loglik)
    else if (predict == "p0")  return(P0_i)
    else if (predict == "mp0") return(mp0)
    else if (dydx)             return(dP_ij)
}



// Counterfactual using loop algorithm
real colvector function alogit_counterfactual_loop(real colvector u,
                                                   real colvector phi,
                                                   real colvector defcons,
                                                   real matrix info)
{
    real scalar i, getdef, sel, PA
    real colvector P_ij, select, nj_i, PY, PY_ijc
    real colvector exb, lphi_c, lphi_nc, lphi_c_ij, lphi_nc_ij
    real colvector u_ij, y_ij, exb_ij
    real rowvector c

    getdef  = (st_local("e(default)") != "") & (st_local("e(consider)") == "")
    getdef  = getdef | ((st_local("default") != "") & (st_local("consider") == ""))
    if (!getdef) phi[selectindex(defcons), ] = J(sum(defcons), 1, 1)
    _editvalue(phi, 1, 1 - epsilon(1))
    _editvalue(phi, 0, epsilon(1))

    select  = selectindex(st_data(., st_local("touse")))
    nj_i    = info[, 2] :- info[, 1] :+ 1
    exb     = exp(u)
    lphi_c  = log(phi)
    lphi_nc = log(1 :- phi)
    P_ij    = J(rows(u), 1, .)

    for (i = 1; i <= rows(info); i++) {
        if(!anyof(select, info[i, 1]) & !anyof(select, info[i, 2])) continue
        u_ij = panelsubmatrix(u, i, info)
        y_ij = max(u_ij) :== u_ij

        if(sum(y_ij) != 1) {
            // if(sum(y_ij) > 1) {
            //     errprintf("dropped group with multiple positive outcomes\n")
            // }
            // else if(sum(y_ij) == 0) {
            //     errprintf("dropped group with no positive outcomes\n")
            // }
            continue
        }

        lphi_c_ij  = panelsubmatrix(lphi_c, i, info)
        lphi_nc_ij = panelsubmatrix(lphi_nc, i, info)
        exb_ij     = panelsubmatrix(exb, i, info)

        // Initialize P(Y)
        // ---------------

        sel    = selectindex(y_ij)
        if (getdef) {
            // If you specify a default good, start with the empty set
            c  = J(1, nj_i[i], 0)
            PY = panelsubmatrix(defcons, i, info)
        }
        else {
            // Otherwise, start with the first set that includes the
            // goods that must be considered.
            c  = panelsubmatrix(defcons, i, info)'
            PY = c' :* exb_ij :/ (c * exb_ij)
        }
        PA     = exp(c * lphi_c_ij + (1 :- c) * lphi_nc_ij)
        PY_ijc = PA :* PY

        // Loop through all consideration sets
        // -----------------------------------

        while (!min(c)) {
             _alogit_binadd_reverse(c)
            PA     = exp(c * lphi_c_ij + (1 :- c) * lphi_nc_ij)
            PY     = c' :* exb_ij :/ (c * exb_ij)
            PY_ijc = PY_ijc :+ PA :* PY
        }

        P_ij[|info[i, ]'|] = PY_ijc
    }
    return(P_ij)
}
end
///////////////////////////////////////////////////////////////////////
//                           DSC functions                           //
///////////////////////////////////////////////////////////////////////

// capture mata: mata drop _alogit_dsc_exact()
// capture mata: mata drop _alogit_dsc_exact_fast()
// capture mata: mata drop alogit_dsc_predict()
// capture mata: mata drop alogit_dsc_counterfactual()

mata

// Loglik for DSC model, only one loop
void function _alogit_dsc_exact(transmorphic M,
                                real scalar todo,
                                real rowvector b,
                                fv, g, H)
{
    // Function variables
    real matrix x, xz, info, xbmat, xzdmat
    real colvector y, def, xb, xzd, d_sel
    real rowvector bb, gd
    real scalar asexp, kx, kxz, ktot

    // Function helpers
    real matrix xdphi_ij
    real colvector exb, pa_ij, pi_ij

    // Loop variables to use for the likelihood
    real scalar i, k, sel, PA
    real colvector y_ij, exb_ij, d_ij, PY, PY_ij, P_ij

    // Loop variables to use for the gradient
    real scalar index, dLdP
    real rowvector dLdb_k, dLdg_k, dLdb_k_aux
    real matrix x_ij, xz_ij, xdphi_i

    // Loop variables to use for the Hessian
    real rowvector PYx_ij
    real matrix hb, hg, hbg
    real matrix dPdb_kl1_1, dPdb_kl2_1, dPdb_kl2_2, dPdb_kl
    real matrix dPdbg_kl
    real matrix dPdgd_aux, dPdgd_kl

    // Variables to use in the computation
    // -----------------------------------

    // Get choice sets, xb, and xg + zd
    x     = moptimize_init_userinfo(M, 5) // utility vars
    xz    = moptimize_init_userinfo(M, 6) // attention vars
    info  = moptimize_init_userinfo(M, 7) // pos of each group in data
    asexp = moptimize_init_userinfo(M, 8) // Estimate exp(b) isntead of b
    def   = moptimize_init_userinfo(M, 9) // default good
    y     = moptimize_util_depvar(M, 1)   // outcome (choice)
    kx    = cols(x)
    kxz   = cols(xz)
    ktot  = kx + kxz

    if ( asexp ) {
        bb      = b[|1 \ kx|]      // beta
        gd      = b[|kx + 1 \ .|]  // gamma, delta
        xbmat   = (bb :^ x)        // exp(b)^x * b
        xzdmat  = (gd :^ xz)       // exp(g, d)^(x, z)
        xb      = xbmat[, 1]       // prod(exp(b)^x over kx) = exp(xb)
        xzd     = xzdmat[, 1]      // prod(exp(g, d)^(x, z) over kx + kz) = exp(xg + zd)
        for (k = 2; k <= kx; k++)  xb  = xb  :* xbmat[, k]
        for (k = 2; k <= kxz; k++) xzd = xzd :* xzdmat[, k]
    }
    else {
        xb  = moptimize_util_xb(M, b, 1) // x * b
        xzd = moptimize_util_xb(M, b, 2) // x * g + z * d
    }

    d_sel = selectindex(def)  // Index of default goods
    xzd   = xzd[d_sel]        // P(A) only depends on default good

    // Get phi for P(attention)
    pa_ij = asexp? (xzd :/ (1 :+ xzd)): invlogit(xzd)

    _editvalue(pa_ij, 1, 1 - epsilon(1)) // Numerical precision foo
    _editvalue(pa_ij, 0, epsilon(1))     // Numerical precision foo

    pi_ij = 1 :- pa_ij         // Get phi for P(inattention)
    exb   = asexp? xb: exp(xb) // Helper for P(Y)

    // Helpers for the Hessian, gradient
    xdphi_ij = ( asexp? (xzd :/ (1 :+ xzd):^2): alogit_logisticden(xzd) ) :* xz[d_sel, ]

    H  = J(ktot, ktot,  0)  // Hessian
    g  = J(1, ktot,  0)     // Gradient
    fv = 0                  // Log likelihood

    // Loop over individuals
    // ---------------------

    for (i = 1; i <= rows(info); i++) {
        // Drop individuals with > 1 or 0 goods chosen
        y_ij = panelsubmatrix(y, i, info)
        if(sum(y_ij) != 1) continue

        // Compute the likelihood
        // ----------------------

        // Probability of each choice set. We want to compute P(Y |
        // c, theta) P(c | theta); we compute P(c | theta) / (total
        // for P(Y)) and then sum over all the choice sets
        sel    = selectindex (y_ij)
        d_ij   = panelsubmatrix(def, i, info)
        exb_ij = panelsubmatrix(exb, i, info)
        PY     = exb_ij :/ sum(exb_ij)
        PA     = pa_ij[i]
        PY_ij  = PA :* PY
        P_ij   = PY_ij + d_ij :* pi_ij[i]

        // Compute the likelihood
        fv = fv + log(P_ij[sel])

        // Compute the gradient and the Hessian
        // ------------------------------------

        if (todo >= 1) {
            index   = selectindex (d_ij)
            dLdP    = 1 / P_ij[sel]
            x_ij    = panelsubmatrix(x, i, info)
            xdphi_i = xdphi_ij[i, ]

            // Gradient
            dLdb_k_aux = (x_ij[sel, ] :- PY' * x_ij)
            dLdb_k = dLdb_k_aux * PY_ij[sel] * dLdP
            dLdg_k = xdphi_i :* (PY[sel] - (sel == index)) * dLdP

            g = g :+ (dLdb_k, dLdg_k)

            if (todo == 2) {
                xz_ij = xz[d_sel[i], ]

                // Hessian, beta
                PYx_ij     = PY' * x_ij
                dPdb_kl1_1 = dLdb_k_aux' * dLdb_k_aux
                dPdb_kl2_1 = PYx_ij' * PYx_ij
                dPdb_kl2_2 = x_ij' * diag(PY) * x_ij
                dPdb_kl    = (dPdb_kl1_1 - (dPdb_kl2_2 - dPdb_kl2_1)) * PY_ij[sel]

                // Hessian, cross
                dPdbg_kl   = PY[sel] * (dLdb_k_aux' * xdphi_i)

                // Hessian, gamma, delta
                dPdgd_aux  = (1 - 2 * PA) * (PY[sel] - (sel == index))
                dPdgd_kl   = dPdgd_aux * (xdphi_i' * xz_ij)

                // Hessian
                hb  = dLdP :* dPdb_kl  :- dLdb_k' * dLdb_k
                hg  = dLdP :* dPdgd_kl :- dLdg_k' * dLdg_k
                hbg = dLdP :* dPdbg_kl :- dLdb_k' * dLdg_k
                H = H :+ ((hb, hbg) \ (hbg', hg))
            }
        }
    }
    if (todo == 2) H = makesymmetric(H')
    if ( asexp ) {
        if (todo >= 1) g = g :/ b
        if (todo == 2) H = H :/ (b' * b) :- diag(g :/ b)
    }
}

void function _alogit_dsc_exact_fast(transmorphic M,
                                     real scalar todo,
                                     real rowvector b,
                                     fv, g, H)
{
    real scalar ktot
    ktot = st_numscalar("__alogit_ktot")
    H    = J(ktot, ktot,  0);
    g    = J(1, ktot,  0);
    st_numscalar("__alogit_ll", .)
    st_matrix("__alogit_b", b)
    st_numscalar("__alogit_todo", todo)
    if (moptimize_init_userinfo(M, 1) == "fast") {
        stata("plugin call aloglik_fast " + st_local("fastcall") + " dsc")
    }
    else {
        stata("plugin call aloglik_faster " + st_local("fastcall") + " dsc")
    }
    if (todo >= 1) g = st_matrix("__alogit_g")
    if (todo == 2) H = makesymmetric(uppertriangle(st_matrix("__alogit_H"))')
    fv = st_numscalar("__alogit_ll")
    if ( moptimize_init_userinfo(M, 8) ) {
        if (todo >= 1) g = g :/ b
        if (todo == 2) H = H :/ (b' * b) :- diag(g :/ b)
    }
}

///////////////////////////////////////////////////////////////////////
//                      DSC prediction helpers                       //
///////////////////////////////////////////////////////////////////////

real function alogit_dsc_predict(transmorphic M,
                                 real matrix y,
                                 real matrix x,
                                 real matrix xz,
                                 real matrix def,
                                 string scalar method,
                                 string scalar predict,
                                 | real rowvector b)
{

    // Function variables
    real matrix info
    real colvector xb, xzd, d_sel
    real scalar kx, ipos, apos, db, dg, dydx
    string scalar afit
    real rowvector b0

    // Function helpers
    real colvector exb, phi_ij, pa_ij, pi_ij

    // Loop variables to use for the likelihood
    real scalar i, PA, PI
    real colvector y_ij, exb_ij, d_ij
    real colvector dPYPC, dPYPC0, PY, PY_ij
    real rowvector dPY_ij

    // Variables to return
    real scalar mp0, mp1, loglik
    real colvector P_ij, PYC_ij, P0_i, dP_ij

    if (method == "importance") {
        errprintf("-method(importance)- is not available for -model(dsc)-\n")
        exit()
    }

    // Variables to use in the computation
    // -----------------------------------

    // Re-initialize info, group if user changed them
    moptimize_init_touse(M, st_local("touse"))
    moptimize_init_by(M, st_local("group"))

    info = panelsetup(*moptimize_util_by(M), 1)

    // Use user coefficients if passed
    if (args() == 7) {
        afit = st_local("afit")
        b0   = moptimize_init_eq_coefs(M, 1), moptimize_init_eq_coefs(M, 2)
        b    = (afit == "loglik")? b0: moptimize_result_coefs(M)
    }
    kx  = cols(x)
    if ( moptimize_init_userinfo(M, 8) ) b = log(b)

    // Helper for derivatives/elasticities
    dydx = (predict == "dpycdpc") | (predict == "dpyc") | (predict == "dpc")
    if (dydx) {
        ipos = strtoreal(st_local("ipos"))
        apos = strtoreal(st_local("apos"))
        db   = (ipos == 0)? 1 : b[ipos]
        dg   = (apos == 0)? 1 : b[kx + apos]
    }

    // New fit!
    xb  = x  * b[|1 \ kx|]'
    if (predict == "u") return(xb)
    xzd = xz * b[|kx + 1 \ .|]'
    if (predict == "a") return(xzd)

    // For the likelihood, P(A), P(Y)
    d_sel = selectindex(def)
    pa_ij = invlogit(xzd[d_sel])
    _editvalue(pa_ij, 1, 1 - epsilon(1))
    _editvalue(pa_ij, 0, epsilon(1))
    pi_ij = 1 :- pa_ij
    exb   = exp(xb)

    // Loop through individuals
    loglik = mp0 = mp1 = 0
    P_ij   = PYC_ij = P0_i = dP_ij = phi_ij = J(rows(y), 1, .)
    for (i = 1; i <= rows(info); i++) {
        y_ij = panelsubmatrix(y, i, info)
        if(sum(y_ij) != 1) {
            // if(sum(y_ij) > 1) {
            //     errprintf("dropped group with multiple positive outcomes\n")
            // }
            // else if(sum(y_ij) == 0) {
            //     errprintf("dropped group with no positive outcomes\n")
            // }
            continue
        }

        exb_ij = panelsubmatrix(exb, i, info)
        PY     = exb_ij :/ sum(exb_ij)
        PYC_ij[|info[i, ]'|] = PY
        if (predict == "pyc") continue

        d_ij = panelsubmatrix(def, i, info)
        PI   = pi_ij[i]
        PA   = pa_ij[i]

        P0_i[|info[i, ]'|] = J(length(y_ij), 1, PI)
        if (predict == "p0") continue

        phi_ij[|info[i, ]'|] = J(rows(y_ij), 1, PA)
        if (predict == "pa") continue

        mp0 = mp0 + PI
        if (predict == "mp0") continue

        mp1 = mp1 + PA
        if (predict == "mp1") continue

        _editmissing(PY, 0)
        PY_ij = PA :* PY + d_ij :* PI
        P_ij[|info[i, ]'|] = PY_ij
        if (predict == "py") continue

        loglik = loglik + log(PY_ij[selectindex(y_ij)])
        if (predict == "ll") continue

        ///////////////
        //  Margins  //
        ///////////////
        if (predict == "dpycdpc") {
            dPYPC0 = (PA * PI * dg) :* (PY :- 1)
            dPYPC  = (1 :- PY) :* db
        }
        else if (predict == "dpyc") {
            dPYPC0 = 0
            dPYPC  = (1 :- PY) :* db
        }
        else if (predict == "dpc") {
            dPYPC0 = (PA * PI * dg) :* (PY :- 1)
            dPYPC  = 0
        }
        else {
            dPYPC0 = 0
            dPYPC  = 1
        }
        ///////////////
        //  Margins  //
        ///////////////

        dPY_ij = (PA :* PY :* dPYPC) + d_ij :* dPYPC0
        dP_ij[|info[i, ]'|]  = dPY_ij
    }
    mp0 = mp0 / (i - 1)
    mp1 = mp1 / (i - 1)

         if (predict == "py")  return(P_ij)
    else if (predict == "pyc") return(PYC_ij)
    else if (predict == "pa")  return(phi_ij)
    else if (predict == "ll")  return(loglik)
    else if (predict == "p0")  return(P0_i)
    else if (predict == "mp0") return(mp0)
    else if (predict == "mp1") return(mp1)
    else if (dydx)             return(dP_ij)
}

real colvector function alogit_dsc_counterfactual(real colvector u,
                                                  real colvector phi,
                                                  real colvector def,
                                                  real matrix info)
{
    real scalar i, PA, PI
    real colvector P_ij, select, PY
    real colvector exb, pa_ij, PY_ij
    real colvector u_ij, y_ij, exb_ij, d_ij

    select = selectindex(st_data(., st_local("touse")))
    exb    = exp(u)
    pa_ij  = phi[selectindex(def)]
    _editvalue(pa_ij, 1, 1 - epsilon(1))
    _editvalue(pa_ij, 0, epsilon(1))
    P_ij   = J(rows(u), 1, .)

    for (i = 1; i <= rows(info); i++) {
        if(!anyof(select, info[i, 1]) & !anyof(select, info[i, 2])) continue
        u_ij = panelsubmatrix(u, i, info)
        y_ij = max(u_ij) :== u_ij

        if(sum(y_ij) != 1) {
            // if(sum(y_ij) > 1) {
            //     errprintf("dropped group with multiple positive outcomes\n")
            // }
            // else if(sum(y_ij) == 0) {
            //     errprintf("dropped group with no positive outcomes\n")
            // }
            continue
        }
        d_ij   = panelsubmatrix(def, i, info)
        PI     = 1 - pa_ij[i]
        PA     = pa_ij[i]
        exb_ij = panelsubmatrix(exb, i, info)
        PY     = exb_ij :/ sum(exb_ij)
        _editmissing(PY, 0)
        PY_ij  = PA :* PY
        P_ij[|info[i, ]'|] = PY_ij + d_ij :* PI
    }
    return(P_ij)
}
end