text:
"""
```
eqcond(m::Model1010depo_alt)
```

Expresses the equilibrium conditions in canonical form using Γ0, Γ1, C, Ψ, and Π matrices.
Using the mappings of states/equations to integers defined in m1010.jl, coefficients are
specified in their proper positions.

### Outputs

* `Γ0` (`n_states` x `n_states`) holds coefficients of current time states.
* `Γ1` (`n_states` x `n_states`) holds coefficients of lagged states.
* `C`  (`n_states` x `1`) is a vector of constants
* `Ψ`  (`n_states` x `n_shocks_exogenous`) holds coefficients of iid shocks.
* `Π`  (`n_states` x `n_states_expectational`) holds coefficients of expectational states.
"""
function eqcond(m::Model1010depo_alt)
    endo = m.endogenous_states
    exo  = m.exogenous_shocks
    ex   = m.expected_shocks
    eq   = m.equilibrium_conditions

    Γ0 = zeros(n_states(m), n_states(m))
    Γ1 = zeros(n_states(m), n_states(m))
    C  = zeros(n_states(m))
    Ψ  = zeros(n_states(m), n_shocks_exogenous(m))
    Π  = zeros(n_states(m), n_shocks_expectational(m))

    ### ENDOGENOUS STATES ###

    ### 1. Consumption Euler Equation

    # Sticky prices and wages
    Γ0[eq[:eq_euler], endo[:c_t]]  = 1.
    Γ0[eq[:eq_euler], endo[:R_t]]  = (1 - m[:h]*exp(-m[:z_star]))/(m[:σ_c]*(1 + m[:h]*exp(-m[:z_star])))
    Γ0[eq[:eq_euler], endo[:b_t]]  = -1.
    Γ0[eq[:eq_euler], endo[:Eπ_t]] = -(1 - m[:h]*exp(-m[:z_star]))/(m[:σ_c]*(1 + m[:h]*exp(-m[:z_star])))
    Γ0[eq[:eq_euler], endo[:z_t]]  = (m[:h]*exp(-m[:z_star]))/(1 + m[:h]*exp(-m[:z_star]))
    Γ0[eq[:eq_euler], endo[:Ec_t]] = -1/(1 + m[:h]*exp(-m[:z_star]))
    Γ0[eq[:eq_euler], endo[:Ez_t]] = -1/(1 + m[:h]*exp(-m[:z_star]))
    Γ0[eq[:eq_euler], endo[:L_t]]  = -(m[:σ_c] - 1)*m[:wl_c]/(m[:σ_c]*(1 + m[:h]*exp(-m[:z_star])))
    Γ0[eq[:eq_euler], endo[:EL_t]] = (m[:σ_c] - 1)*m[:wl_c]/(m[:σ_c]*(1 + m[:h]*exp(-m[:z_star])))
    Γ1[eq[:eq_euler], endo[:c_t]]  = (m[:h]*exp(-m[:z_star]))/(1 + m[:h]*exp(-m[:z_star]))

    # Flexible prices and wages
    Γ0[eq[:eq_euler_f], endo[:c_f_t]]  = 1.
    Γ0[eq[:eq_euler_f], endo[:r_f_t]]  = (1 - m[:h]*exp(-m[:z_star]))/(m[:σ_c]*(1 + m[:h]*exp(-m[:z_star])))
    Γ0[eq[:eq_euler_f], endo[:b_t]]    = -1.
    Γ0[eq[:eq_euler_f], endo[:z_t]]    = (m[:h]*exp(-m[:z_star]))/(1 + m[:h]*exp(-m[:z_star]))
    Γ0[eq[:eq_euler_f], endo[:Ec_f_t]] = -1/(1 + m[:h]*exp(-m[:z_star]))
    Γ0[eq[:eq_euler_f], endo[:Ez_t]]   = -1/(1 + m[:h]*exp(-m[:z_star]))
    Γ0[eq[:eq_euler_f], endo[:L_f_t]]  = -(m[:σ_c] - 1)*m[:wl_c]/(m[:σ_c]*(1 + m[:h]*exp(-m[:z_star])))
    Γ0[eq[:eq_euler_f], endo[:EL_f_t]] = (m[:σ_c] - 1)*m[:wl_c]/(m[:σ_c]*(1 + m[:h]*exp(-m[:z_star])))
    Γ1[eq[:eq_euler_f], endo[:c_f_t]]  = (m[:h]*exp(-m[:z_star]))/(1 + m[:h]*exp(-m[:z_star]))

    ### 2. Investment Euler Equation

    # Sticky prices and wages
    Γ0[eq[:eq_inv], endo[:qk_t]] = -1/(m[:S′′]*exp(2.0*m[:z_star])*(1 + m[:β]*exp((1 - m[:σ_c])*m[:z_star])))
    Γ0[eq[:eq_inv], endo[:i_t]]  = 1.
    Γ0[eq[:eq_inv], endo[:z_t]]  = 1/(1 + m[:β]*exp((1 - m[:σ_c])*m[:z_star]))
    Γ1[eq[:eq_inv], endo[:i_t]]  = 1/(1 + m[:β]*exp((1 - m[:σ_c])*m[:z_star]))
    Γ0[eq[:eq_inv], endo[:Ei_t]] = -m[:β]*exp((1 - m[:σ_c])*m[:z_star])/(1 + m[:β]*exp((1 - m[:σ_c])*m[:z_star]))
    Γ0[eq[:eq_inv], endo[:Ez_t]] = -m[:β]*exp((1 - m[:σ_c])*m[:z_star])/(1 + m[:β]*exp((1 - m[:σ_c])*m[:z_star]))
    Γ0[eq[:eq_inv], endo[:μ_t]]  = -1.

    # Flexible prices and wages
    Γ0[eq[:eq_inv_f], endo[:qk_f_t]] = -1/(m[:S′′]*exp(2*m[:z_star])*(1 + m[:β]*exp((1 - m[:σ_c])*m[:z_star])))
    Γ0[eq[:eq_inv_f], endo[:i_f_t]]  = 1.
    Γ0[eq[:eq_inv_f], endo[:z_t]]    = 1/(1 + m[:β]*exp((1 - m[:σ_c])*m[:z_star]))
    Γ1[eq[:eq_inv_f], endo[:i_f_t]]  = 1/(1 + m[:β]*exp((1 - m[:σ_c])*m[:z_star]))
    Γ0[eq[:eq_inv_f], endo[:Ei_f_t]] = -m[:β]*exp((1 - m[:σ_c])*m[:z_star])/(1 + m[:β]*exp((1 - m[:σ_c])*m[:z_star]))
    Γ0[eq[:eq_inv_f], endo[:Ez_t]]   = -m[:β]*exp((1 - m[:σ_c])*m[:z_star])/(1 + m[:β]*exp((1 - m[:σ_c])*m[:z_star]))
    Γ0[eq[:eq_inv_f], endo[:μ_t]]    = -1.

    ### 3. Financial Friction Block

    # Return to capital
    # Sticky prices and wages
    Γ0[eq[:eq_capval], endo[:Rktil_t]] = 1.
    Γ0[eq[:eq_capval], endo[:π_t]]      = -1.
    Γ0[eq[:eq_capval], endo[:rk_t]]     = -m[:r_k_star]/(1 + m[:r_k_star] - m[:δ])
    Γ0[eq[:eq_capval], endo[:qk_t]]     = -(1 - m[:δ])/(1 + m[:r_k_star] - m[:δ])
    Γ1[eq[:eq_capval], endo[:qk_t]]     = -1.

    # Spreads
    # Sticky prices and wages
    Γ0[eq[:eq_spread], endo[:ERtil_k_t]] = 1.
    Γ0[eq[:eq_spread], endo[:R_t]]       = -1.
    Γ0[eq[:eq_spread], endo[:b_t]]   = (m[:σ_c]*(1 + m[:h]*exp(-m[:z_star])))/(1 - m[:h]*exp(-m[:z_star]))
    Γ0[eq[:eq_spread], endo[:qk_t]]      = -m[:ζ_spb]
    Γ0[eq[:eq_spread], endo[:kbar_t]]    = -m[:ζ_spb]
    Γ0[eq[:eq_spread], endo[:n_t]]       = m[:ζ_spb]
    Γ0[eq[:eq_spread], endo[:σ_ω_t]]     = -1.
    Γ0[eq[:eq_spread], endo[:μ_e_t]]     = -1.

    # Flexible prices and wages
    Γ0[eq[:eq_spread_f], endo[:ERktil_f_t]] = 1.
    Γ0[eq[:eq_spread_f], endo[:r_f_t]]      = -1.
    Γ0[eq[:eq_spread_f], endo[:b_t]]   = (m[:σ_c]*(1 + m[:h]*exp(-m[:z_star])))/(1 - m[:h]*exp(-m[:z_star]))
    Γ0[eq[:eq_spread_f], endo[:qk_f_t]]      = -m[:ζ_spb]
    Γ0[eq[:eq_spread_f], endo[:kbar_f_t]]    = -m[:ζ_spb]
    Γ0[eq[:eq_spread_f], endo[:n_f_t]]       = m[:ζ_spb]
    Γ0[eq[:eq_spread_f], endo[:σ_ω_t]]     = -1.
    Γ0[eq[:eq_spread_f], endo[:μ_e_t]]     = -1.

    # n evol
    # Sticky prices and wages
    Γ0[eq[:eq_nevol], endo[:n_t]]      = 1.
    Γ0[eq[:eq_nevol], endo[:γ_t]]      = -1.
    Γ0[eq[:eq_nevol], endo[:z_t]]      = m[:γ_star]*m[:vstar]/m[:nstar]
    Γ0[eq[:eq_nevol], endo[:Rktil_t]] = -m[:ζ_nRk]
    Γ0[eq[:eq_nevol], endo[:π_t]]      = (m[:ζ_nRk] - m[:ζ_nR])
    Γ1[eq[:eq_nevol], endo[:σ_ω_t]]    = -m[:ζ_nσ_ω]/m[:ζ_spσ_ω]
    Γ1[eq[:eq_nevol], endo[:μ_e_t]]    = -m[:ζ_nμ_e]/m[:ζ_spμ_e]
    Γ1[eq[:eq_nevol], endo[:qk_t]]     = m[:ζ_nqk]
    Γ1[eq[:eq_nevol], endo[:kbar_t]]   = m[:ζ_nqk]
    Γ1[eq[:eq_nevol], endo[:n_t]]      = m[:ζ_nn]
    Γ1[eq[:eq_nevol], endo[:R_t]]      = -m[:ζ_nR]
    Γ1[eq[:eq_nevol], endo[:b_t]]  = m[:ζ_nR]*((m[:σ_c]*(1.0+m[:h]*exp(-m[:z_star])))/(1.0-m[:h]*exp(-m[:z_star])))

    # Flexible prices and wages
    Γ0[eq[:eq_nevol_f], endo[:n_f_t]]      = 1.
    Γ0[eq[:eq_nevol_f], endo[:z_t]]      = m[:γ_star]*m[:vstar]/m[:nstar]
    Γ0[eq[:eq_nevol_f], endo[:rktil_f_t]] = -m[:ζ_nRk]
    Γ1[eq[:eq_nevol_f], endo[:σ_ω_t]]    = -m[:ζ_nσ_ω]/m[:ζ_spσ_ω]
    Γ1[eq[:eq_nevol_f], endo[:μ_e_t]]    = -m[:ζ_nμ_e]/m[:ζ_spμ_e]
    Γ1[eq[:eq_nevol_f], endo[:qk_f_t]]     = m[:ζ_nqk]
    Γ1[eq[:eq_nevol_f], endo[:kbar_f_t]]   = m[:ζ_nqk]
    Γ1[eq[:eq_nevol_f], endo[:n_f_t]]      = m[:ζ_nn]
    Γ1[eq[:eq_nevol_f], endo[:r_f_t]]      = -m[:ζ_nR]
    Γ1[eq[:eq_nevol_f], endo[:b_t]]  = m[:ζ_nR]*((m[:σ_c]*(1.0+m[:h]*exp(-m[:z_star])))/(1.0-m[:h]*exp(-m[:z_star])))

    # Flexible prices and wages - ASSUME NO FINANCIAL FRICTIONS
    Γ0[eq[:eq_capval_f], endo[:rktil_f_t]] = 1.
    Γ0[eq[:eq_capval_f], endo[:rk_f_t]]     = -m[:r_k_star]/(m[:r_k_star]+1-m[:δ])
    Γ0[eq[:eq_capval_f], endo[:qk_f_t]]     = -(1-m[:δ])/(m[:r_k_star]+1-m[:δ])
    Γ1[eq[:eq_capval_f], endo[:qk_f_t]]     = -1.

    ### 4. Aggregate Production Function

    # Sticky prices and wages
    Γ0[eq[:eq_output], endo[:y_t]] =  1.
    Γ0[eq[:eq_output], endo[:k_t]] = -m[:Φ]*m[:α]
    Γ0[eq[:eq_output], endo[:L_t]] = -m[:Φ]*(1 - m[:α])

    # Flexible prices and wages
    Γ0[eq[:eq_output_f], endo[:y_f_t]] =  1.
    Γ0[eq[:eq_output_f], endo[:k_f_t]] = -m[:Φ]*m[:α]
    Γ0[eq[:eq_output_f], endo[:L_f_t]] = -m[:Φ]*(1 - m[:α])

    ### 5. Capital Utilization

    # Sticky prices and wages
    Γ0[eq[:eq_caputl], endo[:k_t]]    =  1.
    Γ1[eq[:eq_caputl], endo[:kbar_t]] =  1.
    Γ0[eq[:eq_caputl], endo[:z_t]]    = 1.
    Γ0[eq[:eq_caputl], endo[:u_t]]    = -1.

    # Flexible prices and wages
    Γ0[eq[:eq_caputl_f], endo[:k_f_t]]    =  1.
    Γ1[eq[:eq_caputl_f], endo[:kbar_f_t]] =  1.
    Γ0[eq[:eq_caputl_f], endo[:z_t]]      = 1.
    Γ0[eq[:eq_caputl_f], endo[:u_f_t]]    = -1.

    ### 6. Rental Rate of Capital

    # Sticky prices and wages
    Γ0[eq[:eq_capsrv], endo[:u_t]]  = 1.
    Γ0[eq[:eq_capsrv], endo[:rk_t]] = -(1 - m[:ppsi])/m[:ppsi]

    # Flexible prices and wages
    Γ0[eq[:eq_capsrv_f], endo[:u_f_t]]  = 1.
    Γ0[eq[:eq_capsrv_f], endo[:rk_f_t]] = -(1 - m[:ppsi])/m[:ppsi]

    ### 7. Evolution of Capital

    # Sticky prices and wages
    Γ0[eq[:eq_capev], endo[:kbar_t]] = 1.
    Γ1[eq[:eq_capev], endo[:kbar_t]] = 1 - m[:istar]/m[:kbarstar]
    Γ0[eq[:eq_capev], endo[:z_t]]    = 1 - m[:istar]/m[:kbarstar]
    Γ0[eq[:eq_capev], endo[:i_t]]    = -m[:istar]/m[:kbarstar]
    Γ0[eq[:eq_capev], endo[:μ_t]]    = -m[:istar]*m[:S′′]*exp(2*m[:z_star])*(1 + m[:β]*exp((1 - m[:σ_c])*m[:z_star]))/m[:kbarstar]

    # Flexible prices and wages
    Γ0[eq[:eq_capev_f], endo[:kbar_f_t]] = 1.
    Γ1[eq[:eq_capev_f], endo[:kbar_f_t]] = 1 - m[:istar]/m[:kbarstar]
    Γ0[eq[:eq_capev_f], endo[:z_t]]      = 1 - m[:istar]/m[:kbarstar]
    Γ0[eq[:eq_capev_f], endo[:i_f_t]]    = -m[:istar]/m[:kbarstar]
    Γ0[eq[:eq_capev_f], endo[:μ_t]]      = -m[:istar]*m[:S′′]*exp(2*m[:z_star])*(1 + m[:β]*exp((1 - m[:σ_c])*m[:z_star]))/m[:kbarstar]

    ### 8. Price Markup

    # Sticky prices and wages
    Γ0[eq[:eq_mkupp], endo[:mc_t]] =  1.
    Γ0[eq[:eq_mkupp], endo[:w_t]]  = -1.
    Γ0[eq[:eq_mkupp], endo[:L_t]]  = -m[:α]
    Γ0[eq[:eq_mkupp], endo[:k_t]]  =  m[:α]

    # Flexible prices and wages
    Γ0[eq[:eq_mkupp_f], endo[:w_f_t]] = 1.
    Γ0[eq[:eq_mkupp_f], endo[:L_f_t]] =  m[:α]
    Γ0[eq[:eq_mkupp_f], endo[:k_f_t]] =  -m[:α]

    ### 9. Phillips Curve

    # Sticky prices and wages
    Γ0[eq[:eq_phlps], endo[:π_t]]  = 1.
    Γ0[eq[:eq_phlps], endo[:mc_t]] =  -((1 - m[:ζ_p]*m[:β]*exp((1 - m[:σ_c])*m[:z_star]))*
        (1 - m[:ζ_p]))/(m[:ζ_p]*((m[:Φ]- 1)*m[:ϵ_p] + 1))/(1 + m[:ι_p]*m[:β]*exp((1 - m[:σ_c])*m[:z_star]))
    Γ1[eq[:eq_phlps], endo[:π_t]]  = m[:ι_p]/(1 + m[:ι_p]*m[:β]*exp((1 - m[:σ_c])*m[:z_star]))
    Γ0[eq[:eq_phlps], endo[:Eπ_t]] = -m[:β]*exp((1 - m[:σ_c])*m[:z_star])/(1 + m[:ι_p]*m[:β]*
        exp((1 - m[:σ_c])*m[:z_star]))

    # Comment out for counterfactual with no price mark up shock
    Γ0[eq[:eq_phlps], endo[:λ_f_t]] = -(1 + m[:ι_p]*m[:β]*exp((1 - m[:σ_c])*m[:z_star]))/
        (1 + m[:ι_p]*m[:β]*exp((1 - m[:σ_c])*m[:z_star]))

    # Flexible prices and wages not necessary

    ### 10. Rental Rate of Capital

    # Sticky prices and wages
    Γ0[eq[:eq_caprnt], endo[:rk_t]] = 1.
    Γ0[eq[:eq_caprnt], endo[:k_t]]  = 1.
    Γ0[eq[:eq_caprnt], endo[:L_t]]  = -1.
    Γ0[eq[:eq_caprnt], endo[:w_t]]  = -1.

    # Flexible prices and wages
    Γ0[eq[:eq_caprnt_f], endo[:rk_f_t]] = 1.
    Γ0[eq[:eq_caprnt_f], endo[:k_f_t]]  = 1.
    Γ0[eq[:eq_caprnt_f], endo[:L_f_t]]  = -1.
    Γ0[eq[:eq_caprnt_f], endo[:w_f_t]]  = -1.

    ### 11. Marginal Substitution

    # Sticky prices and wages
    Γ0[eq[:eq_msub], endo[:μ_ω_t]] = 1.
    Γ0[eq[:eq_msub], endo[:L_t]]   = m[:ν_l]
    Γ0[eq[:eq_msub], endo[:c_t]]   = 1/(1 - m[:h]*exp(-m[:z_star]))
    Γ1[eq[:eq_msub], endo[:c_t]]   = m[:h]*exp(-m[:z_star])/(1 - m[:h]*exp(-m[:z_star]))
    Γ0[eq[:eq_msub], endo[:z_t]]   = m[:h]*exp(-m[:z_star]) /(1 - m[:h]*exp(-m[:z_star]))
    Γ0[eq[:eq_msub], endo[:w_t]]   = -1.

    # Flexible prices and wages
    Γ0[eq[:eq_msub_f], endo[:w_f_t]] = -1.
    Γ0[eq[:eq_msub_f], endo[:L_f_t]] = m[:ν_l]
    Γ0[eq[:eq_msub_f], endo[:c_f_t]] = 1/(1 - m[:h]*exp(-m[:z_star]))
    Γ1[eq[:eq_msub_f], endo[:c_f_t]] = m[:h]*exp(-m[:z_star])/(1 - m[:h]*exp(-m[:z_star]))
    Γ0[eq[:eq_msub_f], endo[:z_t]]   = m[:h]*exp(-m[:z_star])/(1 - m[:h]*exp(-m[:z_star]))

    ### 12. Evolution of Wages

    # Sticky prices and wages
    Γ0[eq[:eq_wage], endo[:w_t]]   = 1
    Γ0[eq[:eq_wage], endo[:μ_ω_t]] = (1 - m[:ζ_w]*m[:β]*exp((1 - m[:σ_c])*m[:z_star]))*
        (1 - m[:ζ_w])/(m[:ζ_w]*((m[:λ_w] - 1)*m[:ϵ_w] + 1))/(1 + m[:β]*exp((1 - m[:σ_c])*m[:z_star]))
    Γ0[eq[:eq_wage], endo[:π_t]]   = (1 + m[:ι_w]*m[:β]*exp((1 - m[:σ_c])*m[:z_star]))/(1 + m[:β]*exp((1 - m[:σ_c])*m[:z_star]))
    Γ1[eq[:eq_wage], endo[:w_t]]   = 1/(1 + m[:β]*exp((1 - m[:σ_c])*m[:z_star]))
    Γ0[eq[:eq_wage], endo[:z_t]]   = 1/(1 + m[:β]*exp((1 - m[:σ_c])*m[:z_star]))
    Γ1[eq[:eq_wage], endo[:π_t]]   = m[:ι_w]/(1 + m[:β]*exp((1 - m[:σ_c])*m[:z_star]))
    Γ0[eq[:eq_wage], endo[:Ew_t]]  = -m[:β]*exp((1 - m[:σ_c])*m[:z_star])/(1 + m[:β]*exp((1 - m[:σ_c])*m[:z_star]))
    Γ0[eq[:eq_wage], endo[:Ez_t]]  = -m[:β]*exp((1 - m[:σ_c])*m[:z_star])/(1 + m[:β]*exp((1 - m[:σ_c])*m[:z_star]))
    Γ0[eq[:eq_wage], endo[:Eπ_t]]  = -m[:β]*exp((1 - m[:σ_c])*m[:z_star])/(1 + m[:β]*exp((1 - m[:σ_c])*m[:z_star]))
    Γ0[eq[:eq_wage], endo[:λ_w_t]] = -1.

    # Flexible prices and wages not necessary

    ### 13. Monetary Policy Rule

    # Sticky prices and wages
    Γ0[eq[:eq_mp], endo[:R_t]]      = 1.
    Γ1[eq[:eq_mp], endo[:R_t]]      = m[:ρ]
    Γ0[eq[:eq_mp], endo[:π_t]]      = -(1 - m[:ρ])*m[:ψ1]
    Γ0[eq[:eq_mp], endo[:π_star_t]] = (1 - m[:ρ])*m[:ψ1]
    Γ0[eq[:eq_mp], endo[:y_t]]      = -(1 - m[:ρ])*m[:ψ2] - m[:ψ3]
    Γ0[eq[:eq_mp], endo[:y_f_t]]    = (1 - m[:ρ])*m[:ψ2] + m[:ψ3]
    Γ1[eq[:eq_mp], endo[:y_t]]      = -m[:ψ3]
    Γ1[eq[:eq_mp], endo[:y_f_t]]    = m[:ψ3]
    Γ0[eq[:eq_mp], endo[:rm_t]]     = -1.

    # Flexible prices and wages not necessary

    ### 14. Resource Constraint

    # Sticky prices and wages
    Γ0[eq[:eq_res], endo[:y_t]] = 1.
    Γ0[eq[:eq_res], endo[:g_t]] = -m[:g_star]
    Γ0[eq[:eq_res], endo[:c_t]] = -m[:cstar]/m[:ystar]
    Γ0[eq[:eq_res], endo[:i_t]] = -m[:istar]/m[:ystar]
    Γ0[eq[:eq_res], endo[:u_t]] = -m[:r_k_star]*m[:kstar]/m[:ystar]

    # Flexible prices and wages
    Γ0[eq[:eq_res_f], endo[:y_f_t]] = 1.
    Γ0[eq[:eq_res_f], endo[:g_t]]   = -m[:g_star]
    Γ0[eq[:eq_res_f], endo[:c_f_t]] = -m[:cstar]/m[:ystar]
    Γ0[eq[:eq_res_f], endo[:i_f_t]] = -m[:istar]/m[:ystar]
    Γ0[eq[:eq_res_f], endo[:u_f_t]] = -m[:r_k_star]*m[:kstar]/m[:ystar]

    ### 15. Extra States
    # These aren't strictly necessary, but they track lags or simplify the equations

    # π_t1
    Γ0[eq[:eq_π1], endo[:π_t1]] = 1.
    Γ1[eq[:eq_π1], endo[:π_t]]  = 1.

    # π_t2
    Γ0[eq[:eq_π2], endo[:π_t2]] = 1.
    Γ1[eq[:eq_π2], endo[:π_t1]] = 1.

    # π_a
    Γ0[eq[:eq_π_a], endo[:π_a_t]] = 1.
    Γ0[eq[:eq_π_a], endo[:π_t]]   = -1.
    Γ0[eq[:eq_π_a], endo[:π_t1]]  = -1.
    Γ0[eq[:eq_π_a], endo[:π_t2]]  = -1.
    Γ1[eq[:eq_π_a], endo[:π_t2]]  = 1.

    # Rt1
    Γ0[eq[:eq_Rt1], endo[:R_t1]] = 1.
    Γ1[eq[:eq_Rt1], endo[:R_t]]  = 1.

    # Ez_t
    Γ0[eq[:eq_Ez], endo[:Ez_t]]   = 1.
    Γ0[eq[:eq_Ez], endo[:ztil_t]] = -(m[:ρ_z]-1)/(1-m[:α])
    Γ0[eq[:eq_Ez], endo[:zp_t]]   = -m[:ρ_z_p]

    # Nominal Deposit rate
    Γ0[eq[:eq_dep], endo[:Rd_t]]   = 1.
    Γ0[eq[:eq_dep], endo[:R_t]]  = -1.
    Γ0[eq[:eq_dep], endo[:b_t]]   = (m[:σ_c]*(1 + m[:h]*exp(-m[:z_star])))/(1 - m[:h]*exp(-m[:z_star]))

    ### EXOGENOUS SHOCKS ###

    # Neutral technology
    Γ0[eq[:eq_z], endo[:z_t]]    = 1.
    Γ1[eq[:eq_z], endo[:ztil_t]] = (m[:ρ_z] - 1)/(1 - m[:α])
    Γ0[eq[:eq_z], endo[:zp_t]]   = -1.
    Ψ[eq[:eq_z], exo[:z_sh]]     = 1/(1 - m[:α])

    Γ0[eq[:eq_ztil], endo[:ztil_t]] = 1.
    Γ1[eq[:eq_ztil], endo[:ztil_t]] = m[:ρ_z]
    Ψ[eq[:eq_ztil], exo[:z_sh]]     = 1.

    # Long-run changes to productivity
    Γ0[eq[:eq_zp], endo[:zp_t]] = 1.
    Γ1[eq[:eq_zp], endo[:zp_t]] = m[:ρ_z_p]
    Ψ[eq[:eq_zp], exo[:zp_sh]]  = 1.

    # Government spending
    Γ0[eq[:eq_g], endo[:g_t]] = 1.
    Γ1[eq[:eq_g], endo[:g_t]] = m[:ρ_g]
    Ψ[eq[:eq_g], exo[:g_sh]]  = 1.
    Ψ[eq[:eq_g], exo[:z_sh]]  = m[:η_gz]

    # Asset shock
    Γ0[eq[:eq_b_til], endo[:b_til_t]] = 1.  # transient component
    Γ1[eq[:eq_b_til], endo[:b_til_t]] = m[:ρ_b_til]
    Ψ[eq[:eq_b_til], exo[:b_til_sh]]   = 1.

    Γ0[eq[:eq_b_p], endo[:b_p_t]] = 1.  # permanent component
    Γ1[eq[:eq_b_p], endo[:b_p_t]] = m[:ρ_b_p]
    Ψ[eq[:eq_b_p], exo[:b_p_sh]]   = 1.

    Γ0[eq[:eq_b], endo[:b_t]] = 1.
    Γ0[eq[:eq_b], endo[:b_til_t]] = -1.
    Γ0[eq[:eq_b], endo[:b_p_t]]   = -1.

    # Investment-specific technology
    Γ0[eq[:eq_μ], endo[:μ_t]] = 1.
    Γ1[eq[:eq_μ], endo[:μ_t]] = m[:ρ_μ]
    Ψ[eq[:eq_μ], exo[:μ_sh]]  = 1.

    # Price mark-up shock
    Γ0[eq[:eq_λ_f], endo[:λ_f_t]]  = 1.
    Γ1[eq[:eq_λ_f], endo[:λ_f_t]]  = m[:ρ_λ_f]
    Γ1[eq[:eq_λ_f], endo[:λ_f_t1]] = -m[:η_λ_f]
    Ψ[eq[:eq_λ_f], exo[:λ_f_sh]]   = 1.

    Γ0[eq[:eq_λ_f1], endo[:λ_f_t1]] = 1.
    Ψ[eq[:eq_λ_f1], exo[:λ_f_sh]]   = 1.

    # Wage mark-up shock
    Γ0[eq[:eq_λ_w], endo[:λ_w_t]]  = 1.
    Γ1[eq[:eq_λ_w], endo[:λ_w_t]]  = m[:ρ_λ_w]
    Γ1[eq[:eq_λ_w], endo[:λ_w_t1]] = -m[:η_λ_w]
    Ψ[eq[:eq_λ_w], exo[:λ_w_sh]]   = 1.

    Γ0[eq[:eq_λ_w1], endo[:λ_w_t1]] = 1.
    Ψ[eq[:eq_λ_w1], exo[:λ_w_sh]]   = 1.

    # Monetary policy shock
    Γ0[eq[:eq_rm], endo[:rm_t]] = 1.
    Γ1[eq[:eq_rm], endo[:rm_t]] = m[:ρ_rm]
    Ψ[eq[:eq_rm], exo[:rm_sh]]  = 1.

    ### Financial frictions

    # Standard deviation of capital shock to entrepreneurs
    Γ0[eq[:eq_σ_ω], endo[:σ_ω_t]] = 1.
    Γ1[eq[:eq_σ_ω], endo[:σ_ω_t]] = m[:ρ_σ_w]
    Ψ[eq[:eq_σ_ω], exo[:σ_ω_sh]]  = 1.

    # Exogenous bankruptcy costs
    Γ0[eq[:eq_μ_e], endo[:μ_e_t]] = 1.
    Γ1[eq[:eq_μ_e], endo[:μ_e_t]] = m[:ρ_μ_e]
    Ψ[eq[:eq_μ_e], exo[:μ_e_sh]]  = 1.

    # Fraction of entrepreneurs surviving period t
    Γ0[eq[:eq_γ], endo[:γ_t]] = 1.
    Γ1[eq[:eq_γ], endo[:γ_t]] = m[:ρ_γ]
    Ψ[eq[:eq_γ], exo[:γ_sh]]  = 1.

    # Long-term inflation expectations
    Γ0[eq[:eq_π_star], endo[:π_star_t]] = 1.
    Γ1[eq[:eq_π_star], endo[:π_star_t]] = m[:ρ_π_star]
    Ψ[eq[:eq_π_star], exo[:π_star_sh]]  = 1.

    # Anticipated policy shocks
    if n_anticipated_shocks(m) > 0

        # This section adds the anticipated shocks. There is one state for all the
        # anticipated shocks that will hit in a given period (i.e. rm_tl2 holds those that
        # will hit in two periods), and the equations are set up so that rm_tl2 last period
        # will feed into rm_tl1 this period (and so on for other numbers), and last period's
        # rm_tl1 will feed into the rm_t process (and affect the Taylor Rule this period).

        Γ1[eq[:eq_rm], endo[:rm_tl1]]   = 1.
        Γ0[eq[:eq_rml1], endo[:rm_tl1]] = 1.
        Ψ[eq[:eq_rml1], exo[:rm_shl1]]  = 1.

        if n_anticipated_shocks(m) > 1
            for i = 2:n_anticipated_shocks(m)
                Γ1[eq[Symbol("eq_rml$(i-1)")], endo[Symbol("rm_tl$i")]] = 1.
                Γ0[eq[Symbol("eq_rml$i")], endo[Symbol("rm_tl$i")]]     = 1.
                Ψ[eq[Symbol("eq_rml$i")], exo[Symbol("rm_shl$i")]]      = 1.
            end
        end
    end

    ### EXPECTATION ERRORS ###

    ### E(c)

    # Sticky prices and wages
    Γ0[eq[:eq_Ec], endo[:c_t]]  = 1.
    Γ1[eq[:eq_Ec], endo[:Ec_t]] = 1.
    Π[eq[:eq_Ec], ex[:Ec_sh]]   = 1.

    # Flexible prices and wages
    Γ0[eq[:eq_Ec_f], endo[:c_f_t]]  = 1.
    Γ1[eq[:eq_Ec_f], endo[:Ec_f_t]] = 1.
    Π[eq[:eq_Ec_f], ex[:Ec_f_sh]]   = 1.

    ### E(q)

    # Sticky prices and wages
    Γ0[eq[:eq_Eqk], endo[:qk_t]]  = 1.
    Γ1[eq[:eq_Eqk], endo[:Eqk_t]] = 1.
    Π[eq[:eq_Eqk], ex[:Eqk_sh]]   = 1.

    # Flexible prices and wages
    Γ0[eq[:eq_Eqk_f], endo[:qk_f_t]]  = 1.
    Γ1[eq[:eq_Eqk_f], endo[:Eqk_f_t]] = 1.
    Π[eq[:eq_Eqk_f], ex[:Eqk_f_sh]]   = 1.

    ### E(i)

    # Sticky prices and wages
    Γ0[eq[:eq_Ei], endo[:i_t]]  = 1.
    Γ1[eq[:eq_Ei], endo[:Ei_t]] = 1.
    Π[eq[:eq_Ei], ex[:Ei_sh]]   = 1.

    # Flexible prices and wages
    Γ0[eq[:eq_Ei_f], endo[:i_f_t]]  = 1.
    Γ1[eq[:eq_Ei_f], endo[:Ei_f_t]] = 1.
    Π[eq[:eq_Ei_f], ex[:Ei_f_sh]]   = 1.

    ### E(π)

    # Sticky prices and wages
    Γ0[eq[:eq_Eπ], endo[:π_t]]  = 1.
    Γ1[eq[:eq_Eπ], endo[:Eπ_t]] = 1.
    Π[eq[:eq_Eπ], ex[:Eπ_sh]]   = 1.

    ### E(l)

    # Sticky prices and wages
    Γ0[eq[:eq_EL], endo[:L_t]]  = 1.
    Γ1[eq[:eq_EL], endo[:EL_t]] = 1.
    Π[eq[:eq_EL], ex[:EL_sh]]   = 1.

    # Flexible prices and wages
    Γ0[eq[:eq_EL_f], endo[:L_f_t]]  = 1.
    Γ1[eq[:eq_EL_f], endo[:EL_f_t]] = 1.
    Π[eq[:eq_EL_f], ex[:EL_f_sh]]   = 1.

    ### E(rk)

    # Sticky prices and wages
    Γ0[eq[:eq_Erk], endo[:rk_t]]  = 1.
    Γ1[eq[:eq_Erk], endo[:Erk_t]] = 1.
    Π[eq[:eq_Erk], ex[:Erk_sh]]   = 1.

    # Flexible prices and wages
    Γ0[eq[:eq_Erktil_f], endo[:rktil_f_t]]  = 1.
    Γ1[eq[:eq_Erktil_f], endo[:ERktil_f_t]] = 1.
    Π[eq[:eq_Erktil_f], ex[:Erktil_f_sh]]   = 1.

    ### E(w)

    # Sticky prices and wages
    Γ0[eq[:eq_Ew], endo[:w_t]]  = 1.
    Γ1[eq[:eq_Ew], endo[:Ew_t]] = 1.
    Π[eq[:eq_Ew], ex[:Ew_sh]]   = 1.

    ### E(Rktil)

    # Sticky prices and wages
    Γ0[eq[:eq_ERktil], endo[:Rktil_t]]  = 1.
    Γ1[eq[:eq_ERktil], endo[:ERtil_k_t]] = 1.
    Π[eq[:eq_ERktil], ex[:ERktil_sh]]    = 1.

    return Γ0, Γ1, C, Ψ, Π
end


text:
using Luna

radius = 13e-6 # HCF core radius
flength = 2 # HCF length

gas = :H2
pressure = 5.0 # gas pressure in bar

λ0 = 800e-9 # central wavelength of the pump pulse
τfwhm = 20e-15 # FWHM duration of the pump pulse
energy = 1e-6 # energy in the pump pulse

rssfs = prop_capillary(radius, flength, gas, pressure; λ0, τfwhm, energy,
                     trange=40e-12, λlims=(400e-9, 2e-6),
                     raman=true, envelope=true, loss=false, plasma=false)

Plotting.prop_2D(rssfs, :λ; modes=:sum, trange=(-0.5e-12, 2e-12), λrange=(500e-9, 1400e-9),
                 dBmin=-30)
Plotting.time_1D(rssfs, range(0,flength,length=5), modes=:sum, trange=(-0.5e-12, 2e-12))
Plotting.spec_1D(rssfs, range(0,flength,length=5), modes=:sum, λrange=(500e-9, 1400e-9))

ext: jl


text:
struct minusone end
import Base:^
^(x::minusone, n::Real) = isodd(n) ? -1 : 1
macro x_str(s::String)
    minusone()
end

max_stars_repo_name: mmikhasenko/ThreeBodyDecay

max_issues_repo_name: mmikhasenko/ThreeBodyDecay

max_forks_repo_name: mmikhasenko/ThreeBodyDecay


text:
@testset "triangle" begin
    r = 2  # radius
    b = Polygon{3}([0,0], r)  # regular triangle inscribed in circle of radius r (triangel side = √3r/2)

    va = [0,r]  # apex
    vl = [-√3r/2,-r/2]  # left vertex
    vr = [√3r/2,-r/2]  # right vertex
    vs = (va, vl, vr)

    hb = [0,-r/2]  # midpoint of base
    hr = [√3r/4,r/4]  # midpoint of right side
    hl = [-√3r/4,r/4]  # midpoint of left side

    @test b == deepcopy(b)
    @test hash(b) == hash(deepcopy(b))
    @test [0,0] ∈ b  # center
    @test (s = one⁻; s.*va ∈ b && s.*vl ∈ b && s.*vr ∈ b)  # inside apex, left, right
    @test (s = one⁺; s.*va ∉ b && s.*vl ∉ b && s.*vr ∉ b)  # outside apex, left, right

    @test (xn = surfpt_nearby([0,0],b); xn≈(hb,normalize(hb)) || xn≈(hr,normalize(hr)) || xn≈(hl,normalize(hl)))  # handle point at center properly

    @test all([(p₀ = hb; p = s.*p₀; surfpt_nearby(p,b) ≈ (p₀, normalize(p₀))) for s = (one⁻⁻,1,one⁺⁺)])  # around base
    @test all([(p₀ = hr; p = s.*p₀; surfpt_nearby(p,b) ≈ (p₀, normalize(p₀))) for s = (one⁻⁻,1,one⁺⁺)])  # around right side
    @test all([(p₀ = hl; p = s.*p₀; surfpt_nearby(p,b) ≈ (p₀, normalize(p₀))) for s = (one⁻⁻,1,one⁺⁺)])  # around left side

    @test all([(p = one⁺⁺.*p₀; surfpt_nearby(p,b) ≈ (p₀, normalize(p₀))) for p₀ = vs])  # outside corners
    @test all([(p = one⁺⁺.*[p₀[1]*0.98, p₀[2]*1.02]; surfpt_nearby(p,b) ≈ (p₀, normalize(p-p₀))) for p₀ = vs])  # outside corners
    @test all([(p = p₀; surfpt_nearby(p,b) ≈ (p₀, normalize(p₀))) for p₀ = vs])  # at corners

    @test bounds(b) ≈ (vl, [vr[1],va[2]])
    @test checkbounds(b)

    @test (∆ = rand(2); translate(b,∆) ≈ Polygon{3}(∆,r,0.0))

    @test b ≈ Polygon([va vl vr])
    @test Polygon([va vl vr]) == Polygon([va vr vl]) == Polygon([vl va vr]) == Polygon([vl vr va]) == Polygon([vr va vl]) == Polygon([vr vl va])
    @test b ≈ Isosceles((vl,vr), 1.5r)
    @test b ≈ Isosceles([0,-r/2], √3r, 1.5r)
end  # @testset "triangle"

@testset "trapezoid" begin
    b2 = 1
    h = 1
    ϕ = π/3
    t = Trapezoid(([-b2,0],[b2,0]), h, ϕ)

    @test t ≈ Trapezoid(([-b2,0],[b2,0]), h, (ϕ,ϕ))
    @test t ≈ Trapezoid([0,0], 2b2, h, ϕ)
end  # @test "trapezoid"


text:
export geam!, gemm_batched!, gemm_batched

for (f,T) in (
    (:(:cublasSgeam),:Float32),
    (:(:cublasDgeam),:Float64))
    @eval begin
        function geam!(tA::Char, tB::Char, alpha::$T, A::CuMatrix{$T},
            beta::$T, B::CuMatrix{$T}, C::CuMatrix{$T})

            m = size(A, tA == 'N' ? 1 : 2)
            n = size(A, tA == 'N' ? 2 : 1)
            if m != size(C,1) || n != size(C,2)
                throw(DimensionMismatch(""))
            end

            @cublas($f, (
                Ptr{Cvoid},Cint,Cint,Cint,Cint,
                Ptr{$T},Ptr{$T},Cint,
                Ptr{$T},Ptr{$T},Cint,
                Ptr{$T},Cint),
                gethandle(), cublasop(tA), cublasop(tB), m, n,
                [alpha], A, stride(A,2),
                [beta], B, stride(B,2),
                C, stride(C,2))
            C
        end
    end
end

function Base.transpose(x::CuMatrix{T}) where T
    t = similar(x, size(x,2), size(x,1))
    geam!('T', 'N', T(1), x, T(0), t, t)
    t
end
Base.transpose(x::CuVector) = transpose(reshape(x,length(x),1))

for (f,T) in (
        (:(:cublasDgemmBatched),:Float64),
        (:(:cublasSgemmBatched),:Float32))
    @eval begin
        function gemm_batched!(tA::Char, tB::Char, m::Int, n::Int, k::Int,
            alpha::$T, As::Vector{Ptr{$T}}, lda::Int, Bs::Vector{Ptr{$T}}, ldb::Int,
            beta::$T, Cs::Vector{Ptr{$T}}, ldc::Int, batchcount::Int)

            if (length(As) != length(Bs) || length(As) != length(Cs))
                throw(DimensionMismatch(""))
            end

            @cublas($f, (
                Ptr{Cvoid},Cint,Cint,Cint,Cint,Cint,
                Ptr{$T},Ptr{Ptr{$T}},Cint,
                Ptr{Ptr{$T}},Cint,
                Ptr{$T},Ptr{Ptr{$T}},Cint,Cint),
                gethandle(), cublasop(tA), cublasop(tB), m, n, k,
                [alpha], CuArray(As), lda, CuArray(Bs), ldb,
                [beta], CuArray(Cs), ldc, batchcount)
            Cs
        end
    end
end

function gemm_batched!(tA::Char, tB::Char,
    alpha::T, As::Vector{CuMatrix{T}}, Bs::Vector{CuMatrix{T}},
    beta::T, Cs::Vector{CuMatrix{T}}) where T

    for (A,B,C) in zip(As,Bs,Cs)
        m = size(A, tA == 'N' ? 1 : 2)
        kA = size(A, tA == 'N' ? 2 : 1)
        kB = size(B, tB == 'N' ? 1 : 2)
        n = size(B, tB == 'N' ? 2 : 1)
        if m != size(C,1) || n != size(C,2) || kA != kB
            throw(DimensionMismatch(""))
        end
    end
    m = size(As[1], tA == 'N' ? 1 : 2)
    k = size(As[1], tA == 'N' ? 2 : 1)
    n = size(Bs[1], tB == 'N' ? 2 : 1)
    ptrAs = map(pointer, As)
    ptrBs = map(pointer, Bs)
    ptrCs = map(pointer, Cs)
    lda = stride(As[1], 2)
    ldb = stride(Bs[1], 2)
    ldc = stride(Cs[1], 2)
    gemm_batched!(tA, tB, m, n, k, alpha, ptrAs, lda, ptrBs, ldb, beta, ptrCs, ldc, length(As))
    Cs
end

function gemm_batched!(tA::Char, tB::Char, alpha::T, A::CuArray{T,3}, B::CuArray{T,3}, beta::T, C::CuArray{T,3}) where T
    @assert size(A,3) == size(B,3) == size(C,3)
    m = size(A, tA == 'N' ? 1 : 2)
    kA = size(A, tA == 'N' ? 2 : 1)
    kB = size(B, tB == 'N' ? 1 : 2)
    n = size(B, tB == 'N' ? 2 : 1)
    if m != size(C,1) || n != size(C,2) || kA != kB
        throw(DimensionMismatch(""))
    end
    k = kA
    ptrAs = [pointer(A,(i-1)*size(A,1)*size(A,2)+1) for i=1:size(A,3)]
    ptrBs = [pointer(B,(i-1)*size(B,1)*size(B,2)+1) for i=1:size(B,3)]
    ptrCs = [pointer(C,(i-1)*size(C,1)*size(C,2)+1) for i=1:size(C,3)]
    lda = stride(A, 2)
    ldb = stride(B, 2)
    ldc = stride(C, 2)
    gemm_batched!(tA, tB, m, n, k, alpha, ptrAs, lda, ptrBs, ldb, beta, ptrCs, ldc, size(A,3))
    C
end

function gemm_batched(tA::Char, tB::Char, alpha::T, As::Vector{CuMatrix{T}}, Bs::Vector{CuMatrix{T}}) where T
    m = size(As[1], tA == 'N' ? 1 : 2)
    n = size(Bs[1], tB == 'N' ? 2 : 1)
    Cs = [CuArray{T}(m,n) for i=1:length(As)]
    gemm_batched!(tA, tB, alpha, As, Bs, T(0), Cs)
end

function gemm_batched(tA::Char, tB::Char, alpha::T, A::CuArray{T,3}, B::CuArray{T,3}) where T
    m = size(A, tA == 'N' ? 1 : 2)
    n = size(B, tB == 'N' ? 2 : 1)
    C = CuArray{T}(m, n, size(A,3))
    gemm_batched!(tA, tB, alpha, A, B, T(0), C)
end


text:
while true
       start = 0
       space = 110
       x = fill('.', space)
       i = 0

       d = Dict{String,Char}()
       open("day12.txt") do f
       for (k,v) in enumerate(split(readline(f), ""))
       x[1, k + 5] = v[1]
       end
       for line in eachline(f)
       d[line[1:5]] = line[7]
       end
       end
       y = deepcopy(x)

       while i < 50000000000
       x = deepcopy(y)
       y = fill('.', space)
       for j in 1:146
       b = string(x[j:j + 4]...)
       y[j + 2] = (b in keys(d) ? d[b] : '.')
       end

       if any(==('#'), y[1:4])
       space += 2; start -= 2
       else
       space -= 2; start += 2
       any(==('#'), y[end - 4:end]) && (space += 2)

       end; s = 0
       for i in 1:150
       x[21,i] == '#' && (s += i - 26)
       end; s
       i += 1; i > 50000000000 && break

max_stars_repo_name: AnAverageHuman/competitive

max_issues_repo_name: AnAverageHuman/competitive

max_forks_repo_name: AnAverageHuman/competitive


text:
using Tabben

using Test

include("utils.jl")

@testset "unit tests" begin

    include("test_dataset_functions.jl")
    include("test_dataset_loading.jl")

end;


text:
"""
    Sequential quadratic programming with trust region
"""
abstract type AbstractSqpTrOptimizer <: AbstractSqpOptimizer end

mutable struct SqpTR{T,TD,TI} <: AbstractSqpTrOptimizer
    @sqp_fields

    # directions for multipliers
    p_lambda::TD
    p_mult_x_L::TD
    p_mult_x_U::TD

    E_soc::TD # constraint evaluation for SOC
    soc::TD # second-order correction direction

    phi::T # merit function value
    μ::T # penalty parameter

    Δ::T # current trust region size
    Δ_min::T # minimum trust region size allowed
    Δ_max::T # maximum trust region size allowed

    step_acceptance::Bool

    function SqpTR(problem::Model{T,TD}, TI = Vector{Int}) where {T,TD<:AbstractArray{T}}
        sqp = new{T,TD,TI}()
        sqp.problem = problem
        sqp.x = deepcopy(problem.x)
        sqp.p = zeros(T, problem.n)
        sqp.p_soc = zeros(T, problem.n)
        sqp.p_slack = Dict()
        sqp.lambda = zeros(T, problem.m)
        sqp.mult_x_L = zeros(T, problem.n)
        sqp.mult_x_U = zeros(T, problem.n)
        sqp.df = TD(undef, problem.n)
        sqp.E = TD(undef, problem.m)
        sqp.dE = TD(undef, length(problem.j_str))

        # FIXME: Replace Vector{Int} with TI?
        sqp.j_row = TI(undef, length(problem.j_str))
        sqp.j_col = TI(undef, length(problem.j_str))
        for i = 1:length(problem.j_str)
            sqp.j_row[i] = Int(problem.j_str[i][1])
            sqp.j_col[i] = Int(problem.j_str[i][2])
        end
        sqp.Jacobian =
            sparse(sqp.j_row, sqp.j_col, ones(length(sqp.j_row)), problem.m, problem.n)
        sqp.h_row = TI(undef, length(problem.h_str))
        sqp.h_col = TI(undef, length(problem.h_str))
        for i = 1:length(problem.h_str)
            sqp.h_row[i] = Int(problem.h_str[i][1])
            sqp.h_col[i] = Int(problem.h_str[i][2])
        end
        sqp.h_val = TD(undef, length(problem.h_str))
        sqp.Hessian =
            sparse(sqp.h_row, sqp.h_col, ones(length(sqp.h_row)), problem.n, problem.n)

        sqp.p_lambda = zeros(T, problem.m)
        sqp.p_mult_x_L = zeros(T, problem.n)
        sqp.p_mult_x_U = zeros(T, problem.n)

        sqp.E_soc = TD(undef, problem.m)
        sqp.soc = zeros(T, problem.n)

        sqp.phi = 1.0e+20
        sqp.μ = 1.0e+4
        sqp.Δ = 10.0
        sqp.Δ_min = 1.0e-4
        sqp.Δ_max = 1.0e+4
        sqp.step_acceptance = true

        sqp.prim_infeas = Inf
        sqp.dual_infeas = Inf
        sqp.compl = Inf

        sqp.options = problem.parameters
        sqp.optimizer = nothing
        sqp.sub_status = nothing

        sqp.feasibility_restoration = false
        sqp.iter = 1
        sqp.ret = -5
        sqp.start_time = 0.0
        sqp.start_iter_time = 0.0

        sqp.tmpx = TD(undef, problem.n)
        sqp.tmpE = TD(undef, problem.m)

        return sqp
    end
end

"""
    run!

Run the line-search SQP algorithm
"""
function run!(sqp::AbstractSqpTrOptimizer)

    sqp.μ = sqp.options.init_mu
    sqp.Δ = sqp.options.tr_size

    sqp.start_time = time()

    if sqp.options.OutputFlag == 0
        Logging.disable_logging(Logging.Info)
    end

    print_header(sqp)

    # Find the initial point feasible to linear and bound constraints
    lpviol = violation_of_linear_constraints(sqp, sqp.x)
    if lpviol > sqp.options.tol_infeas
        @info "Initial point not feasible to linear constraints..."
        sub_optimize_lp!(sqp)

        print(sqp, "LP")
    else
        @info "Initial point feasible to linear constraints..." lpviol
    end

    while true

        # Iteration counter limit
        if terminate_by_iterlimit(sqp)
            break
        end

        sqp.start_iter_time = time()

        # evaluate function, constraints, gradient, Jacobian
        if sqp.step_acceptance
            eval_functions!(sqp)
            sqp.prim_infeas = norm_violations(sqp, 1)
            sqp.dual_infeas = KT_residuals(sqp)
            sqp.compl = norm_complementarity(sqp)
        end

        # solve QP subproblem
        QP_time = @elapsed compute_step!(sqp)
        add_statistics(sqp.problem, "QP_time", QP_time)

        if sqp.sub_status ∈ [MOI.OPTIMAL, MOI.ALMOST_OPTIMAL, MOI.ALMOST_LOCALLY_SOLVED, MOI.LOCALLY_SOLVED]
            # do nothing
        elseif sqp.sub_status ∈ [MOI.INFEASIBLE, MOI.LOCALLY_INFEASIBLE, MOI.DUAL_INFEASIBLE, MOI.NORM_LIMIT, MOI.OBJECTIVE_LIMIT]
            if sqp.feasibility_restoration == true
                @info "Failed to find a feasible direction"
                if sqp.prim_infeas <= sqp.options.tol_infeas
                    sqp.ret = 6
                else
                    sqp.ret = 2
                end
                break
            else
                @info "Feasibility restoration starts... (status: $(sqp.sub_status))"
                # println("Feasibility restoration ($(sqp.sub_status), |p| = $(norm(sqp.p, Inf))) begins.")
                sqp.feasibility_restoration = true
                print(sqp)
                collect_statistics(sqp)
                sqp.iter += 1
                continue
            end
        else
            sqp.ret == -3
            if sqp.prim_infeas <= sqp.options.tol_infeas
                sqp.ret = 6
            end
            break
        end

        if sqp.step_acceptance
            sqp.phi = compute_phi(sqp, sqp.x, 0.0, sqp.p)
        end

        print(sqp)
        collect_statistics(sqp)

        if norm(sqp.p, Inf) <= sqp.options.tol_direction
            if sqp.feasibility_restoration
                sqp.feasibility_restoration = false
                sqp.iter += 1
                continue
            else
                sqp.ret = 0
                break
            end
        end

        if sqp.prim_infeas <= sqp.options.tol_infeas &&
           sqp.compl <= sqp.options.tol_residual #&&
        #    norm(sqp.p, Inf) <= sqp.options.tol_direction
            if sqp.feasibility_restoration
                sqp.feasibility_restoration = false
                sqp.iter += 1
                continue
            elseif sqp.dual_infeas <= sqp.options.tol_residual
                sqp.ret = 0
                break
            end
        end

        do_step!(sqp)

        # NOTE: This is based on the algorithm of filterSQP.
        if sqp.feasibility_restoration && sqp.step_acceptance
            sqp.feasibility_restoration = false
        end

        sqp.iter += 1
    end
    sqp.problem.obj_val = sqp.problem.eval_f(sqp.x)
    sqp.problem.status = Int(sqp.ret)
    sqp.problem.x .= sqp.x
    sqp.problem.g .= sqp.E
    sqp.problem.mult_g .= sqp.lambda
    sqp.problem.mult_x_U .= sqp.mult_x_U
    sqp.problem.mult_x_L .= sqp.mult_x_L
    add_statistic(sqp.problem, "iter", sqp.iter)
end

"""
    violation_of_linear_constraints

Compute the violation of linear constraints at a given point `x`

# Arguments
- `sqp`: SQP model struct
- `x`: solution to evaluate the violations

# Note
This function assumes that the first `sqp.problem.num_linear_constraints` constraints are linear.
"""
function violation_of_linear_constraints(sqp::AbstractSqpTrOptimizer, x::TD)::T where {T, TD <: AbstractVector{T}}
    # evaluate constraints
    sqp.problem.eval_g(x, sqp.E)

    lpviol = 0.0
    for i = 1:sqp.problem.num_linear_constraints
        lpviol += max(0.0, sqp.problem.g_L[i] - sqp.E[i])
        lpviol -= min(0.0, sqp.problem.g_U[i] - sqp.E[i])
    end
    for i = 1:sqp.problem.n
        lpviol += max(0.0, sqp.problem.x_L[i] - x[i])
        lpviol -= min(0.0, sqp.problem.x_U[i] - x[i])
    end
    return lpviol
end

"""
    sub_optimize_lp!

Compute the initial point that is feasible to linear constraints and variable bounds.

# Arguments
- `sqp`: SQP model struct
"""
function sub_optimize_lp!(sqp::AbstractSqpTrOptimizer)
    sqp.f = sqp.problem.eval_f(sqp.x)
    sqp.problem.eval_grad_f(sqp.x, sqp.df)
    eval_Jacobian!(sqp)
    if 1 == 1
        sqp.x, sqp.lambda, sqp.mult_x_U, sqp.mult_x_L, sqp.sub_status = sub_optimize_lp(
            sqp.options.external_optimizer,
            sqp.Jacobian, sqp.problem.g_L, sqp.problem.g_U,
            sqp.problem.x_L, sqp.problem.x_U, sqp.x,
            sqp.problem.num_linear_constraints, sqp.problem.m
        )
    else
        fill!(sqp.E, 0.0)

        sqp.optimizer =
            SubOptimizer(
                JuMP.Model(sqp.options.external_optimizer),
                QpData(
                    MOI.MIN_SENSE,
                    nothing,
                    sqp.df,
                    sqp.Jacobian,
                    sqp.E,
                    sqp.problem.g_L,
                    sqp.problem.g_U,
                    sqp.problem.x_L,
                    sqp.problem.x_U,
                    sqp.problem.num_linear_constraints
                )
            )
        create_model!(sqp.optimizer, sqp.Δ)
        sqp.x, sqp.lambda, sqp.mult_x_U, sqp.mult_x_L, sqp.sub_status = sub_optimize_lp(sqp.optimizer, sqp.x)
    end

    # TODO: Do we need to discard small numbers?
    dropzeros!(sqp.x)
    dropzeros!(sqp.lambda)
    dropzeros!(sqp.mult_x_U)
    dropzeros!(sqp.mult_x_L)
    return
end

"""
    sub_optimize!

Solve trust-region QP subproblem. If in feasibility restoration phase, the feasibility restoration subproblem is solved.

# Arguments
- `sqp`: SQP model struct
"""
function sub_optimize!(sqp::AbstractSqpTrOptimizer)
    if isnothing(sqp.optimizer)
        sqp.optimizer = SubOptimizer(
            JuMP.Model(sqp.options.external_optimizer),
            QpData(sqp),
        )
        create_model!(sqp.optimizer, sqp.Δ)
    else
        sqp.optimizer.data = QpData(sqp)
    end
    # TODO: This can be modified to Sl1QP.
    if sqp.feasibility_restoration
        return sub_optimize_FR!(sqp.optimizer, sqp.x, sqp.Δ)
    else
        return sub_optimize!(sqp.optimizer, sqp.x, sqp.Δ)
        # return sub_optimize_L1QP!(sqp.optimizer, sqp.x, sqp.Δ, sqp.μ)
    end
end

"""
    sub_optimize_soc!

Solve second-order correction QP subproblem.

# Arguments
- `sqp`: SQP model struct
"""
function sub_optimize_soc!(sqp::AbstractSqpTrOptimizer)
    sqp.problem.eval_g(sqp.x + sqp.p, sqp.E_soc)
    sqp.E_soc -= sqp.Jacobian * sqp.p
    sqp.optimizer.data = QpData(
        MOI.MIN_SENSE,
        sqp.Hessian,
        sqp.df,
        sqp.Jacobian,
        sqp.E_soc,
        sqp.problem.g_L,
        sqp.problem.g_U,
        sqp.problem.x_L,
        sqp.problem.x_U,
        sqp.problem.num_linear_constraints
    )
    p, _, _, _, _, _ = sub_optimize!(sqp.optimizer, sqp.x, sqp.Δ)
    sqp.p_soc .= sqp.p .+ p
    return nothing
    # return sub_optimize_L1QP!(sqp.optimizer, sqp.x, sqp.Δ, sqp.μ)
end

"""
    compute_step!

Compute the step direction with respect to priaml and dual variables by solving QP subproblem and also updates the penalty parameter μ.

# Arguments
- `sqp`: SQP model struct
"""
function compute_step!(sqp::AbstractSqpTrOptimizer)

    @info "solve QP subproblem..."
    sqp.p, lambda, mult_x_U, mult_x_L, sqp.p_slack, sqp.sub_status = sub_optimize!(sqp)

    sqp.p_lambda = lambda - sqp.lambda
    sqp.p_mult_x_L = mult_x_L - sqp.mult_x_L
    sqp.p_mult_x_U = mult_x_U - sqp.mult_x_U
    sqp.μ = max(sqp.μ, norm(sqp.lambda, Inf))
    @info "...found a direction"
end

"""
    compute_step_Sl1QP!

Compute the step direction with respect to priaml and dual variables by solving an elastic-mode QP subproblem and also updates the penalty parameter μ.

# Arguments
- `sqp`: SQP model struct

# Note
This is not currently used.
"""
function compute_step_Sl1QP!(sqp::AbstractSqpTrOptimizer)

    ϵ_1 = 0.9
    ϵ_2 = 0.1

    sqp.p, lambda, mult_x_U, mult_x_L, sqp.p_slack, sqp.sub_status = sub_optimize!(sqp)

    if sqp.sub_status ∈ [MOI.OPTIMAL, MOI.ALMOST_LOCALLY_SOLVED, MOI.LOCALLY_SOLVED]
        # compute the constraint violation
        m_0 = norm_violations(sqp, 1)
        m_μ = 0.0
        for (_, slacks) in sqp.p_slack
            m_μ += sum(slacks)
        end

        if m_μ > 1.0e-8
            p, infeasibility = sub_optimize_infeas(sqp.optimizer, sqp.x, sqp.Δ)
            # @show m_μ, infeasibility
            if infeasibility < 1.0e-8
                while m_μ > 1.0e-8 && sqp.μ < sqp.options.max_mu
                    sqp.μ = min(10.0 * sqp.μ, sqp.options.max_mu)

                    sqp.p, lambda, mult_x_U, mult_x_L, sqp.p_slack, sqp.sub_status = sub_optimize_L1QP!(sqp.optimizer, sqp.x, sqp.Δ, sqp.μ)

                    m_μ = 0.0
                    for (_, slacks) in sqp.p_slack, s in slacks
                        m_μ += s
                    end
                    @info "L1QP solve for feasible QP" infeasibility sqp.μ sqp.sub_status m_μ
                end
            else
                m_inf = norm_violations(
                    sqp.E + sqp.Jacobian * p,
                    sqp.problem.g_L,
                    sqp.problem.g_U,
                    sqp.x + p,
                    sqp.problem.x_L,
                    sqp.problem.x_U,
                    1,
                )
                while m_0 - m_μ < ϵ_1 * (m_0 - m_inf) && sqp.μ < sqp.options.max_mu
                    sqp.μ = min(10.0 * sqp.μ, sqp.options.max_mu)
                    sqp.p, lambda, mult_x_U, mult_x_L, sqp.p_slack, sqp.sub_status = sub_optimize_L1QP!(sqp.optimizer, sqp.x, sqp.Δ, sqp.μ)

                    m_μ = 0.0
                    for (_, slacks) in sqp.p_slack, s in slacks
                        m_μ += s
                    end
                    @info "L1QP solve for infeasible QP" infeasibility sqp.μ m_0 m_μ
                end
            end
        end

        # q_0 = compute_qmodel(sqp, false)
        # q_k = compute_qmodel(sqp, true)
        # @info "L1QP solve for μ+" q_0 q_k m_0 m_μ
        # while q_0 - q_k < ϵ_2 * sqp.μ * (m_0 - m_μ)
        #     sqp.μ = min(2.0 * sqp.μ, sqp.options.max_mu)
        #     sqp.p, lambda, mult_x_U, mult_x_L, sqp.p_slack, sqp.sub_status = sub_optimize_L1QP!(sqp.optimizer, sqp.x, sqp.Δ, sqp.μ)

        #     m_μ = 0.0
        #     for (_, slacks) in sqp.p_slack, s in slacks
        #         m_μ += s
        #     end
        #     q_k = compute_qmodel(sqp, true)
        #     @info "L1QP solve for μ+" q_0 q_k m_0 m_μ
        # end

    else
        @error "Unexpected QP subproblem status $(sqp.sub_status)"
    end

    @info "...solved QP subproblem"

    sqp.p_lambda = lambda - sqp.lambda
    sqp.p_mult_x_L = mult_x_L - sqp.mult_x_L
    sqp.p_mult_x_U = mult_x_U - sqp.mult_x_U
    sqp.μ = max(sqp.μ, norm(sqp.lambda, Inf))
end

"""
    compute_qmodel

Evaluate the quadratic model q(p) with ℓ₁ penalty term, which is given by
q(p) = fₖ + ∇fₖᵀp + 0.5 pᵀ ∇ₓₓ²Lₖ p + μ ∑ᵢ|cᵢ(xₖ) + ∇cᵢ(xₖ)ᵀp| + μ ∑ᵢ[cᵢ(xₖ) + ∇cᵢ(xₖ)ᵀp]⁻

# Arguments
- `sqp::SqpTR`: SQP model struct
- `p::TD`: direction vector
- `with_step::Bool`: `true` for q(p); `false` for `q(0)`

# Note
For p=0, the model is simplified to q(0) = μ ∑ᵢ|cᵢ(xₖ)| + μ ∑ᵢ[cᵢ(xₖ)]⁻
"""
function compute_qmodel(sqp::AbstractSqpTrOptimizer, p::TD, with_step::Bool = false) where {T, TD<:AbstractArray{T}}
    qval = 0.0
    if with_step
        qval += sqp.df' * p + 0.5 * p' * sqp.Hessian * p
        sqp.tmpx .= sqp.x .+ p
        sqp.tmpE .= sqp.E .+ sqp.Jacobian * p
    else
        sqp.tmpx .= sqp.x
        sqp.tmpE .= sqp.E
    end
    qval += sqp.μ * norm_violations(
        sqp.tmpE,
        sqp.problem.g_L,
        sqp.problem.g_U,
        sqp.tmpx,
        sqp.problem.x_L,
        sqp.problem.x_U,
        1,
    )
    return qval
end
compute_qmodel(sqp::AbstractSqpTrOptimizer, with_step::Bool = false) = compute_qmodel(sqp, sqp.p, with_step)

"""
    do_step!

Test the step `p` whether to accept or reject.
"""
function do_step!(sqp::AbstractSqpTrOptimizer)

    ϕ_k = compute_phi(sqp, sqp.x, 1.0, sqp.p)
    ared = sqp.phi - ϕ_k
    # @show sqp.phi, ϕ_k

    pred = 1.0
    if !sqp.feasibility_restoration
        q_0 = compute_qmodel(sqp, false)
        q_k = compute_qmodel(sqp, true)
        pred = q_0 - q_k
        # @show q_0, q_k
    end

    ρ = ared / pred
    if ared > 0 && ρ > 0
        sqp.x .+= sqp.p
        sqp.lambda .+= sqp.p_lambda
        sqp.mult_x_L .+= sqp.p_mult_x_L
        sqp.mult_x_U .+= sqp.p_mult_x_U
        if sqp.Δ == norm(sqp.p, Inf)
            sqp.Δ = min(2 * sqp.Δ, sqp.Δ_max)
        end
        sqp.step_acceptance = true
    else
        sqp.tmpx .= sqp.x .+ sqp.p
        c_k = norm_violations(sqp, sqp.tmpx)

        perform_soc = false
        if sqp.options.use_soc
            if c_k > 0 && sqp.feasibility_restoration == false
                @info "Try second-order correction..."

                # sqp.p should be adjusted inside sub_optimize_soc!
                sub_optimize_soc!(sqp)

                ϕ_soc = compute_phi(sqp, sqp.x, 1.0, sqp.p_soc)
                ared = sqp.phi - ϕ_soc
                pred = 1.0
                if !sqp.feasibility_restoration
                    q_soc = compute_qmodel(sqp, sqp.p_soc, true)
                    pred = q_0 - q_soc
                end
                ρ_soc = ared / pred
                if ared > 0 && ρ_soc > 0
                    @info "SOC" ϕ_k ϕ_soc ared pred ρ_soc
                    @info "...second-order correction added"
                    sqp.x .+= sqp.p_soc
                    sqp.lambda .+= sqp.p_lambda
                    sqp.mult_x_L .+= sqp.p_mult_x_L
                    sqp.mult_x_U .+= sqp.p_mult_x_U
                    sqp.step_acceptance = true
                    perform_soc = true
                else
                    @info "...second-order correction discarded"
                end
            end
        end

        if !perform_soc
            sqp.Δ = 0.5 * min(sqp.Δ, norm(sqp.p, Inf))
            sqp.step_acceptance = false
        end
    end
end

function print_header(sqp::AbstractSqpTrOptimizer)
    if sqp.options.OutputFlag == 0
        return
    end
    @printf("  %6s", "iter")
    @printf(" ")
    @printf("  %15s", "f(x_k)")
    @printf("  %15s", "ϕ(x_k)")
    @printf("  %15s", "μ")
    @printf("  %15s", "|λ|∞")
    @printf("  %14s", "Δ")
    @printf("  %14s", "|p|")
    @printf("  %14s", "inf_pr")
    @printf("  %14s", "inf_du")
    @printf("  %14s", "compl")
    @printf("  %10s", "time")
    @printf("\n")
end

"""
    print

Print iteration information.
"""
function print(sqp::AbstractSqpTrOptimizer, status_mark = "  ")
    if sqp.options.OutputFlag == 0
        return
    end
    if sqp.iter > 1 && (sqp.iter - 1) % 25 == 0
        print_header(sqp)
    end
    st = ifelse(sqp.feasibility_restoration, "FR", status_mark)
    @printf("%2s%6d", st, sqp.iter)
    @printf("%1s", ifelse(sqp.step_acceptance, "a", "r"))
    @printf("  %+6.8e", sqp.f)
    @printf("  %+6.8e", sqp.phi)
    @printf("  %+6.8e", sqp.μ)
    @printf("  %+6.8e", norm(sqp.lambda,Inf))
    @printf("  %6.8e", sqp.Δ)
    @printf("  %6.8e", norm(sqp.p, Inf))
    if isinf(sqp.prim_infeas)
        @printf("  %14s", "Inf")
    else
        @printf("  %6.8e", sqp.prim_infeas)
    end
    if isinf(sqp.dual_infeas)
        @printf("  %14s", "Inf")
    else
        @printf("  %6.8e", sqp.dual_infeas)
    end
    if isinf(sqp.compl)
        @printf("  %14s", "Inf")
    else
        @printf("  %6.8e", sqp.compl)
    end
    @printf("  %10.2f", time() - sqp.start_time)
    @printf("\n")
end

"""
    collect_statistics

Collect iteration information.
"""
function collect_statistics(sqp::AbstractSqpTrOptimizer)
    if sqp.options.StatisticsFlag == 0
        return
    end
    add_statistics(sqp.problem, "f(x)", sqp.f)
    add_statistics(sqp.problem, "ϕ(x_k))", sqp.phi)
    add_statistics(sqp.problem, "D(ϕ,p)", sqp.directional_derivative)
    add_statistics(sqp.problem, "|p|", norm(sqp.p, Inf))
    add_statistics(sqp.problem, "|J|2", norm(sqp.dE, 2))
    add_statistics(sqp.problem, "|J|inf", norm(sqp.dE, Inf))
    add_statistics(sqp.problem, "inf_pr", sqp.prim_infeas)
    # add_statistics(sqp.problem, "inf_du", dual_infeas)
    add_statistics(sqp.problem, "compl", sqp.compl)
    add_statistics(sqp.problem, "alpha", sqp.alpha)
    add_statistics(sqp.problem, "iter_time", time() - sqp.start_iter_time)
    add_statistics(sqp.problem, "time_elapsed", time() - sqp.start_time)
end


text:
function parse_input(lines)
    function parse_bag_content(content)
        if content == "no other bags."
            []
        elseif !occursin(",", content)
            child_bag = join(split(content)[2:3], " ")
            [child_bag]
        else
            map(x -> join(split(x)[2:3], " "), split(content, ", "))
        end
    end

    parents = Dict{AbstractString, Set{AbstractString}}()
    for line in lines
        first_part, second_part = split(line, " contain ")
        parent_bag = join(split(first_part)[1:2], " ")
        children = parse_bag_content(second_part)
        for child_bag in children
            vals = get(parents, child_bag, Set{AbstractString}())
            push!(vals, parent_bag)
            parents[child_bag] = vals
        end
    end
    parents
end

function find_possible_outer_bags(parent_bags, start_bag)
    visited_bags = Set{AbstractString}()
    current_bags = Set([start_bag])
    while !isempty(current_bags)
        new_parents = Set{AbstractString}()
        for bag in current_bags
            parents = get(parent_bags, bag, Set())
            for parent in parents
                if parent ∉ visited_bags
                    push!(new_parents, parent)
                end
            end
        end
        current_bags = new_parents
        union!(visited_bags, new_parents)
    end
    visited_bags
end

function main()
    lines = readlines("../resources/input_07.txt")
    parents = parse_input(lines)
    result = find_possible_outer_bags(parents, "shiny gold")
    println(length(result))
end

main()


text:
@testset "1631.path-with-minimum-effort.jl" begin
    @test minimum_effort_path([[1, 2, 2], [3, 8, 2], [5, 3, 5]]) == 2
    @test minimum_effort_path([[1, 2, 3], [3, 8, 4], [5, 3, 5]]) == 1
    @test minimum_effort_path([
        [1, 2, 1, 1, 1], [1, 2, 1, 2, 1], [1, 2, 1, 2, 1], [1, 2, 1, 2, 1], [1, 1, 1, 2, 1]
    ]) == 0
    @test minimum_effort_path([[1, 10, 6, 7, 9, 10, 4, 9]]) == 9
end


text:
# Calculate and plot in and out degrees of generated Kronecker Graph
using JLD,PyPlot

# Associative arrays to accumulate in and out degrees
Aall = Assoc("","","")

# Iterate through files
(@isdefined Nfile) || (Nfile = 8)
for i = 1:Nfile
    global Aall
    ptime = @elapsed begin
        fname = joinpath(Base.source_dir(),"data", string(i))

        # Load associative array
        A = load(fname*".A.jld")["A"]
        #A = loadassoc(fname*".A.jld")

        # Create associative arrays and accumulate degrees
        Aall = A + Aall
    end

    println("Sum Time: ", ptime)
    println("Edges/sec: ", string(nnz(A)/ptime))
end

# Plot out degree distribution
figure()
loglog(convert(Array,OutDegree(Aall)'),"o")
xlabel("out degree")

# Plot out degree distribution
figure()
loglog(convert(Array,InDegree(Aall)'),"o")
xlabel("in degree")


text:
"""
Calculates the optimal feedback gains K,d as well as the 2nd Order approximation of the
Cost-to-Go, using a backward Riccati-style recursion. (non-allocating)
"""
function backwardpass!(solver::iLQRSolver{T,QUAD,L,O,n,n̄,m}) where {T,QUAD<:QuadratureRule,L,O,n,n̄,m}
	N = solver.N

    # Objective
    obj = solver.obj
    model = solver.model

    # Extract variables
    Z = solver.Z; K = solver.K; d = solver.d;
    G = solver.G
    S = solver.S
	Quu_reg = solver.Quu_reg
	Qux_reg = solver.Qux_reg

    # Terminal cost-to-go
	Q = solver.E[N]
    S[N].Q .= Q.Q
    S[N].q .= Q.q

    # Initialize expecte change in cost-to-go
    ΔV = @SVector zeros(2)

	k = N-1
    while k > 0
        ix = Z[k]._x
        iu = Z[k]._u

		# Get error state expanions
		fdx,fdu = TO.error_expansion(solver.D[k], model)
		cost_exp = solver.E[k]
		Q = solver.Q_tmp

		# Calculate action-value expansion
		_calc_Q!(Q, cost_exp, S[k+1], fdx, fdu, S[k])

		# Regularization
		Quu_reg .= Q.R #+ solver.ρ[1]*I
		Quu_reg .+= solver.ρ[1]*Diagonal(@SVector ones(m))
		Qux_reg .= Q.H

	    if solver.opts.bp_reg
	        vals = eigvals(Hermitian(Quu_reg))
	        if minimum(vals) <= 0
	            @warn "Backward pass regularized"
	            regularization_update!(solver, :increase)
	            k = N-1
	            ΔV = @SVector zeros(2)
	            continue
	        end
	    end

        # Compute gains
		_calc_gains!(K[k], d[k], Quu_reg, Qux_reg, Q.r)

		# Calculate cost-to-go (using unregularized Quu and Qux)
		ΔV += _calc_ctg!(S[k], Q, K[k], d[k])

        k -= 1
    end

    regularization_update!(solver, :decrease)

    return ΔV

end

function static_backwardpass!(solver::iLQRSolver{T,QUAD,L,O,n,n̄,m}, grad_only=false) where {T,QUAD<:QuadratureRule,L,O,n,n̄,m}
	N = solver.N

    # Objective
    obj = solver.obj
    model = solver.model

    # Extract variables
    Z = solver.Z; K = solver.K; d = solver.d;
    G = solver.G
    S = solver.S
	Quu_reg = SMatrix(solver.Quu_reg)
	Qux_reg = SMatrix(solver.Qux_reg)

    # Terminal cost-to-go
	# Q = error_expansion(solver.Q[N], model)
	Q = solver.E[N]
	Sxx = SMatrix(Q.Q)
	Sx = SVector(Q.q)

	if solver.opts.save_S
		S[end].Q .= Sxx
		S[end].q .= Sx
	end

    # Initialize expected change in cost-to-go
	ΔV = @SVector zeros(2)

	k = N-1
    while k > 0
        ix = Z[k]._x
        iu = Z[k]._u

		# Get error state expanions
		fdx,fdu = TO.error_expansion(solver.D[k], model)
		fdx,fdu = SMatrix(fdx), SMatrix(fdu)
		Q = TO.static_expansion(solver.E[k])

		# Calculate action-value expansion
		Q = _calc_Q!(Q, Sxx, Sx, fdx, fdu, grad_only)

		# Save Q
		solver.Q[k].Q .= Q.xx
		solver.Q[k].R .= Q.uu
		solver.Q[k].H .= Q.ux
		solver.Q[k].q .= Q.x
		solver.Q[k].r .= Q.u

		# Regularization
		Quu_reg, Qux_reg = _bp_reg!(Q, fdx, fdu, solver.ρ[1], solver.opts.bp_reg_type)

	    if solver.opts.bp_reg
	        vals = eigvals(Hermitian(Quu_reg))
	        if minimum(vals) <= 0
	            @warn "Backward pass regularized"
	            regularization_update!(solver, :increase)
	            k = N-1
	            ΔV = @SVector zeros(2)
	            continue
	        end
	    end

        # Compute gains
		K_, d_ = _calc_gains!(K[k], d[k], Quu_reg, Qux_reg, Q.u, grad_only)

		# Calculate cost-to-go (using unregularized Quu and Qux)
		Sxx, Sx, ΔV_ = _calc_ctg!(Q, K_, d_, grad_only)
		# k >= N-2 && println(diag(Sxx))
		if solver.opts.save_S
			S[k].Q .= Sxx
			S[k].q .= Sx
			# S[k].c = ΔV_[1]
		end
		ΔV += ΔV_
        k -= 1
    end

    regularization_update!(solver, :decrease)

    return ΔV
end

function _bp_reg!(Quu_reg::SizedMatrix{m,m}, Qux_reg, Q, fdx, fdu, ρ, ver=:control) where {m}
    if ver == :state
        Quu_reg .= Q.uu #+ solver.ρ[1]*fdu'fdu
		mul!(Quu_reg, Transpose(fdu), fdu, ρ, 1.0)
        Qux_reg .= Q.ux #+ solver.ρ[1]*fdu'fdx
		mul!(Qux_reg, fdu', fdx, ρ, 1.0)
    elseif ver == :control
        Quu_reg .= Q.uu #+ solver.ρ[1]*I
		Quu_reg .+= ρ*Diagonal(@SVector ones(m))
        Qux_reg .= Q.ux
    end
end

function _bp_reg!(Q, fdx, fdu, ρ, ver=:control)
    if ver == :state
		Quu_reg = Q.uu + ρ * fdu'fdu
		Qux_reg = Q.ux + ρ * fdu'fdx
    elseif ver == :control
		Quu_reg = Q.uu + ρ * I
        Qux_reg = Q.ux
    end

	Quu_reg, Qux_reg
end

function _calc_Q!(Q, cost_exp, S1, fdx, fdu, Q_tmp)
	# Compute the cost-to-go, stashing temporary variables in S[k]
	# Qx =  Q.x[k] + fdx'S.x[k+1]
	mul!(Q.q, Transpose(fdx), S1.q)
	Q.q .+= cost_exp.q

    # Qu =  Q.u[k] + fdu'S.x[k+1]
	mul!(Q.r, Transpose(fdu), S1.q)
	Q.r .+= cost_exp.r

    # Qxx = Q.xx[k] + fdx'S.xx[k+1]*fdx
	mul!(Q_tmp.Q, Transpose(fdx), S1.Q)
	mul!(Q.Q, Q_tmp.Q, fdx)
	Q.Q .+= cost_exp.Q

    # Quu = Q.uu[k] + fdu'S.xx[k+1]*fdu
	mul!(Q_tmp.H, Transpose(fdu), S1.Q)
	mul!(Q.R, Q_tmp.H, fdu)
	Q.R .+= cost_exp.R

    # Qux = Q.ux[k] + fdu'S.xx[k+1]*fdx
	mul!(Q_tmp.H, Transpose(fdu), S1.Q)
	mul!(Q.H, Q_tmp.H, fdx)
	Q.H .+= cost_exp.H

	return nothing
end

function _calc_Q!(Q::TO.StaticExpansion, Sxx, Sx, fdx::SMatrix, fdu::SMatrix, grad_only=false)
	Qx = Q.x + fdx'Sx
	Qu = Q.u + fdu'Sx
	if grad_only
		Qxx = Q.xx
		Quu = Q.uu
		Qux = Q.ux
	else
		Qxx = Q.xx + fdx'Sxx*fdx
		Quu = Q.uu + fdu'Sxx*fdu
		Qux = Q.ux + fdu'Sxx*fdx
	end
	TO.StaticExpansion(Qx,Qxx,Qu,Quu,Qux)
end

function _calc_gains!(K::SizedArray, d::SizedArray, Quu::SizedArray, Qux::SizedArray, Qu)
	LAPACK.potrf!('U',Quu.data)
	K .= Qux
	d .= Qu
	LAPACK.potrs!('U', Quu.data, K.data)
	LAPACK.potrs!('U', Quu.data, d.data)
	K .*= -1
	d .*= -1
	# return K,d
end

function _calc_gains!(K, d, Quu::SMatrix, Qux::SMatrix, Qu::SVector, grad_only=false)
	if grad_only
		K_ = SMatrix(K)
	else
		K_ = -Quu\Qux
		K .= K_
	end
	d_ = -Quu\Qu
	d .= d_
	return K_,d_
end

function _calc_ctg!(S, Q, K, d)
	# S.x[k]  =  Qx + K[k]'*Quu*d[k] + K[k]'* Qu + Qux'd[k]
	tmp1 = S.r
	S.q .= Q.q
	mul!(tmp1, Q.R, d)
	mul!(S.q, Transpose(K), tmp1, 1.0, 1.0)
	mul!(S.q, Transpose(K), Q.r, 1.0, 1.0)
	mul!(S.q, Transpose(Q.H), d, 1.0, 1.0)

	# S.xx[k] = Qxx + K[k]'*Quu*K[k] + K[k]'*Qux + Qux'K[k]
	tmp2 = S.H
	S.Q .= Q.Q
	mul!(tmp2, Q.R, K)
	mul!(S.Q, Transpose(K), tmp2, 1.0, 1.0)
	mul!(S.Q, Transpose(K), Q.H, 1.0, 1.0)
	mul!(S.Q, Transpose(Q.H), K, 1.0, 1.0)
	transpose!(Q.Q, S.Q)
	S.Q .+= Q.Q
	S.Q .*= 0.5

    # calculated change is cost-to-go over entire trajectory
	t1 = dot(d, Q.r)
	mul!(Q.r, Q.R, d)
	t2 = 0.5*dot(d, Q.r)
    return @SVector [t1, t2]
end

function _calc_ctg!(Q::TO.StaticExpansion, K::SMatrix, d::SVector, grad_only::Bool=false)
	Sx = Q.x + K'Q.uu*d + K'Q.u + Q.ux'd
	if grad_only
		Sxx = Q.xx
	else
		Sxx = Q.xx + K'Q.uu*K + K'Q.ux + Q.ux'K
		Sxx = 0.5*(Sxx + Sxx')
	end
	t1 = d'Q.u
	t2 = 0.5*d'Q.uu*d
	return Sxx, Sx, @SVector [t1, t2]
end


text:
using SatisfiabilityInterface, Symbolics

"Constraints representing different colours for each edge of a graph"
different_colours(E, c) = [c[i] ≠ c[j] for (i, j) in E]

"k is the number of colours"
function graph_colouring_problem(V, E, k=3)

    colours = [:red, :green, :yellow, :blue, :black][1:k]
    # c = [Num(Variable(:c, i)) for i in 1:length(V)]   # colour variables

    @variables c[1:length(V)]

    constraints =
    [
        [c[i] ∈ colours for i in 1:length(V)]

        [c[i] ≠ c[j] for (i, j) in E]  #  different_colours(E, c)
    ]

    return DiscreteCSP(constraints)
end

# cs = [Num(Variable(:c, i)) for i in 1:length(V)]   # colour variables

# constraints = [
#     [c ∈ colours for c in cs]
#     different_colours(E, cs)
# ]

V = [1, 2, 3]  # vertices
E = [(1, 2), (2, 3)]  # edges

prob = graph_colouring_problem(V, E, 2)
status, results = solve(prob)

# final_colours = [second(results[k]) for k in prob.original_vars]

# # check that it satisfies the constraint:
# all(different_neighbours(E, final_colours))

## Ring graph

function ring_graph(n=11)
    V = 1:n
    E = [(i, mod1((i+1), n)) for i in 1:n]

    return V, E
end

V, E = ring_graph(11)

prob = graph_colouring_problem(V, E, 2)
status, results = solve(prob)

prob = graph_colouring_problem(V, E, 3)
status, results = solve(prob)
status==:sat

# final_colours = [results[k] for k in prob.variables]
# all(different_neighbours(E, final_colours))

k = 3
colours = [:red, :green, :yellow, :blue, :black][1:k]

c = [Num(Variable(:c, i)) for i in 1:length(V)]   # colour variables

constraints =
    [
        [c[i] ∈ colours for i in 1:length(V)]

        [c[i] ≠ c[j] for (i, j) in E]  #  different_colours(E, c)
    ]

constraints

prob = ConstraintSatisfactionProblem(constraints)
prob2 = DiscreteCSP(prob)

solve(prob2)

k = 3


text:
module JLLPrefixes
using Pkg, Pkg.Artifacts, Base.BinaryPlatforms

export collect_artifact_paths, symlink_artifact_paths, unsymlink_artifact_paths, copy_artifact_paths

# Bring in helpers for git repositories
include("libgit2_utils.jl")

# Bring in helpers to deal with JLL packages
include("pkg_utils.jl")

# Bring in helpers to deal with symlink nests
include("symlink_utils.jl")

# Only update the registry once per session, by default
const _registry_updated = Ref(false)

"""
    collect_artifact_paths(dependencies::Vector;
                           platform = HostPlatform(),
                           verbose = false)

Collect all (recursive) JLL dependency artifact paths for the given `platform`.  Returns
a dictionary mapping each (recursive) dependency to its set of artifact paths, which
can then be flattened and given to other tools, such as `symlink_artifact_paths()` or
`copy_artifact_paths()`.

Because the dependencies can be specified as a `PkgSpec`, it is possible to request
particular versions of a package just as you would with `Pkg.add()`.

The `platform` keyword argument allows for collecting artifacts for a foreign platform,
as well as a different Julia version.  This is especially useful for packages that are
stdlibs, and thus locked to a single version based on the Julia version.
"""
function collect_artifact_paths(dependencies::Vector{PkgSpec};
                                platform::AbstractPlatform = HostPlatform(),
                                project_dir::AbstractString = mktempdir(),
                                update_registry::Bool = _registry_updated[],
                                verbose::Bool = false)
    # We occasionally generate "illegal" package specs, where we provide both version and tree hash.
    # we trust the treehash over the version, so drop the version for any that exists here:
    function filter_redundant_version(p::PkgSpec)
        if p.version !== nothing && p.tree_hash !== nothing
            return Pkg.Types.PackageSpec(;name=p.name, tree_hash=p.tree_hash, repo=p.repo)
        end
        return p
    end
    dependencies = filter_redundant_version.(dependencies)
    dependencies_names = [d.name for d in dependencies]

    # Get julia version specificity, if it exists, from the `Platform` object
    julia_version = nothing
    if haskey(platform, "julia_version")
        julia_version = VersionNumber(platform["julia_version"])
    end

    # This is what we will eventually return
    artifact_paths = Dict{PkgSpec, Vector{String}}()

    # We're going to create a project and install all dependent packages within
    # it, then create symlinks from those installed products to our build prefix
    deps_project = joinpath(project_dir, "Project.toml")
    Pkg.activate(deps_project) do
        ctx = Pkg.Types.Context(;julia_version)
        pkg_io = verbose ? stdout : devnull

        # Update registry first, in case the jll packages we're looking for have just been registered/updated
        if update_registry
            Pkg.Registry.update(
                [Pkg.RegistrySpec(uuid = "23338594-aafe-5451-b93e-139f81909106")];
                io=pkg_io,
            )
            _registry_updated[] = true
        end

        # Add all dependencies to our project
        Pkg.add(ctx, dependencies; platform=platform, io=pkg_io)

        # Ony Julia v1.6, `Pkg.add()` doesn't mutate `dependencies`, so we can't use the `UUID`
        # that was found during resolution there.  Instead, we'll make use of `ctx.env` to figure
        # out the UUIDs of all our packages.
        dependency_uuids = Set([uuid for (uuid, pkg) in ctx.env.manifest if pkg.name ∈ dependencies_names])

        # Some JLLs are also standard libraries that may be present in the manifest because
        # they were pulled by other stdlibs (e.g. through dependence on `Pkg`), not beacuse
        # they were actually required for this package. Filter them out if they're present
        # in the manifest but aren't direct dependencies or dependencies of other JLLS.
        installed_jll_uuids = collect_jll_uuids(ctx.env.manifest, dependency_uuids)
        installed_jlls = [
            PkgSpec(;
                name=pkg.name,
                uuid,
                tree_hash=pkg.tree_hash,
            ) for (uuid, pkg) in ctx.env.manifest if uuid ∈ installed_jll_uuids
        ]

        # Check for stdlibs lurking in the installed JLLs
        stdlib_pkgspecs = PkgSpec[]
        for dep in installed_jlls
            # If the `tree_hash` is `nothing`, then this JLL was treated as an stdlib
            if dep.tree_hash === nothing
                # Figure out what version this stdlib _should_ be at for this version
                dep.version = stdlib_version(dep.uuid, julia_version)

                # Interrogate the registry to determine the correct treehash
                Pkg.Operations.load_tree_hash!(ctx.registries, dep, nothing)

                # We'll still use `Pkg.add()` to install the version we want, even though
                # we've used the above two lines to figure out the treehash, so construct
                # an addable spec that will get the correct bits down on disk.
                push!(stdlib_pkgspecs, get_addable_spec(dep.name, dep.version; verbose))
            end
        end

        # Re-install stdlib dependencies, but this time with `julia_version = nothing`
        if !isempty(stdlib_pkgspecs)
            Pkg.add(ctx, stdlib_pkgspecs; io=pkg_io, julia_version=nothing)
        end

        # Load their Artifacts.toml files
        for dep in installed_jlls
            dep_path = Pkg.Operations.find_installed(dep.name, dep.uuid, dep.tree_hash)

            # Skip dependencies that didn't get installed, but warn as this should never happen
            if dep_path === nothing
                @warn("Dependency $(dep.name) not installed, despite our best efforts!")
                continue
            end

            # Load the Artifacts.toml file
            artifacts_toml = joinpath(dep_path, "Artifacts.toml")
            if !isfile(artifacts_toml)
                # Try `StdlibArtifacts.toml` instead
                artifacts_toml = joinpath(dep_path, "StdlibArtifacts.toml")
                if !isfile(artifacts_toml)
                    @warn("Dependency $(dep.name) does not have an (Stdlib)Artifacts.toml in $(dep_path)!")
                    continue
                end
            end

            # If the artifact is available for the given platform, make sure it
            # is also installed.  It may not be the case for lazy artifacts.
            meta = artifact_meta(dep.name[1:end-4], artifacts_toml; platform=platform)
            if meta === nothing
                # This only gets printed if we're verbose, as this can be kind of common
                if verbose
                    @warn("Dependency $(dep.name) does not have a mapping for artifact $(dep.name[1:end-4]) for platform $(triplet(platform))")
                end
                continue
            end

            # NOTE: We are here assuming that each dependency should download only a single
            # artifact with a specific name (e.g. `Zlib_jll` should download the `Zlib`
            # artifact).  This is currently how JLLs work, but it may not always be the case!
            # We need to come up with some way of allowing the user to specify _which_ artifacts
            # they want downloaded!
            ensure_artifact_installed(dep.name[1:end-4], meta, artifacts_toml; platform=platform)

            # Copy the artifact from the global installation location into this build-specific artifacts collection
            src_path = Pkg.Artifacts.artifact_path(Base.SHA1(meta["git-tree-sha1"]))

            # Keep track of our dep paths for later symlinking
            if !haskey(artifact_paths, dep)
                artifact_paths[dep] = String[]
            end
            push!(artifact_paths[dep], src_path)
        end
    end

    return artifact_paths
end
function collect_artifact_paths(pkg_names::Vector{<:AbstractString}; kwargs...)
    return collect_artifact_paths(
        [PackageSpec(; name) for name in pkg_names];
        kwargs...
    )
end

# Helper function for throwing away the `PkgSpec` information
function flatten_artifact_paths(d::Dict{PkgSpec, Vector{String}})
    return vcat(values(d)...)
end

"""
    symlink_artifact_paths(dest::String, artifact_paths)

Symlinks all files from the given `artifact_paths` to the given `dest`.
Provides a merged symlink nest for smashing multiple JLLs together into a single
prefix, as a cheaper alternative to `copy_artifact_paths`.
"""
function symlink_artifact_paths(dest::AbstractString, artifact_paths::Vector{String}; verbose::Bool = true)
    for artifact_path in artifact_paths
        symlink_tree(dest, artifact_path; verbose)
    end
end
function symlink_artifact_paths(dest::AbstractString, artifact_paths::Dict{PkgSpec, Vector{String}}; kwargs...)
    return symlink_artifact_paths(dest, flatten_artifact_paths(artifact_paths); kwargs...)
end

"""
    unsymlink_artifact_paths(dest::String, artifact_paths)

Removes all symlinks in `dest` that originate from `artifact_paths`.  Useful if
there is a mixture of symlinked artifacts and other files, and you want to remove
all symlinks from the artifact paths.
"""
function unsymlink_artifact_paths(dest::AbstractString, artifact_paths::Vector{String})
    for artifact_path in artifact_paths
        unsymlink_tree(dest, artifact_path)
    end
end
function unsymlink_artifact_paths(dest::AbstractString, artifact_paths::Dict{PkgSpec, Vector{String}})
    return unsymlink_artifact_paths(dest, flatten_artifact_paths(artifact_paths))
end

"""
    copy_artifact_paths(dest::String, artifact_paths)
"""
function copy_artifact_paths(dest::AbstractString, artifact_paths::Vector{String})
    for artifact_path in artifact_paths
        copy_tree(dest, artifact_path)
    end
end
function copy_artifact_paths(dest::AbstractString, artifact_paths::Dict{PkgSpec, Vector{String}})
    return copy_artifact_paths(dest, flatten_artifact_paths(artifact_paths))
end

end # module JLLPrefixes


text:
function HaloOrbit(arg0::CR3BPSystem, arg1::PVCoordinates, arg2::jdouble)
    return HaloOrbit((CR3BPSystem, PVCoordinates, jdouble), arg0, arg1, arg2)
end

function HaloOrbit(arg0::RichardsonExpansion, arg1::jdouble, arg2::LibrationOrbitFamily)
    return HaloOrbit((RichardsonExpansion, jdouble, LibrationOrbitFamily), arg0, arg1, arg2)
end

function apply_differential_correction(obj::LibrationOrbit)
    return jcall(obj, "applyDifferentialCorrection", void, ())
end

function apply_differential_correction(obj::LibrationOrbit, arg0::AttitudeProvider, arg1::TimeScale)
    return jcall(obj, "applyDifferentialCorrection", void, (AttitudeProvider, TimeScale), arg0, arg1)
end

function get_initial_pv(obj::LibrationOrbit)
    return jcall(obj, "getInitialPV", PVCoordinates, ())
end

function get_manifolds(obj::LibrationOrbit, arg0::SpacecraftState, arg1::jboolean)
    return jcall(obj, "getManifolds", PVCoordinates, (SpacecraftState, jboolean), arg0, arg1)
end

function get_orbital_period(obj::LibrationOrbit)
    return jcall(obj, "getOrbitalPeriod", jdouble, ())
end


text:
# This file is auto-generated by AWSMetadata.jl
using AWS
using AWS.AWSServices: elastic_inference
using AWS.Compat
using AWS.UUIDs

"""
    describe_accelerator_offerings(location_type)
    describe_accelerator_offerings(location_type, params::Dict{String,<:Any})

 Describes the locations in which a given accelerator type or set of types is present in a
given region.

# Arguments
- `location_type`:  The location type that you want to describe accelerator type offerings
  for. It can assume the following values: region: will return the accelerator type offering
  at the regional level. availability-zone: will return the accelerator type offering at the
  availability zone level. availability-zone-id: will return the accelerator type offering at
  the availability zone level returning the availability zone id.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"acceleratorTypes"`:  The list of accelerator types to describe.
"""
describe_accelerator_offerings(locationType; aws_config::AbstractAWSConfig=global_aws_config()) = elastic_inference("POST", "/describe-accelerator-offerings", Dict{String, Any}("locationType"=>locationType); aws_config=aws_config)
describe_accelerator_offerings(locationType, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = elastic_inference("POST", "/describe-accelerator-offerings", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("locationType"=>locationType), params)); aws_config=aws_config)

"""
    describe_accelerator_types()
    describe_accelerator_types(params::Dict{String,<:Any})

 Describes the accelerator types available in a given region, as well as their
characteristics, such as memory and throughput.

"""
describe_accelerator_types(; aws_config::AbstractAWSConfig=global_aws_config()) = elastic_inference("GET", "/describe-accelerator-types"; aws_config=aws_config)
describe_accelerator_types(params::AbstractDict{String, Any}; aws_config::AbstractAWSConfig=global_aws_config()) = elastic_inference("GET", "/describe-accelerator-types", params; aws_config=aws_config)

"""
    describe_accelerators()
    describe_accelerators(params::Dict{String,<:Any})

 Describes information over a provided set of accelerators belonging to an account.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"acceleratorIds"`:  The IDs of the accelerators to describe.
- `"filters"`:  One or more filters. Filter names and values are case-sensitive. Valid
  filter names are: accelerator-types: can provide a list of accelerator type names to filter
  for. instance-id: can provide a list of EC2 instance ids to filter for.
- `"maxResults"`:  The total number of items to return in the command's output. If the
  total number of items available is more than the value specified, a NextToken is provided
  in the command's output. To resume pagination, provide the NextToken value in the
  starting-token argument of a subsequent command. Do not use the NextToken response element
  directly outside of the AWS CLI.
- `"nextToken"`:  A token to specify where to start paginating. This is the NextToken from
  a previously truncated response.
"""
describe_accelerators(; aws_config::AbstractAWSConfig=global_aws_config()) = elastic_inference("POST", "/describe-accelerators"; aws_config=aws_config)
describe_accelerators(params::AbstractDict{String, Any}; aws_config::AbstractAWSConfig=global_aws_config()) = elastic_inference("POST", "/describe-accelerators", params; aws_config=aws_config)

"""
    list_tags_for_resource(resource_arn)
    list_tags_for_resource(resource_arn, params::Dict{String,<:Any})

 Returns all tags of an Elastic Inference Accelerator.

# Arguments
- `resource_arn`:  The ARN of the Elastic Inference Accelerator to list the tags for.

"""
list_tags_for_resource(resourceArn; aws_config::AbstractAWSConfig=global_aws_config()) = elastic_inference("GET", "/tags/$(resourceArn)"; aws_config=aws_config)
list_tags_for_resource(resourceArn, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = elastic_inference("GET", "/tags/$(resourceArn)", params; aws_config=aws_config)

"""
    tag_resource(resource_arn, tags)
    tag_resource(resource_arn, tags, params::Dict{String,<:Any})

 Adds the specified tags to an Elastic Inference Accelerator.

# Arguments
- `resource_arn`:  The ARN of the Elastic Inference Accelerator to tag.
- `tags`:  The tags to add to the Elastic Inference Accelerator.

"""
tag_resource(resourceArn, tags; aws_config::AbstractAWSConfig=global_aws_config()) = elastic_inference("POST", "/tags/$(resourceArn)", Dict{String, Any}("tags"=>tags); aws_config=aws_config)
tag_resource(resourceArn, tags, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = elastic_inference("POST", "/tags/$(resourceArn)", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("tags"=>tags), params)); aws_config=aws_config)

"""
    untag_resource(resource_arn, tag_keys)
    untag_resource(resource_arn, tag_keys, params::Dict{String,<:Any})

 Removes the specified tags from an Elastic Inference Accelerator.

# Arguments
- `resource_arn`:  The ARN of the Elastic Inference Accelerator to untag.
- `tag_keys`:  The list of tags to remove from the Elastic Inference Accelerator.

"""
untag_resource(resourceArn, tagKeys; aws_config::AbstractAWSConfig=global_aws_config()) = elastic_inference("DELETE", "/tags/$(resourceArn)", Dict{String, Any}("tagKeys"=>tagKeys); aws_config=aws_config)
untag_resource(resourceArn, tagKeys, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = elastic_inference("DELETE", "/tags/$(resourceArn)", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("tagKeys"=>tagKeys), params)); aws_config=aws_config)


text:
export testLossFunction

function testLossFunction(S::SoftmaxLoss{T};nf=10,nex=20,nc=7,out=false) where {T<:Number}

Y   = randn(T,nf,nex)/10
nw2 = (S.addBias) ? nf+1 : nf
W   = randn(T,nc,nw2)/10
C   = full(sparse(full(rand(1:nc,nex)),collect(1:nex),ones(nex),nc,nex))
C   = convert.(T,C)

@testset "test dW" begin
function testdWFun(x,v=[])
    x = reshape(x,nc,nw2)
    if !(isempty(v))
        F, para, dWF, = getMisfit(S,x,Y,C,false,true)
        return F, dot(dWF,v)
    else
        return getMisfit(S,x,Y,C,false,false)[1]
    end
end
chkDer, = checkDerivative(testdWFun,copy(vec(W)),out=out)
@test chkDer
end

@testset "test d2W" begin
function testd2WFun(x,v=[])
    x = reshape(x,nc,nw2)
    if !(isempty(v))
        F, para, dWF,d2F = getMisfit(S,x,Y,C,false,true)
        return vec(dWF),d2F*vec(v)
    else
        return vec(getMisfit(S,x,Y,C,false,true)[3])
    end
end
chkDer, = checkDerivative(testd2WFun,copy(vec(W)),out=out)
@test chkDer
end

@testset "test dY" begin
function testdYFun(x,v=[])
    x = reshape(x,nf,nex)
    if !(isempty(v))
        F, para, dWF,d2WF,dYF = getMisfit(S,W,x,C,true,false)
        return F, dot(dYF,v)
    else
        return getMisfit(S,W,x,C,false,false)[1]
    end
end
chkDer, = checkDerivative(testdYFun,copy(vec(Y)),out=out)
@test chkDer
end

@testset "test d2Y" begin
function testd2YFun(x,v=[])
    x = reshape(x,nf,nex)
    if !(isempty(v))
        F, para, dWF,d2WF,dYF,d2YF = getMisfit(S,W,x,C,true,false)
        return vec(dYF),d2YF*vec(v)
    else
        return vec(getMisfit(S,W,x,C,true,false)[5])
    end
end
chkDer, = checkDerivative(testd2YFun,copy(vec(Y)),out=out)
@test chkDer
end
end


text:
using Test
using BinaryBuilderBase
using BinaryBuilderBase: dlext, exeext

@testset "Wrappers utilities" begin
    @test nbits(Linux(:i686)) == 32
    @test nbits(ExtendedPlatform(Linux(:x86_64); march="avx")) == 64
    @test nbits(Linux(:armv7l)) == 32
    @test nbits(ExtendedPlatform(Linux(:aarch64); cuda="10.1")) == 64
    @test nbits(Linux(:powerpc64le)) == 64
    @test nbits(AnyPlatform()) == 64
    @test_throws ErrorException nbits(UnknownPlatform())

    @test proc_family(Linux(:i686)) == :intel
    @test proc_family(ExtendedPlatform(Linux(:x86_64); march="avx")) == :intel
    @test proc_family(Linux(:armv7l)) == :arm
    @test proc_family(ExtendedPlatform(Linux(:aarch64); cuda="10.1")) == :arm
    @test proc_family(Linux(:powerpc64le)) == :power
    @test proc_family(AnyPlatform()) == :intel
    @test_throws ErrorException proc_family(UnknownPlatform())

    @test dlext(Linux(:i686)) == "so"
    @test dlext(FreeBSD(:x86_64)) == "so"
    @test dlext(MacOS(:x86_64)) == "dylib"
    @test dlext(Windows(:i686)) == "dll"
    @test dlext(ExtendedPlatform(Linux(:x86_64); march="avx512")) == "so"
    @test_throws ErrorException dlext(AnyPlatform())
    @test_throws ErrorException dlext(UnknownPlatform())

    @test exeext(Linux(:i686)) == ""
    @test exeext(FreeBSD(:x86_64)) == ""
    @test exeext(MacOS(:x86_64)) == ""
    @test exeext(Windows(:i686)) == ".exe"
    @test exeext(ExtendedPlatform(Linux(:x86_64); march="avx512")) == ""
    @test_throws ErrorException exeext(AnyPlatform())
    @test_throws ErrorException exeext(UnknownPlatform())

    @test aatriplet(ExtendedPlatform(Linux(:x86_64, libc=:glibc, compiler_abi=CompilerABI(libgfortran_version=v"4.0.0")); march="avx", cuda="9.2")) == "x86_64-linux-gnu"
end

# Are we using docker? If so, test that the docker runner works...
@testset "Runner utilities" begin
    # Test that is_ecryptfs works for something we're certain isn't encrypted
    if isdir("/proc")
        isecfs = (false, "/proc/")
        @test_logs (:info, "Checking to see if /proc/ is encrypted...") @test BinaryBuilderBase.is_ecryptfs("/proc"; verbose=true) == isecfs
        @test_logs (:info, "Checking to see if /proc/ is encrypted...") @test BinaryBuilderBase.is_ecryptfs("/proc/"; verbose=true) == isecfs
        @test_logs (:info, "Checking to see if /proc/not_a_file is encrypted...") @test BinaryBuilderBase.is_ecryptfs("/proc/not_a_file"; verbose=true) == isecfs
    else
        @test_logs (:info, "Checking to see if /proc/ is encrypted...") @test BinaryBuilderBase.is_ecryptfs("/proc"; verbose=true) == (false, "/proc")
        @test_logs (:info, "Checking to see if /proc/ is encrypted...") @test BinaryBuilderBase.is_ecryptfs("/proc/"; verbose=true) == (false, "/proc/")
        @test_logs (:info, "Checking to see if /proc/not_a_file is encrypted...") @test BinaryBuilderBase.is_ecryptfs("/proc/not_a_file"; verbose=true) == (false, "/proc/not_a_file")
    end

    if isa(preferred_runner(), BinaryBuilderBase.DockerRunner)
        @testset "Docker image importing" begin
            # First, delete the docker image, in case it already existed
            BinaryBuilderBase.delete_docker_image()

            # Next, import it and ensure that doesn't throw
            rootfs = first(BinaryBuilderBase.choose_shards(platform))
            mktempdir() do dir
                @test BinaryBuilderBase.import_docker_image(rootfs, dir; verbose=true) === nothing
            end

            # Test that deleting the docker image suceeds, now that we know
            # it exists
            @test BinaryBuilderBase.delete_docker_image()
        end
    end

    @testset "hello world" begin
        mktempdir() do dir
            ur = preferred_runner()(dir; platform=Linux(:x86_64; libc=:musl))
            iobuff = IOBuffer()
            @test run(ur, `/bin/bash -c "echo test"`, iobuff)
            seek(iobuff, 0)
            # Test that we get the output we expect (e.g. the second line is `test`)
            @test split(String(read(iobuff)), "\n")[2] == "test"
        end
    end

    @testset "Compilation" begin
        mktempdir() do dir
            ur = preferred_runner()(dir; platform=Linux(:x86_64; libc=:musl))
            iobuff = IOBuffer()
            @test run(ur, `/bin/bash -c "echo 'int main() {return 0;}' | cc -x c -"`, iobuff; tee_stream=devnull)
            seekstart(iobuff)
            @test split(String(read(iobuff)), "\n")[2] == ""
        end
    end

    @testset "Locking microarchitecture" begin
        mktempdir() do dir
            platform = Linux(:x86_64; libc=:musl)
            cmd = `/bin/bash -c "echo 'int main() {return 0;}' | cc -x c -march=native -"`
            ur = preferred_runner()(dir; platform=platform, lock_microarchitecture=true)
            iobuff = IOBuffer()
            @test !run(ur, cmd, iobuff; tee_stream=devnull)
            seekstart(iobuff)
            @test split(String(read(iobuff)), "\n")[2] == "Cannot force an architecture"

            ur = preferred_runner()(dir; platform=platform, lock_microarchitecture=false)
            iobuff = IOBuffer()
            @test run(ur, cmd, iobuff; tee_stream=devnull)
            seekstart(iobuff)
            @test split(String(read(iobuff)), "\n")[2] == ""
        end
    end

    @testset "Unsafe flags" begin
        mktempdir() do dir
            platform = Linux(:x86_64; libc=:musl)
            cmd = `/bin/bash -c "echo 'int main() {return 0;}' | cc -x c -Ofast -"`
            ur = preferred_runner()(dir; platform=platform, allow_unsafe_flags=false)
            iobuff = IOBuffer()
            @test !run(ur, cmd, iobuff; tee_stream=devnull)
            seekstart(iobuff)
            lines = split(String(read(iobuff)), "\n")
            @test lines[2] == "You used one or more of the unsafe flags: -Ofast, -ffast-math, -funsafe-math-optimizations"
            @test lines[3] == "Please repent."

            ur = preferred_runner()(dir; platform=platform, allow_unsafe_flags=true)
            iobuff = IOBuffer()
            @test run(ur, cmd, iobuff; tee_stream=devnull)
            seekstart(iobuff)
            @test split(String(read(iobuff)), "\n")[2] == ""
        end
    end
end


text:
struct ClickHouseClient
    url::URI
end

"""
    ClickHouseClient(url::AbstractString; params...)

# Example
```julia
client = ClickHouseClient("http://localhost:8123")
client = ClickHouseClient("http://localhost:8123", user="...", password="...")
```
"""
function ClickHouseClient(url::AbstractString; params...)
    # TODO: Cleanup
    url = URI(url)
    url = URI(
        scheme = url.scheme,
        userinfo = url.userinfo,
        host = url.host,
        port = url.port,
        path = url.path,
        query = params,
        fragment = url.fragment
    )
    ClickHouseClient(URI(url))
end

# ClickHouseClient()

# TODO: Auth, database, ...
# TODO: Error and exceptions handling

function query(client::ClickHouseClient, query_str)
    query_str = strip(query_str)
    query_str = strip(query_str, [';'])

    if !startswith(query_str, "INSERT INTO")
        query_str = "$(query_str) FORMAT CSVWithNames"
    end

    @debug query_str

    response = HTTP.post(string(client.url), [], query_str)
    resp_str = String(response.body)
    # TODO: Check status code

    @debug response

    CSV.read(IOBuffer(resp_str))
end


text:
############################################################################
### Types

struct ImplySymbolic{
    F <: Any,
    F3 <: Any,
    A <: AbstractArray,
    A3 <: AbstractArray,
    S <: Array{Float64},
    F2 <: Any,
    A2 <: Union{Nothing, AbstractArray}} <: Imply
    imp_fun::F
    gradient_fun::F3
    imp_cov::A
    ∇Σ::A3
    start_val::S
    imp_fun_mean::F2
    imp_mean::A2
end

# for forward diff
mutable struct ImplySymbolicForward{
    F <: Any,
    S <: Array{Float64}} <: Imply
    imp_fun::F
    imp_cov
    start_val::S
end

# for reverse diff and others
mutable struct ImplySymbolicAlloc{
    F <: Any,
    S <: Array{Float64}} <: Imply
    imp_fun::F
    imp_cov
    start_val::S
end

struct ImplySymbolicDefinition{
    F <: Any,
    A <: AbstractArray,
    S <: Array{Float64},
    F2 <: Any,
    A2 <: Union{Nothing, AbstractArray},
    I <: Int64,
    P <: Int64,
    R <: AbstractArray,
    D <: AbstractArray} <: ImplyDefinition
imp_fun::F
imp_cov::A # Array of matrices
start_val::S
imp_fun_mean::F2
imp_mean::A2 # Array of matrices of meanstructure
n_obs::I
n_patterns::P
rows::R
data_def::D
end

struct ImplySymbolicWLS{
    F <: Any,
    F2 <: Any,
    F3 <: Any,
    A <: AbstractArray,
    A2 <: AbstractArray,
    A3 <: Any,
    S <: Array{Float64}} <: Imply
    imp_fun::F
    gradient_fun::F2
    hessian_fun::F3
    imp_cov::A
    ∇Σ::A2
    ∇²Σ::A3
    start_val::S
end

struct ImplySymbolicSWLS{
    F <: Any,
    A <: AbstractArray,
    S <: Array{Float64}} <: Imply
    imp_fun::F
    G::A
    start_val::S
end
############################################################################
### Constructors

function ImplySymbolic(
        A::Spa1,
        S::Spa2,
        F::Spa3,
        parameters,
        start_val;
        M::Spa4 = nothing
            ) where {
            Spa1 <: SparseMatrixCSC,
            Spa2 <: SparseMatrixCSC,
            Spa3 <: SparseMatrixCSC,
            Spa4 <: Union{Nothing, AbstractArray}
            }

    #Model-implied covmat

    invia = neumann_series(A)

    imp_cov_sym = F*invia*S*permutedims(invia)*permutedims(F)

    imp_cov_sym = Array(imp_cov_sym)
    imp_cov_sym = ModelingToolkit.simplify.(imp_cov_sym)

    imp_fun =
        eval(ModelingToolkit.build_function(
            imp_cov_sym,
            parameters
        )[2])

    imp_cov = zeros(size(F)[1], size(F)[1])

    ∇Σ_sym = ModelingToolkit.jacobian(vec(imp_cov_sym), parameters)
    gradient_fun =
        eval(ModelingToolkit.build_function(
                ∇Σ_sym,
                parameters
            )[2])
    ∇Σ = zeros(size(F, 1)^2, size(parameters, 1))
    #imp_cov = Base.invokelatest(imp_fun_, start_val)
    #Model implied mean
    if !isnothing(M)
        imp_mean_sym = F*invia*M
        imp_mean_sym = Array(imp_mean_sym)
        imp_mean_sym = ModelingToolkit.simplify.(imp_mean_sym)

        imp_fun_mean =
            eval(ModelingToolkit.build_function(
                imp_mean_sym,
                parameters
            )[2])

        imp_mean = zeros(size(F)[1])
    else
        imp_fun_mean = nothing
        imp_mean = nothing
    end

    return ImplySymbolic(
        imp_fun,
        gradient_fun,
        imp_cov,
        ∇Σ,
        copy(start_val),
        imp_fun_mean,
        imp_mean
    )
end

function ImplySymbolicWLS(
        A::Spa1,
        S::Spa2,
        F::Spa3,
        parameters::A2,
        start_val::B;
        hessian = false
            ) where {
            Spa1 <: SparseMatrixCSC,
            Spa2 <: SparseMatrixCSC,
            Spa3 <: SparseMatrixCSC,
            A2 <: AbstractArray,
            B <: Array{Float64}
            }

    #Model-implied covmat

    invia = neumann_series(A)

    imp_cov_sym = F*invia*S*permutedims(invia)*permutedims(F)

    imp_cov_sym = Array(imp_cov_sym)
    imp_cov_sym = ModelingToolkit.simplify.(imp_cov_sym)
    imp_cov_sym = imp_cov_sym[tril(trues(size(F, 1), size(F, 1)))]
    parvec = [parameters[i] for i in 1:size(parameters, 1)]
    ∇Σ_sym = ModelingToolkit.jacobian(imp_cov_sym, parvec)

    imp_fun =
        eval(ModelingToolkit.build_function(
            imp_cov_sym,
            parameters
        )[2])

    imp_cov = zeros(Int64(0.5*size(F)[1]*(size(F)[1] + 1)))

    gradient_fun =
        eval(ModelingToolkit.build_function(
            ∇Σ_sym,
            parameters
        )[2])
    ∇Σ = zeros(size(imp_cov_sym, 1), size(parameters, 1))

    if hessian
        n_lower = size(imp_cov_sym, 1)
        n_par = size(parvec, 1)
        ∇²Σ_sym_vec = [ModelingToolkit.sparsehessian(imp_cov_sym[i], parvec) for i = 1:n_lower]
        @variables J[1:n_lower]
        # ∇²Σ = similar_sparse_float.(∇²Σ_sym)
        ∇²Σ_sym = zeros(Num, n_par, n_par)
        for i in 1:n_lower
            ∇²Σ_sym += J[i]*∇²Σ_sym_vec[i]
        end
        # ∇²Σ_sym = simplify.(∇²Σ_sym)
        hessian_fun = eval(ModelingToolkit.build_function(∇²Σ_sym, J, parvec)[2])
        # print(ModelingToolkit.build_function(∇²Σ_sym, J, parvec)[2])
        ∇²Σ = zeros(n_par, n_par)
    else
        hessian_fun = nothing
        ∇²Σ = nothing
    end
    return ImplySymbolicWLS(
        imp_fun,
        gradient_fun,
        hessian_fun,
        imp_cov,
        ∇Σ,
        ∇²Σ,
        copy(start_val)
    )
end

function ImplySymbolicSWLS(
        A::Spa1,
        S::Spa2,
        F::Spa3,
        parameters,
        start_val
            ) where {
            Spa1 <: SparseMatrixCSC,
            Spa2 <: SparseMatrixCSC,
            Spa3 <: SparseMatrixCSC,
            Spa4 <: Union{Nothing, AbstractArray}
            }

    nobs, ntot = size(F)
    nlat = ntot - nobs

    # Selection matrices
    ind = findall(tril(trues(nobs, nobs)))
    Es = []
    for i in ind
        E = zeros(Bool, nobs, nobs)
        E[i] = true
        E = vec(E)
        push!(Es, E)
    end
    L = hcat(Es...)
    K = L*inv(transpose(L)*L)
    K = sparse(K)

    ind = findall(!iszero, S)
    ind = filter(x -> (x[1] >= x[2]), ind)
    Es = []
    for i in ind
        E = zeros(Bool, nobs+nlat, nobs+nlat)
        E[i] = true
        E = vec(E)
        push!(Es, E)
    end
    L_Ω = hcat(Es...)
    L_Ω = sparse(L_Ω)

    #Model-implied covmat
    A = neumann_series(A)

    G = transpose(K)*(kron(F, F)*kron(A, A))*L_Ω
    G = Array(G)
    G = ModelingToolkit.simplify.(G)

    imp_fun =
        eval(ModelingToolkit.build_function(
            G,
            parameters
        )[2])

    G = zeros(size(G))

    return ImplySymbolicSWLS(
        imp_fun,
        G,
        copy(start_val)
    )
end

function ImplySymbolicAlloc(
        A::Spa1,
        S::Spa2,
        F::Spa3,
        parameters,
        start_val
            ) where {
            Spa1 <: SparseMatrixCSC,
            Spa2 <: SparseMatrixCSC,
            Spa3 <: SparseMatrixCSC
            }

    invia = neumann_series(A)

    imp_cov_sym = F*invia*S*invia'*F'

    imp_cov_sym = Array(imp_cov_sym)
    imp_cov_sym .= ModelingToolkit.simplify.(imp_cov_sym)

    imp_fun =
        eval(ModelingToolkit.build_function(
            imp_cov_sym,
            parameters
        )[1])

    imp_cov = rand(size(F)[1], size(F)[1])
    #imp_cov = Base.invokelatest(imp_fun_, start_val)
    return ImplySymbolicAlloc(imp_fun, imp_cov, copy(start_val))
end

function ImplySymbolicForward(
        A::Spa1,
        S::Spa2,
        F::Spa3,
        parameters,
        start_val
            ) where {
            Spa1 <: SparseMatrixCSC,
            Spa2 <: SparseMatrixCSC,
            Spa3 <: SparseMatrixCSC
            }

    invia = neumann_series(A)

    imp_cov_sym = F*invia*S*invia'*F'

    imp_cov_sym = Array(imp_cov_sym)
    imp_cov_sym .= ModelingToolkit.simplify.(imp_cov_sym)

    #imp_cov_sym = DiffEqBase.dualcache(imp_cov_sym)

    imp_fun =
        eval(ModelingToolkit.build_function(
            imp_cov_sym,
            parameters
        )[1])

    imp_cov = DiffEqBase.dualcache(zeros(size(F)[1], size(F)[1]))
    #imp_cov = Base.invokelatest(imp_fun_, start_val)
    return ImplySymbolicForward(imp_fun, imp_cov, copy(start_val))
end

function ImplySymbolicDefinition(
    A::Spa1,
    S::Spa2,
    F::Spa3,
    M::Spa4,
    parameters,
    def_vars,
    start_val,
    data_def
        ) where {
        Spa1 <: SparseMatrixCSC,
        Spa2 <: SparseMatrixCSC,
        Spa3 <: SparseMatrixCSC,
        Spa4 <: Union{Nothing, AbstractArray}
        }

    n_obs = size(data_def, 1)
    n_def_vars = Float64(size(data_def, 2))

    patterns = [data_def[i, :] for i = 1:n_obs]
    remember = Vector{Vector{Float64}}()
    rows = [Vector{Int64}(undef, 0) for i = 1:n_obs]

    for i = 1:size(patterns, 1)
        unknown = true
        for j = 1:size(remember, 1)
            if patterns[i] == remember[j]
                push!(rows[j], i)
                unknown = false
            end
        end
        if unknown
            push!(remember, patterns[i])
            push!(rows[size(remember, 1)], i)
        end
    end

    rows = rows[1:length(remember)]
    n_patterns = size(rows, 1)

    pattern_n_obs = size.(rows, 1)

    #############################################
    #Model-implied covmat
    invia = sem.neumann_series(A)

    imp_cov_sym = F*invia*S*permutedims(invia)*permutedims(F)

    imp_cov_sym = Array(imp_cov_sym)
    imp_cov_sym = ModelingToolkit.simplify.(imp_cov_sym)

    imp_fun =
        eval(ModelingToolkit.build_function(
            imp_cov_sym,
            parameters,
            def_vars
        )[2])

    imp_cov = [zeros(size(F)[1], size(F)[1]) for i = 1:n_patterns]

    #Model implied mean
    imp_mean_sym = F*invia*M
    imp_mean_sym = Array(imp_mean_sym)
    imp_mean_sym = ModelingToolkit.simplify.(imp_mean_sym)

    imp_fun_mean =
        eval(ModelingToolkit.build_function(
            imp_mean_sym,
            parameters,
            def_vars
        )[2])

    imp_mean = [zeros(size(F)[1]) for i = 1:n_patterns]

    data_def = remember

    return ImplySymbolicDefinition(
        imp_fun,
        imp_cov,
        copy(start_val),
        imp_fun_mean,
        imp_mean,
        n_obs,
        n_patterns,
        rows,
        data_def
    )
end

############################################################################
### loss functions

function (imply::ImplySymbolic)(parameters, model)
    imply.imp_fun(imply.imp_cov, parameters)
    if !isnothing(imply.imp_mean)
        imply.imp_fun_mean(imply.imp_mean, parameters)
    end
end

function (imply::ImplySymbolicWLS)(parameters, model)
    imply.imp_fun(imply.imp_cov, parameters)
end

function (imply::ImplySymbolicSWLS)(parameters, model)
    imply.imp_fun(imply.G, parameters)
end

function (imply::ImplySymbolicAlloc)(parameters, model)
    imply.imp_cov = imply.imp_fun(parameters)
end

function (imply::ImplySymbolicForward)(parameters, model)
    imply.imp_cov = DiffEqBase.get_tmp(imply.imp_cov, parameters)
    imply.imp_fun(imply.imp_cov, parameters)
end

function (imply::ImplySymbolicDefinition)(parameters, model)
    for i = 1:imply.n_patterns
        let (cov,
            mean,
            def_vars) =
                (imply.imp_cov[i],
                imply.imp_mean[i],
                imply.data_def[i])
            imply.imp_fun(cov, parameters, def_vars)
            imply.imp_fun_mean(mean, parameters, def_vars)
        end
    end
end


text:
@testset "parser" begin
    b = Bobby.fen_to_bitboard(
        "rnbqkbnr/pppppppp/8/8/8/8/PPPPPPPP/RNBQKBNR w KQkq - 0 1")
    @test b.free == ~b.taken

    @test_throws ArgumentError Bobby.fen_to_bitboard(
        "rnbqkbnr/pppppppp/8/8/8/8/PPPPPPPP/RNBQKBN")
end


text:
struct Cconst{T, S}
	mem::NTuple{S, UInt8}

	Cconst{T}(x) where {T} = new{T, sizeof(T)}(x)
end
Cconst{T, S}(args...; kwargs...) where {T, S} = Cconst{T}(args...; kwargs...)
Cconst{T}(; kwargs...) where {T} = Cconst(T(; kwargs...))
Cconst(::Type{T}) where {T} = Cconst{nonconst(T), sizeof(nonconst(T))}
Cconst(::Type{CA}) where {T, N, CA<:Carray{T, N}} = Carray(Cconst(nonconst(T)), Val(N))
Cconst(x) = x
Cconst(cc::Cconst) = cc
Cconst(ca::Caggregate) = Cconst{typeof(ca)}(getfield(ca, :mem))

nonconst(::Type{T}) where {T} = T
nonconst(::Type{CA}) where {T<:Cconst, N, CA<:Carray{T, N}} = Carray(nonconst(T), Val(N))
nonconst(::Type{CC}) where {T, CC<:Cconst{T}} = T

Base.convert(::Type{T}, cc::Cconst{T}) where {T} = reinterpret(T, getfield(cc, :mem)[1])
Base.convert(::Type{T}, cc::Cconst{T}) where {T<:Caggregate} = T(cc)
Base.sizeof(::Type{CC}) where {T, CC<:Cconst{T}} = sizeof(T)

function Ctypespec(::Type{CC}) where {CC<:Cconst}
	_fix(::Type{Tuple{typ}}) where {typ} = Tuple{Cconst(typ)}
	_fix(::Type{Tuple{typ, bits}}) where {typ, bits} = Tuple{Cconst(typ), bits}
	_fix(::Type{spec}) where {spec<:Ctypespec} = Ctypespec{Cconst(type(spec)), opaque(spec), strategy(spec), specification(spec)}
	return _fix(Ctypespec(nonconst(CC)))
end