text:
###
### **acb_mat.h** -- matrices over the complex numbers
###

### Types, macros and constants

### Memory management
arbcall"void acb_mat_init(acb_mat_t mat, slong r, slong c)"
arbcall"void acb_mat_clear(acb_mat_t mat)"
arbcall"slong acb_mat_allocated_bytes(const acb_mat_t x)"
arbcall"void acb_mat_window_init(acb_mat_t window, const acb_mat_t mat, slong r1, slong c1, slong r2, slong c2)"
arbcall"void acb_mat_window_clear(acb_mat_t window)"

### Conversions
arbcall"void acb_mat_set(acb_mat_t dest, const acb_mat_t src)"
#ni arbcall"void acb_mat_set_fmpz_mat(acb_mat_t dest, const fmpz_mat_t src)"
#ni arbcall"void acb_mat_set_round_fmpz_mat(acb_mat_t dest, const fmpz_mat_t src, slong prec)"
#ni arbcall"void acb_mat_set_fmpq_mat(acb_mat_t dest, const fmpq_mat_t src, slong prec)"
arbcall"void acb_mat_set_arb_mat(acb_mat_t dest, const arb_mat_t src)"
arbcall"void acb_mat_set_round_arb_mat(acb_mat_t dest, const arb_mat_t src, slong prec)"

### Random generation
#ns arbcall"void acb_mat_randtest(acb_mat_t mat, flint_rand_t state, slong prec, slong mag_bits)"
#ns arbcall"void acb_mat_randtest_eig(acb_mat_t mat, flint_rand_t state, acb_srcptr E, slong prec)"

### Input and output
arbcall"void acb_mat_printd(const acb_mat_t mat, slong digits)"
#ns arbcall"void acb_mat_fprintd(FILE * file, const acb_mat_t mat, slong digits)"

### Comparisons
arbcall"int acb_mat_equal(const acb_mat_t mat1, const acb_mat_t mat2)"
arbcall"int acb_mat_overlaps(const acb_mat_t mat1, const acb_mat_t mat2)"
arbcall"int acb_mat_contains(const acb_mat_t mat1, const acb_mat_t mat2)"
#ni arbcall"int acb_mat_contains_fmpz_mat(const acb_mat_t mat1, const fmpz_mat_t mat2)"
#ni arbcall"int acb_mat_contains_fmpq_mat(const acb_mat_t mat1, const fmpq_mat_t mat2)"
arbcall"int acb_mat_eq(const acb_mat_t mat1, const acb_mat_t mat2)"
arbcall"int acb_mat_ne(const acb_mat_t mat1, const acb_mat_t mat2)"
arbcall"int acb_mat_is_real(const acb_mat_t mat)"
arbcall"int acb_mat_is_empty(const acb_mat_t mat)"
arbcall"int acb_mat_is_square(const acb_mat_t mat)"
arbcall"int acb_mat_is_exact(const acb_mat_t mat)"
arbcall"int acb_mat_is_zero(const acb_mat_t mat)"
arbcall"int acb_mat_is_finite(const acb_mat_t mat)"
arbcall"int acb_mat_is_triu(const acb_mat_t mat)"
arbcall"int acb_mat_is_tril(const acb_mat_t mat)"
arbcall"int acb_mat_is_diag(const acb_mat_t mat)"

### Special matrices
arbcall"void acb_mat_zero(acb_mat_t mat)"
arbcall"void acb_mat_one(acb_mat_t mat)"
arbcall"void acb_mat_ones(acb_mat_t mat)"
arbcall"void acb_mat_indeterminate(acb_mat_t mat)"
arbcall"void acb_mat_dft(acb_mat_t mat, int type, slong prec)"

### Transpose
arbcall"void acb_mat_transpose(acb_mat_t dest, const acb_mat_t src)"
arbcall"void acb_mat_conjugate_transpose(acb_mat_t dest, const acb_mat_t src)"
arbcall"void acb_mat_conjugate(acb_mat_t dest, const acb_mat_t src)"

### Norms
arbcall"void acb_mat_bound_inf_norm(mag_t b, const acb_mat_t A)"
arbcall"void acb_mat_frobenius_norm(acb_t res, const acb_mat_t A, slong prec)"
arbcall"void acb_mat_bound_frobenius_norm(mag_t res, const acb_mat_t A)"

### Arithmetic
arbcall"void acb_mat_neg(acb_mat_t dest, const acb_mat_t src)"
arbcall"void acb_mat_add(acb_mat_t res, const acb_mat_t mat1, const acb_mat_t mat2, slong prec)"
arbcall"void acb_mat_sub(acb_mat_t res, const acb_mat_t mat1, const acb_mat_t mat2, slong prec)"
arbcall"void acb_mat_mul_classical(acb_mat_t res, const acb_mat_t mat1, const acb_mat_t mat2, slong prec)"
arbcall"void acb_mat_mul_threaded(acb_mat_t res, const acb_mat_t mat1, const acb_mat_t mat2, slong prec)"
arbcall"void acb_mat_mul_reorder(acb_mat_t res, const acb_mat_t mat1, const acb_mat_t mat2, slong prec)"
arbcall"void acb_mat_mul(acb_mat_t res, const acb_mat_t mat1, const acb_mat_t mat2, slong prec)"
arbcall"void acb_mat_mul_entrywise(acb_mat_t res, const acb_mat_t mat1, const acb_mat_t mat2, slong prec)"
arbcall"void acb_mat_sqr_classical(acb_mat_t res, const acb_mat_t mat, slong prec)"
arbcall"void acb_mat_sqr(acb_mat_t res, const acb_mat_t mat, slong prec)"
arbcall"void acb_mat_pow_ui(acb_mat_t res, const acb_mat_t mat, ulong exp, slong prec)"
arbcall"void acb_mat_approx_mul(acb_mat_t res, const acb_mat_t mat1, const acb_mat_t mat2, slong prec)"

### Scalar arithmetic
arbcall"void acb_mat_scalar_mul_2exp_si(acb_mat_t B, const acb_mat_t A, slong c)"
arbcall"void acb_mat_scalar_addmul_si(acb_mat_t B, const acb_mat_t A, slong c, slong prec)"
#ni arbcall"void acb_mat_scalar_addmul_fmpz(acb_mat_t B, const acb_mat_t A, const fmpz_t c, slong prec)"
arbcall"void acb_mat_scalar_addmul_arb(acb_mat_t B, const acb_mat_t A, const arb_t c, slong prec)"
arbcall"void acb_mat_scalar_addmul_acb(acb_mat_t B, const acb_mat_t A, const acb_t c, slong prec)"
arbcall"void acb_mat_scalar_mul_si(acb_mat_t B, const acb_mat_t A, slong c, slong prec)"
#ni arbcall"void acb_mat_scalar_mul_fmpz(acb_mat_t B, const acb_mat_t A, const fmpz_t c, slong prec)"
arbcall"void acb_mat_scalar_mul_arb(acb_mat_t B, const acb_mat_t A, const arb_t c, slong prec)"
arbcall"void acb_mat_scalar_mul_acb(acb_mat_t B, const acb_mat_t A, const acb_t c, slong prec)"
arbcall"void acb_mat_scalar_div_si(acb_mat_t B, const acb_mat_t A, slong c, slong prec)"
#ni arbcall"void acb_mat_scalar_div_fmpz(acb_mat_t B, const acb_mat_t A, const fmpz_t c, slong prec)"
arbcall"void acb_mat_scalar_div_arb(acb_mat_t B, const acb_mat_t A, const arb_t c, slong prec)"
arbcall"void acb_mat_scalar_div_acb(acb_mat_t B, const acb_mat_t A, const acb_t c, slong prec)"

### Gaussian elimination and solving
arbcall"int acb_mat_lu_classical(slong * perm, acb_mat_t LU, const acb_mat_t A, slong prec)"
arbcall"int acb_mat_lu_recursive(slong * perm, acb_mat_t LU, const acb_mat_t A, slong prec)"
arbcall"int acb_mat_lu(slong * perm, acb_mat_t LU, const acb_mat_t A, slong prec)"
arbcall"void acb_mat_solve_tril_classical(acb_mat_t X, const acb_mat_t L, const acb_mat_t B, int unit, slong prec)"
arbcall"void acb_mat_solve_tril_recursive(acb_mat_t X, const acb_mat_t L, const acb_mat_t B, int unit, slong prec)"
arbcall"void acb_mat_solve_tril(acb_mat_t X, const acb_mat_t L, const acb_mat_t B, int unit, slong prec)"
arbcall"void acb_mat_solve_triu_classical(acb_mat_t X, const acb_mat_t U, const acb_mat_t B, int unit, slong prec)"
arbcall"void acb_mat_solve_triu_recursive(acb_mat_t X, const acb_mat_t U, const acb_mat_t B, int unit, slong prec)"
arbcall"void acb_mat_solve_triu(acb_mat_t X, const acb_mat_t U, const acb_mat_t B, int unit, slong prec)"
arbcall"void acb_mat_solve_lu_precomp(acb_mat_t X, const slong * perm, const acb_mat_t LU, const acb_mat_t B, slong prec)"
arbcall"int acb_mat_solve(acb_mat_t X, const acb_mat_t A, const acb_mat_t B, slong prec)"
arbcall"int acb_mat_solve_lu(acb_mat_t X, const acb_mat_t A, const acb_mat_t B, slong prec)"
arbcall"int acb_mat_solve_precond(acb_mat_t X, const acb_mat_t A, const acb_mat_t B, slong prec)"
arbcall"int acb_mat_inv(acb_mat_t X, const acb_mat_t A, slong prec)"
arbcall"void acb_mat_det_lu(acb_t det, const acb_mat_t A, slong prec)"
arbcall"void acb_mat_det_precond(acb_t det, const acb_mat_t A, slong prec)"
arbcall"void acb_mat_det(acb_t det, const acb_mat_t A, slong prec)"
arbcall"void acb_mat_approx_solve_triu(acb_mat_t X, const acb_mat_t U, const acb_mat_t B, int unit, slong prec)"
arbcall"void acb_mat_approx_solve_tril(acb_mat_t X, const acb_mat_t L, const acb_mat_t B, int unit, slong prec)"
arbcall"int acb_mat_approx_lu(slong * P, acb_mat_t LU, const acb_mat_t A, slong prec)"
arbcall"void acb_mat_approx_solve_lu_precomp(acb_mat_t X, const slong * perm, const acb_mat_t A, const acb_mat_t B, slong prec)"
arbcall"int acb_mat_approx_solve(acb_mat_t X, const acb_mat_t A, const acb_mat_t B, slong prec)"
arbcall"int acb_mat_approx_inv(acb_mat_t X, const acb_mat_t A, slong prec)"

### Characteristic polynomial and companion matrix
arbcall"void _acb_mat_charpoly(acb_ptr poly, const acb_mat_t mat, slong prec)"
arbcall"void acb_mat_charpoly(acb_poly_t poly, const acb_mat_t mat, slong prec)"
arbcall"void _acb_mat_companion(acb_mat_t mat, acb_srcptr poly, slong prec)"
arbcall"void acb_mat_companion(acb_mat_t mat, const acb_poly_t poly, slong prec)"

### Special functions
arbcall"void acb_mat_exp_taylor_sum(acb_mat_t S, const acb_mat_t A, slong N, slong prec)"
arbcall"void acb_mat_exp(acb_mat_t B, const acb_mat_t A, slong prec)"
arbcall"void acb_mat_trace(acb_t trace, const acb_mat_t mat, slong prec)"
arbcall"void _acb_mat_diag_prod(acb_t res, const acb_mat_t mat, slong a, slong b, slong prec)"
arbcall"void acb_mat_diag_prod(acb_t res, const acb_mat_t mat, slong prec)"

### Component and error operations
arbcall"void acb_mat_get_mid(acb_mat_t B, const acb_mat_t A)"
arbcall"void acb_mat_add_error_mag(acb_mat_t mat, const mag_t err)"

### Eigenvalues and eigenvectors
arbcall"int acb_mat_approx_eig_qr(acb_ptr E, acb_mat_t L, acb_mat_t R, const acb_mat_t A, const mag_t tol, slong maxiter, slong prec)"
arbcall"void acb_mat_eig_global_enclosure(mag_t eps, const acb_mat_t A, acb_srcptr E, const acb_mat_t R, slong prec)"
arbcall"void acb_mat_eig_enclosure_rump(acb_t lambda, acb_mat_t J, acb_mat_t R, const acb_mat_t A, const acb_t lambda_approx, const acb_mat_t R_approx, slong prec)"
arbcall"int acb_mat_eig_simple_rump(acb_ptr E, acb_mat_t L, acb_mat_t R, const acb_mat_t A, acb_srcptr E_approx, const acb_mat_t R_approx, slong prec)"
arbcall"int acb_mat_eig_simple_vdhoeven_mourrain(acb_ptr E, acb_mat_t L, acb_mat_t R, const acb_mat_t A, acb_srcptr E_approx, const acb_mat_t R_approx, slong prec)"
arbcall"int acb_mat_eig_simple(acb_ptr E, acb_mat_t L, acb_mat_t R, const acb_mat_t A, acb_srcptr E_approx, const acb_mat_t R_approx, slong prec)"
arbcall"int acb_mat_eig_multiple_rump(acb_ptr E, const acb_mat_t A, acb_srcptr E_approx, const acb_mat_t R_approx, slong prec)"
arbcall"int acb_mat_eig_multiple(acb_ptr E, const acb_mat_t A, acb_srcptr E_approx, const acb_mat_t R_approx, slong prec)"


text:
function sesam3P(;name, sN = sesam3N(name=:sN))
    @parameters t
    D = Differential(t)

    sts = @variables (begin
        L_P(t),  R_P(t),  I_P(t),
        dL_P(t),  dR_P(t),  dI_P(t),
        Φ_P(t), Φ_Pu(t), Φ_PB(t), Φ_Ptvr(t),
        u_PlantP(t), u_POM(t), u_immPPot(t),
        u_PPot(t), P_synBP(t), M_ImbP(t),
        β_PL(t), β_PR(t),
        leach_P(t),
        resorp_P(t), β_PBtvr(t),
        α_LP(t),
        α_RP(t),
        dec_LP(t), dec_RP(t), dec_LPPlant(t), dec_RPPlant(t),
        # need to be specified by coupled system:
        β_Pi(t), i_IP(t),
        u_PlantPmax(t), k_PlantP(t),
        s_EP(t), pL_sEP(t), # synthesis of E_LP adn E_RP enzymes by plants
        SOP(t), β_PSOM(t), β_NPSOM(t), β_NPB(t), β_NPi(t), β_NPL(t), β_NPR(t)
    end)
    ps = @parameters β_PEnz β_PB l_P  ν_P i_BN i_BP k_LP k_RP ρ_PBtvr=0.0

    @unpack L, R, B, SOC, SON, dec_L, dec_R, i_L, ϵ_tvr, tvr_B, syn_B, tvr_B0, syn_Enz, tvr_Enz, r_tvr, κ_E = sN
    @unpack β_NB, β_Ni, L_N, R_N = sN
    @unpack k_mN_L, k_mN_R = sN

    eqs = [
        β_PL ~ L/L_P, β_PR ~ R/R_P,
        D(L_P) ~ dL_P, dL_P ~ -dec_L/β_PL - dec_LP + i_L/β_Pi,
        D(R_P) ~ dR_P, dR_P ~ -dec_R/β_PR - dec_RP + ϵ_tvr*tvr_B/β_PBtvr + (1-κ_E)*tvr_Enz/β_PEnz,
        D(I_P) ~ dI_P,
        # assume that proportion of s_EP (plant P enzyme production flux) to biomineralizing L is
        # proportional to potential biomineralization flux
        pL_sEP ~ k_LP*L_P/(k_LP*L_P + k_RP*R_P),
        dec_LP ~ k_LP * L_P * (α_LP*syn_Enz + pL_sEP*s_EP)/(k_mN_L + α_LP*syn_Enz + pL_sEP*s_EP),
        dec_LPPlant ~ k_LP * L_P * (pL_sEP*s_EP)/(k_mN_L + pL_sEP*s_EP),
        # while_RP is decomp due to the sum of enzmyes of mic and plant (assumed to be the same with same k_M)
        # dec_RPPlant is the decomp that would be by enzymes levels by solely plant production, used in revenue calculation
        dec_RP ~ k_RP * R_P * (α_RP*syn_Enz + (1-pL_sEP)*s_EP)/(k_mN_R + α_RP*syn_Enz + (1-pL_sEP)*s_EP),
        dec_RPPlant ~ k_RP * R_P * ((1-pL_sEP)*s_EP)/(k_mN_R + (1-pL_sEP)*s_EP),
        u_PlantP ~ min(u_PlantPmax, k_PlantP*I_P),
        dI_P ~ i_IP - u_PlantP - leach_P  + dec_RP + dec_LP + Φ_P,
        leach_P ~ l_P*I_P,
        Φ_P ~ Φ_Pu + Φ_PB + Φ_Ptvr,
        Φ_Ptvr ~ r_tvr/β_PBtvr,
        Φ_Pu ~ (1-ν_P) * u_POM,
        # dec_RP is biomineralization - cleaved P is all mineralized
        # u_POM is from depolymerizing - only a fraction ν_P is mineralized
        u_PPot ~ ν_P * u_POM + u_immPPot + resorp_P,
        u_POM ~ dec_L/β_PL + dec_R/β_PR + κ_E*tvr_Enz/β_PEnz,
        u_immPPot ~ i_BP * I_P,
        P_synBP ~ u_PPot - syn_Enz/β_PEnz,
        M_ImbP ~ u_PPot - (syn_B/β_PB + syn_Enz/β_PEnz),
        Φ_PB ~ M_ImbP - u_immPPot,
        # make sure ρ_PBtvr >= 0
        resorp_P ~ ρ_PBtvr * tvr_B0/β_PB,
        β_PBtvr ~ tvr_B / (tvr_B0/β_PB - resorp_P),
        # observables for diagnostic output
        SOP ~ L_P + R_P + B/β_PB,
        β_PSOM ~ SOC/SOP,
        β_NPSOM ~ SON/SOP,
        β_NPB ~ β_PB/β_NB, # N:P of microbial biomass
        β_NPi ~ β_Pi/β_Ni, # N:P of litter inputs
        β_NPL ~ L_N/L_P,
        β_NPR ~ R_N/R_P,
        ]
    extend(ODESystem(eqs, t, sts, ps; name), sN)
end

function get_revenue_eq_sesam3CNP(sP)
    @parameters t
    @unpack α_L, α_R, dec_L, dec_R, β_NL, β_NR, β_NEnz, syn_Enz = sP
    @unpack β_PL, β_PR, β_PEnz = sP
    @unpack α_LP, α_RP, dec_LP, dec_RP, dec_RPPlant, dec_LPPlant = sP
    @unpack u_immPPot, u_PlantP, u_immNPot, u_PlantN, ν_N, ν_P = sP
    sts = @variables (begin
        α_LT(t), α_RT(t),
        α_LPT(t), α_RPT(t),
        syn_Enz_w(t),
        p_uPmic(t), p_uNmic(t), p_oPmic(t), p_oNmic(t),
        return_L(t), return_R(t), revenue_L(t), revenue_R(t),
        return_LP(t), return_RP(t), revenue_LP(t), revenue_RP(t),
        #invest_Ln(t), invest_Rn(t), return_Ln(t), return_Rn(t),
        revenue_sum(t)
    end)
    # need to be defined in coupled component:
    @variables w_C(t), w_N(t), w_P(t)
    eqs = [
        syn_Enz_w ~ syn_Enz*(w_C + w_N/β_NEnz + w_P/β_PEnz),
        # return_L ~ dec_L * (w_C + w_N/β_NL + w_P/β_PL),
        # return_R ~ dec_R * (w_C + w_N/β_NR + w_P/β_PR),
        return_L ~ dec_L * (w_C + w_N/β_NL*p_oNmic + w_P/β_PL*p_oPmic),
        return_R ~ dec_R * (w_C + w_N/β_NR*p_oNmic + w_P/β_PR*p_oPmic),
        # return of enzymes produced in addition to that of plants
        # only proportion of the biomineralization flux ends up in microbes: part it taken up by plant
        # dec_LP is already in P units, no need to devide by β_P
        p_uPmic ~ u_immPPot/(u_immPPot + u_PlantP),
        p_uNmic ~ u_immNPot/(u_immNPot + u_PlantN),
        p_oPmic ~ ν_P+(1-ν_P)*p_uPmic,
        p_oNmic ~ ν_N+(1-ν_N)*p_uNmic,
        return_LP ~ (dec_LP - dec_LPPlant) * p_uPmic * w_P,
        return_RP ~ (dec_RP - dec_RPPlant) * p_uPmic * w_P,
        revenue_L ~ return_L / (α_L * syn_Enz_w),
        revenue_R ~ return_R / (α_R * syn_Enz_w),
        # avoid deviding by zero. Checking > 0 fails with Derivative need to check > small_number
        revenue_LP ~ IfElse.ifelse(α_LP * syn_Enz_w > 1e-16, return_LP / (α_LP * syn_Enz_w), 0.0),
        revenue_RP ~ IfElse.ifelse(α_RP * syn_Enz_w > 1e-16, return_RP / (α_RP * syn_Enz_w), 0.0),
        revenue_sum ~ revenue_L + revenue_R + revenue_LP + revenue_RP,
        α_LT ~ revenue_L/revenue_sum,
        α_RT ~ revenue_R/revenue_sum,
        α_LPT ~ revenue_LP/revenue_sum,
        α_RPT ~ revenue_RP/revenue_sum,
        # auxilary for plotting
        # invest_Ln ~ invest_L/(invest_L + invest_R),
        # invest_Rn ~ invest_R/(invest_L + invest_R),
        # return_Ln ~ return_L/(return_L + return_R),
        # return_Rn ~ return_R/(return_L + return_R),
        ]
    (;eqs, sts)
end

function sesam3CNP(;name, δ=40.0, max_w=12, use_seam_revenue=false, sP=sesam3P(name=:sP))
    @parameters t
    D = Differential(t)
    @unpack α_L, α_R, syn_B, B, C_synBC, β_NB, N_synBN, tvr_B, τ = sP
    @unpack β_PB, P_synBP, α_LP, α_RP = sP
    sts = @variables (begin
        C_synBmC(t),
        C_synBN(t), C_synBmN(t),
        C_synBP(t), C_synBmP(t),
        #dα_L(t),
        dα_R(t),
        dα_LP(t), dα_RP(t),
        w_C(t), w_N(t), lim_C(t), lim_N(t),
        sum_w(t),
        w_P(t), lim_P(t)
    end)
    ps = @parameters δ=δ max_w=max_w
    eqs_rev, sts_rev = get_revenue_eq_sesam3CNP(sP)
    @variables α_LT(t) α_RT(t)
    @variables α_LPT(t) α_RPT(t)
    eqs = [
        C_synBN ~ β_NB*N_synBN,
        C_synBP ~ β_PB*P_synBP,
        syn_B ~ min(C_synBC, C_synBN, C_synBP),
        # C_synBmC ~ min(C_synBN, C_synBP),
        # C_synBmN ~ min(C_synBC, C_synBP),
        # C_synBmP ~ min(C_synBC, C_synBN),
        # need minimum, otherwise danger of Inf and NaN -> instability
        # w_C ~ min(max_w, exp(δ/tvr_B*(C_synBmC - C_synBC))),
        # w_N ~ min(max_w, exp(δ/tvr_B*(C_synBmN - C_synBN))),
        # w_P ~ min(max_w, exp(δ/tvr_B*(C_synBmP - C_synBP))),
        # need minimum inisde exp, otherwise computation of Duals may cause instability
        # w_C ~ exp(min(max_w, δ/tvr_B*(C_synBmC - C_synBC))),
        # w_N ~ exp(min(max_w, δ/tvr_B*(C_synBmN - C_synBN))),
        # w_P ~ exp(min(max_w, δ/tvr_B*(C_synBmP - C_synBP))),
        # with new formulation available flux cannot be too much smaller than syn_B
        # w_C ~ exp(-δ/tvr_B*(C_synBC - syn_B)),
        # w_N ~ exp(-δ/tvr_B*(C_synBN - syn_B)),
        # w_P ~ exp(-δ/tvr_B*(C_synBP - syn_B)),
        w_C ~ exp(min(max_w, -δ/tvr_B*(C_synBC - syn_B))),
        w_N ~ exp(min(max_w, -δ/tvr_B*(C_synBN - syn_B))),
        w_P ~ exp(min(max_w, -δ/tvr_B*(C_synBP - syn_B))),
        sum_w ~ w_C + w_N + w_P,
        lim_C ~ w_C/sum_w, lim_N ~ w_N/sum_w, lim_P ~ w_P/sum_w,# normalized for plot
        # α_LT, α_RT by get_revenue_eq_X
        #D(α_L) ~ dα_L, dα_L ~ (α_LT - α_L)*(τ + abs(syn_B)/B),
        α_L ~ 1 - (α_R + α_LP + α_RP),
        D(α_R) ~ dα_R, dα_R ~ (α_RT - α_R)*(τ + abs(syn_B)/B),
        D(α_LP) ~ dα_LP, dα_LP ~ (α_LPT - α_LP)*(τ + abs(syn_B)/B),
        D(α_RP) ~ dα_RP, dα_RP ~ (α_RPT - α_RP)*(τ + abs(syn_B)/B),
        ]
    extend(ODESystem(vcat(eqs,eqs_rev), t, vcat(sts, sts_rev), ps; name), sP)
end

sesam3(args...;kwargs...) = sesam3CNP(args...;kwargs...)

"""
R pool is a mixture of microbial turnover and enzymes.
Here the stable C/N ratio is computed based on given parameterization.
"""
function calculate_β_PR_sesam3(p,s)
    w_REnz = p[s.a_E]*(1-p[s.κ_E]) # a_E (1-κ_E) B
    w_RBtvr = p[s.ϵ_tvr]*p[s.τ]       # τ ϵ_tvr B
    # if ρ is not in parameter dictionary, its default is zero
    ρ_CBtvr = haskey(p, s.ρ_CBtvr) ? p[s.ρ_CBtvr] : 0.0
    ρ_PBtvr = haskey(p, s.ρ_PBtvr) ? p[s.ρ_PBtvr] : 0.0
    β_PBtvr = p[s.β_PB] * (1-ρ_CBtvr)/(1-ρ_PBtvr)
    β_PR = 1/((w_REnz/p[s.β_PEnz] + w_RBtvr/β_PBtvr)/(w_REnz + w_RBtvr))
end

ext: jl


text:
#-------------verlet.jl--------------------------------------------------------#
#
# Purpose: To implement the verlet force integrator
#
#------------------------------------------------------------------------------#

using Plots

# For this force integrator, we need the particle's current position and it's
#     previous position. We'll assume that it can only move in 1D for now.
type Particle
    pos::Float64
    prev_pos::Float64
    acc::Float64
end

function drop_ball(height::Float64, timestep::Float64)
    time = 0.0

    ball = Particle(height, height, -10.0)
    while ball.pos > 0.0
        curr_pos = ball.pos
        ball.pos = 2*ball.pos - ball.prev_pos +ball.acc*timestep*timestep
        ball.prev_pos = curr_pos
        time += timestep
        println(ball.pos)
    end

    println(time)
end

drop_ball(5.0, 0.1)

max_stars_repo_name: vchuravy/skillpill-julia

max_issues_repo_name: oist/skillpill-julia

max_forks_repo_name: oist/skillpill-julia


text:
# This file was generated, do not modify it.

x = exp(im*π) + 1

x ≈ 0

abs(x) < eps()

max_stars_repo_name: mcognetta/franklindocs

max_issues_repo_name: mcognetta/franklindocs

max_forks_repo_name: mcognetta/franklindocs


text:
using Distributions
using Test, Random

@test rand(Bernoulli()) isa Bool
@test rand(Bernoulli(), 10) isa Vector{Bool}


text:
using PyPlot
using Peacock
using LinearAlgebra
using SparseArrays
using BenchmarkTools
using CUDA

get_k = Peacock.get_k
DiagonalMatrix = Peacock.DiagonalMatrix

# Permittivity
function epf(x,y)
    # equation of a circle with radius 0.2a
    if x^2+y^2 <= 0.2^2
        # dielectric inside the circle
        return 8.9
    else
        # air outside the circle
        return 1
    end
end

# Permeability is unity everywhere
function muf(x,y)
    return 1
end

a1 = [1, 0]  # first lattice vector
a2 = [0, 1]  # second lattice vector
d1 = 0.01  # resolution along first lattice vector
d2 = 0.01  # resolution along second lattice vector
geometry = Geometry(epf, muf, a1, a2, d1, d2)
fourier_space_cutoff = 17
solver = Solver(geometry, fourier_space_cutoff)
G = BrillouinZoneCoordinate(  0,   0, "Γ")
X = BrillouinZoneCoordinate(1/2,   0, "X")
M = BrillouinZoneCoordinate(1/2, 1/2, "M")
ks = [G,X,M,G]
ks = [typeof(x)==BrillouinZoneCoordinate ? get_k(x,solver.basis) : x for x in ks]

polarisation = TE

k = ks[1]
bands = 1:2

basis = solver.basis
epc, muc = solver.epc, solver.muc
Kx = DiagonalMatrix(basis.kxs) + k[1]*I
Ky = DiagonalMatrix(basis.kys) + k[2]*I
Kx_d = CuArray(Kx); Ky_d = CuArray(Ky); epc_d = CuArray(epc); muc_d = CuArray(muc);
if polarisation == TE
    LHS = Kx/epc*Kx + Ky/epc*Ky
    RHS = muc
    LHS_d = CuArray(LHS)
    RHS_d = CuArray(RHS)
    LHS_dd = Kx_d/epc_d*Kx_d + Ky_d/epc_d*Ky_d
    RHS_dd = muc_d
    label = L"H_z"
elseif polarisation == TM
    LHS = Kx/muc*Kx + Ky/muc*Ky
    RHS = epc
    label = L"E_z"
end
# Sometimes the generalised eigenvalue problem solver
# fails near Γ when the crystals are symmetric.
# In these cases, rewrite as a regular eigenvalue problem
@time freqs_squared, modes_data = eigen(RHS \ LHS);
@time freqs_squared, modes_data = eigen(Array(RHS_d \ LHS_d));


text:
module ReducedBasisMethods

    using HDF5
    using HDF5: H5DataStore
    using Particles
    using TypedTables

    include("utils.jl")

    include("parameter.jl")

    export Parameter, hassamples

    include("parametersampler.jl")

    export ParameterSampler, CartesianParameterSampler, sample

    include("parameterspace.jl")

    export ParameterSpace

    include("regression.jl")

    export get_regression_αβ

    include("trainingset.jl")

    export TrainingSet

    include("time_marching.jl")

    export IntegratorParameters, IntegratorCache, ReducedIntegratorCache
    export integrate_vp, reduced_integrate_vp

    include("eigen.jl")

    export sorteigen

    include("projection.jl")

    export get_Ψ

    include("deim.jl")

    export deim_get_Π

    include("h5routines.jl")

    export h5save, h5load, read_sampling_parameters

end


text:
srand(1)

context("walks.jl") do
  landscapes = [
    NK.NKLandscape(4, 1),
    NK.NKqLandscape(4, 1, 2),
    NK.NKpLandscape(4, 1, 0.90),
  ]

  for ls = landscapes
    # The ith locus is linked to the locus i places to its right, assuming a
    # periodic boundary.
    ls.links = [
      0b1000,
      0b0100,
      0b0010,
      0b0001
    ]

    # Get the min fitness genotype so we know we can improve.
    gmin = NK.Genotype(0b000, ls)
    for a::NK.AlleleString = 0b000:0b1111
      g = NK.Genotype(a, ls)
      if NK.fitness(g) < NK.fitness(gmin)
        gmin = g
      end
    end

    context("$(ls)") do
      context("Adaptive walks") do
        function test_adaptive_walk(walk_function, genotype)
          walk = walk_function(genotype)
          @fact walk.length --> greater_than_or_equal(0)
          for i = 2:walk.length
            @fact NK.fitness(walk.history_list[i]) --> greater_than(NK.fitness(genotype))
          end
        end

        context("Random adaptive walk should terminate and move uphill") do
          test_adaptive_walk(NK.random_adaptive_walk, gmin)
        end

        context("Greedy adaptive walk should terminate and move uphill") do
          test_adaptive_walk(NK.greedy_adaptive_walk, gmin)
        end

        context("Reluctant adaptive walk should terminate and move uphill") do
          test_adaptive_walk(NK.reluctant_adaptive_walk, gmin)
        end
      end

      context("Neutral walks") do
        function test_neutral_walk(walk_function, genotype)
          walk = walk_function(genotype)
          @fact walk.length --> greater_than_or_equal(0)
          f = NK.fitness(genotype)
          for i = 1:walk.length
            f0 = NK.fitness(walk.history_list[i])
            @fact f0 --> roughly(f)
          end
        end

        context("Neutral walk should terminate and move sideways") do
          test_neutral_walk(NK.neutral_walk, gmin)
        end
      end
    end
  end
end


text:
using LabelledGraphs

export
    AbstractGibbsNetwork,
    network_graph,
    vertex_map,
    local_energy,
    interaction_energy,
    connecting_tensor,
    normalize_probability,
    boundary_state,
    local_state_for_node,
    iteration_order,
    fuse_projectors,
    update_gauges!

# S: type of the vertex of network
# T: type of the vertex of underlying factor graph
abstract type AbstractGibbsNetwork{S, T} end

struct NotImplementedError{M} <: Exception
    m::M
    NotImplementedError(m::M) where {M} = new{M}(m)
end

not_implemented(name) = throw(NotImplementedError(name))

factor_graph(network::AbstractGibbsNetwork{S, T}) where {S, T} = network.factor_graph

network_graph(network::AbstractGibbsNetwork{S, T}) where {S, T} = network.network_graph

vertex_map(network::AbstractGibbsNetwork{S, T}) where {S, T} = network.vertex_map

boundary(network::AbstractGibbsNetwork{S, T}, node::S) where {S, T} = not_implemented("boundary")

node_index(network::AbstractGibbsNetwork{S, T}, node::S) where {S, T} = not_implemented("node_index")

update_energy(network::AbstractGibbsNetwork{S, T}, σ::Vector{Int}) where {S, T} = not_implemented("update_energy")

conditional_probability(network::AbstractGibbsNetwork{S, T}, v::Vector{Int}) where {S, T} = not_implemented("conditional_probability")

iteration_order(peps::AbstractGibbsNetwork) = [(i, j) for i ∈ 1:peps.nrows for j ∈ 1:peps.ncols]

function projector(network::AbstractGibbsNetwork{S, T}, v::S, w::S) where {S, T}
    fg = factor_graph(network)
    vmap = vertex_map(network)
    fg_v, fg_w = vmap(v), vmap(w)

    if has_edge(fg, fg_w, fg_v)
        get_prop(fg, fg_w, fg_v, :pr)'
    elseif has_edge(fg, fg_v, fg_w)
        get_prop(fg, fg_v, fg_w, :pl)
    else
        loc_dim = fg_v ∈ vertices(fg) ? length(local_energy(network, v)) : 1
        ones(loc_dim, 1)
    end
end

function projector(network::AbstractGibbsNetwork{S, T}, v::S, W::NTuple{N, S}) where {S, T, N}
    first(fuse_projectors(
        [projector(network, v, w) for w ∈ W]
    ))
end

function fuse_projectors(projectors::Union{Vector{T}, NTuple{N, T}}) where {N, T}
    fused, energy = rank_reveal(hcat(projectors...), :PE)
    i0 = 1
    transitions = []
    for proj ∈ projectors
        ie = i0 + size(proj, 2) - 1
        push!(transitions, energy[:, i0:ie])
        i0 = ie + 1
    end
    fused, transitions
end

function spectrum(network::AbstractGibbsNetwork{S, T}, v::S) where {S, T}
    get_prop(factor_graph(network), vertex_map(network)(v), :spectrum)
end

function local_energy(network::AbstractGibbsNetwork{S, T}, v::S) where {S, T}
    spectrum(network, v).energies
end

function interaction_energy(network::AbstractGibbsNetwork{S, T}, v::S, w::S) where {S, T}
    fg = factor_graph(network)
    vmap = vertex_map(network)
    fg_v, fg_w = vmap(v), vmap(w)
    if has_edge(fg, fg_w, fg_v)
        get_prop(fg, fg_w, fg_v, :en)'
    elseif has_edge(fg, fg_v, fg_w)
        get_prop(fg, fg_v, fg_w, :en)
    else
        zeros(1, 1)
    end
end
#TODO: W T F
function _boundary_index(
    network::AbstractGibbsNetwork{S, T},
    v::S,
    w::Union{S, NTuple{N, S}},
    σ::Vector{Int}
) where {S, T, N}
    state = local_state_for_node(network, σ, v)
    if v ∉ vertices(network.network_graph) return one(typeof(state)) end
    pv = projector(network, v, w)
    [findfirst(x -> x > 0, pv[i, :]) for i ∈ 1:size(pv)[1]][state]
end

#TODO: W T F ^ 2
function _boundary_index(
    network::AbstractGibbsNetwork{S, T},
    v::S, w::S, k::S, l::S,
    σ::Vector{Int}
) where {S, T}
    pv = projector(network, v, w)
    i = _boundary_index(network, v, w, σ)
    j = _boundary_index(network, k, l, σ)
    (j - 1) * size(pv, 2) + i
end

function boundary_state(
    network::AbstractGibbsNetwork{S, T}, σ::Vector{Int}, node::S
) where {S, T}
     [
        _boundary_index(network, x..., σ)
        for x ∈ boundary(network, node)
    ]
end

function local_state_for_node(
    network::AbstractGibbsNetwork{S, T}, σ::Vector{Int}, w::S
) where {S, T}
    k = node_index(network, w)
    0 < k <= length(σ) ? σ[k] : 1
end

function is_compatible(
    factor_graph::LabelledGraph{T, NTuple{2, Int}},
    network_graph::LabelledGraph{S, NTuple{2, Int}}
) where {T, S}
    all(
        has_edge(network_graph, src(edge), dst(edge))
        for edge ∈ edges(factor_graph)
    )
end

function update_gauges!(
    network::AbstractGibbsNetwork{S, T},
    type::Symbol=:rand
) where {S, T}
    @assert type ∈ (:id, :rand)
    for i ∈ 1:network.nrows-1, k ∈ 1:1//2:network.ncols
        j = denominator(k) == 1 ? numerator(k) : k
        _, u, _, d = tensor_size(network, (i+1//2, j))
        Y = type == :id ? ones(u) : rand(u) .+ 0.1 #TODO: maybe 0.42?
        push!(network.gauges, (i + 1//6, j) => Y, (i + 2//6, j) => 1 ./ Y)
        Z = type == :id ? ones(d) : rand(d) .+ 0.1 #TODO: maybe 0.42?
        push!(network.gauges, (i + 4//6, j) => Z, (i + 5//6, j) => 1 ./ Z)
    end
end

function normalize_probability(values::Vector{R}) where {R<:Number}
    # exceptions (negative pdo, etc)
    minp = minimum(values)
    if minp < 0
        amp = abs(minp)
        for (i, p) ∈ enumerate(values)
            p < amp ? values[i] = amp : values[i]
        end
    end
    values / sum(values)
end


text:
abstract type AbstractWebPage end
abstract type AbstractWebPageAspect{T} <: AbstractWebPage{T} end
abstract type AbstractWebPageRole{T}   <: AbstractWebPage{T} end

struct WebPageAspect{T} <: AbstractWebPageAspect{T}
    tag::Symbol
end
struct WebPageRole{T}   <: AbstractWebPageRole{T}
    tag::Symbol
end

const WebPageAspects = Set([])
const WebPageRoles   = Set([])

@quasiabstract mutable struct HtmlShared <: AbstractWebPage
    tag::Symbol           # tag  === Symbol(name)
    name::String          # name === string(tag)
end

tag(x::HtmlShared)  = x.tag
name(x::HtmlShared) = x.name

function HtmlShared(tag::Symbol)
    name = string(tag)
    return HtmlShared(tag, name)
end
function HtmlShared(name::String)
    tag = @sym(name)
    return HtmlShared{T}(tag, name)
end
function HtmlShared(;tag::Symbol, name::String)
    return HtmlShared(tag, name)
end

@quasiabstract mutable struct HtmlItemShared <: HtmlShared
    id::String   # an identifier for the realized item that is unique throughout website
end

id(x::HtmlItemShared) = x.id

function HtmlItemShared(id::String)
    return HtmlShared(tag, name)
end

# Forward methods from `HtmlItemShared` to `HtmlShared`
@forward (HtmlItemShared, :htmlshared) RegularPolygon
abstract type AbstractWebPageElement{Aspect, Role} <: AbstractWebPage{T} end

abstract type AbstractWebPageAttribute{T} <: AbstractWebPageElement{T} end
abstract type AbstractWebPageAttributeValue{T} <: AbstractWebPageElement{T} end

abstract type AbstractHtml{T} <: AbstractWebPage{T} end
abstract type AbstractHtmlElement{T, V} <: AbstractWebPageElement{T} end
abstract type AbstractHtmlAttribute{T, V} <: AbstractWebPageAttribute{T} end
abstract type AbstractHtmlAttributeValue{T, V} <: AbstractWebPageAttributeValue{T} end

abstract type AbstractXHtml{T} <: AbstractWebPage{T} end
abstract type AbstractXHtmlElement{T, V} <: AbstractWebPageElement{T} end
abstract type AbstractXHtmlAttribute{T, V} <: AbstractWebPageAttribute{T} end
abstract type AbstractXHtmlAttributeValue{T, V} <: AbstractWebPageAttributeValue{T} end

# AbstractWeb params

struct HtmlAttribute{T} end
struct HtmlAttributeValue{T} end

const  web_boolean_value = WebAttributeValue{Bool}()
const  web_integer_value = WebAttributeValue{Int}()
const  web_float_value = WebAttributeValue{Float64}()
const  web_string_value = WebAttributeValue{String}()

struct HtmlAttribute{V} <: AbstractHtmlAttribute{T, V}
    name::String
    valuation::V
end

abstract type AbstractWebPage{Kind} end
abstract type AbstractWebConstituent{Sort, Kind} end
abstract type AbstractWebElement{Sort}   <: AbstractWebConstituent{Sort, Kind} end
abstract type AbstractWebAttribute{Sort} <: AbstractWebConstituent{Sort, Kind} end
abstract type AbstractWebProperty{Sort}  <: AbstractWebConstituent{Sort, Kind} end
abstract type AbstractWebValuation{Sort} <: AbstractWebConstituent{Sort, Kind} end

# AbstractWebPage Kinds
struct HTML5 end
struct XHTML end

# AbstractWeb Value Constituents
struct BooleanValued end
struct EnumValued end
struct IntegerValued end
struct NonNegativeIntegerValued end
struct PositiveIntegerValued end
struct FloatValued end

abstract type AbstractHTML5{URI} <: AbstractWebPage{HTML5}
abstract type AbstractXHTML{URI} <: AbstractWebPage{XHTML}


text:
#Requires latest Sobol master.

function sobol_seq(::Val{N}) where N
    isa(N, Integer) || error("$N is not an integer.")
    (N < 0 || N > (length(Sobol.sobol_a) + 1)) && error("invalid Sobol dimension")

    m = ones(UInt32, (N, 32))

    #special cases
    N == 0 && return(SobolSeq{0})(m,UInt32[],UInt32[],zero(UInt32))
    #special cases 1
    N == 1 && return(SobolSeq{N}(m,UInt32[0],UInt32[0],zero(UInt32)))

    for i = 2:N
        a = Sobol.sobol_a[i-1]
        d = floor(Int, log2(a)) #degree of poly

        #set initial values of m from table
        m[i, 1:d] = Sobol.sobol_minit[1:d, i - 1]
        #fill in remaining values using recurrence
        for j = (d+1):32
            ac = a
            m[i,j] = m[i,j-d]
            for k = 0:d-1
                @inbounds m[i,j] = m[i,j] ⊻ (((ac & one(UInt32)) * m[i, j-d+k]) << (d-k))
                ac >>= 1
            end
        end
    end
    SobolSeq{N}(m,zeros(UInt32,N),zeros(UInt32,N),zero(UInt32))
end

function sobol_vec(::Val{N}, npoints::Int = 64, ::Type{T} = Float64) where {N,T}
    s = sobol_seq(Val{N}())
    x = Vector{T}(undef, N)
    Sobol.skip!(s, npoints, x)

    out = Vector{NTuple{N,T}}(undef, npoints)

    for i ∈ 1:npoints
        if s.n == typemax(s.n)
            return rand!(x)
        end

        s.n += one(s.n)
        c = UInt32(trailing_zeros(s.n))

        out[i] = ntuple(i -> gen_sobol_x!(i, c, s.b, s.x, s.m) , Val{N}())

    end
    out
end

function gen_sobol_x!(i, c, sb, sx, sm)
    @inbounds b = sb[i]
    if b >= c
        @inbounds sx[i] = sx[i] ⊻ (sm[i,c+1] << (b-c))
        @inbounds x = sx[i] * Sobol.scale2m[b+1]
    else
        @inbounds sx[i] = (sx[i] << (c-b)) ⊻ sm[i,c+1]
        @inbounds sb[i] = c
        @inbounds x = sx[i] * Sobol.scale2m[c+1]
    end
    #logit(x)
    √2 * erfinv(2x-1)
end

function sobol_vec_reinterp_ver(::Val{N}, npoints::Int = 64, ::Type{T} = Float64) where {N,T}
    s = sobol_seq(Val{N}())
    x = Vector{T}(undef, N)
    Sobol.skip!(s, npoints, x)

    out = Vector{NTuple{N,T}}(undef, npoints)
    out2 = reshape(reinterpret(T, out),N,npoints)

    for i ∈ 1:npoints
        next!(s, @view(out2[:,i]))
    end
    out
end

logit(x) = log(x / (1-x))


text:
@discord_object struct ActivityAssets
    large_image::String
    large_text::String
    small_image::String
    small_text::String
end


text:
"""
The `ViscousFlow` module is here

"""
module ViscousFlow

using Reexport

#== Imports/Exports ==#

include("Utils.jl")
@reexport using .Utils

include("Fields.jl")

@reexport using .Fields

include("Bodies.jl")

@reexport using .Bodies

include("RigidBodyMotions.jl")

@reexport using .RigidBodyMotions

include("Layers.jl")

@reexport using .Layers

include("SaddlePointSystems.jl")

@reexport using .SaddlePointSystems

include("TimeMarching.jl")
@reexport using .TimeMarching

include("Systems.jl")
@reexport using .Systems

#== Plot Recipes ==#

include("plot_recipes.jl")

end


text:
function FixedTransformProvider(arg0::Transform)
    return FixedTransformProvider((Transform,), arg0)
end


text:
module Dispatch
using ..PartialOrder
using ..PartialOrder: Node, insert!

export DNode, Decision, MethodCall, NoMethodNode
export domainof, signatureof, make_namer
export addmethod!, simplify, build_dtree

import ..Patterns
const INode = Patterns.Node # todo: remove!

## Method interface: ##

signatureof(method)  = error("unimplemented!")
domainof(method)     = error("unimplemented!")
is_empty_domain(dom) = error("unimplemented!")
hullof(method)       = error("unimplemented!")
make_namer(methods::Vector) = error("unimplemented!")

# ---- Method DAG manipulation ------------------------------------------------

methodsof{M}(top::Node{M}) = [node.value for node in subDAGof(top)]

function addmethod!{M}(top::Node{M}, name::Symbol, m::M)
    insert!(top, Node(m))

    methods = methodsof(top)
    for mk in methods
        if mk === m;  continue  end
        lb = domainof(m) & domainof(mk)
        if is_empty_domain(lb); continue; end
        if any([domainof(ml) == lb for ml in methods]) continue; end

        sig1 = signatureof(m,  "_A")
        sig2 = signatureof(mk, "_B")

        println("Warning: New @pattern method ", name, sig1)
        println("         is ambiguous with   ", name, sig2, '.')
        println("         Make sure ", name, sig1 & sig2," is defined first.")
    end
end

function simplify{M}(top::Node{M}, hull)
    top = copyDAG(top)

    # filter out too specific methods
    keep = Set{M}(filter(m->!(hullof(m) < hull), methodsof(top)))
    @assert !isempty(keep)
    raw_filter!(top, keep)

    # filter out non-questions
    top = simplify!(top, hull)
end

simplify!{T}(n::Node{T}, domain) = simplify!(Dict{Node{T},Node{T}}(), n, domain)
function simplify!{T}(subs::Dict{Node{T},Node{T}}, node::Node{T}, domain)
    if haskey(subs, node);  return subs[node]  end
    ndom = node.value & domain
    for child in node.gt
        if ndom == child.value & domain
            return subs[node] = simplify!(child, domain)
        end
    end

    node.gt = Set{Node{T}}([simplify!(child, domain) for child in node.gt])
    subs[node] = node
end

# ---- Decision Tree ----------------------------------------------------------

abstract DNode{M}

type Decision{M} <: DNode{M}
    domain
    pass::DNode
    fail::DNode
    methods::Vector{M}
    pre::Vector{INode}
    seq::Vector{INode}

    Decision(domain, pass, fail, ms) = new(domain, pass, fail, ms)
end
Decision{M}(dom, pass, fail, ms::Vector{M}) = Decision{M}(dom, pass, fail, ms)

type MethodCall{M} <: DNode{M}
    m::M
    bind_seq::Vector
    bindings::Vector{INode}

    MethodCall(m::M) = new(m)
end
MethodCall{M}(m::M) = MethodCall{M}(m)

type NoMethodNode <: DNode; end
const nomethodnode = NoMethodNode()

# ---- create decision tree ---------------------------------------------------

function choose_pivot{M}(top::Node{M}, ms::Set{Node{M}})
    nmethods = length(ms)
    p_opt = nothing
    n_opt = nmethods+1
    for pivot in intersect(top.gt, ms)
        below = intersect(ms, subDAGof(pivot))
        npass = length(below)
        nfail = nmethods - npass
        n = max(npass, nfail)
        if n < n_opt
            p_opt, n_opt = pivot, n
        end
    end
    p_opt::Node{M}
end

build_dtree(top::Node) = build_dtree(top, subDAGof(top))
function build_dtree{M}(top::Node{M}, ms::Set{Node{M}})
    if isempty(top.gt) || length(ms) == 1
        top.value.body === nothing ? nomethodnode : MethodCall(top.value)
    else
        pivot = choose_pivot(top, ms)
        below = subDAGof(pivot)
        pass = build_dtree(pivot, intersect(ms, below))
        fail = build_dtree(top,   setdiff(ms, below))

        methods = M[node.value for node in filter(node->(node in ms),
                                                  ordered_subDAGof(top))]
        Decision(domainof(pivot.value), pass, fail, methods)
    end
end

end # module


text:
# Autogenerated wrapper script for libsass_jll for i686-w64-mingw32-cxx11
export libsass_so

JLLWrappers.@generate_wrapper_header("libsass")
JLLWrappers.@declare_library_product(libsass_so, "libsass-1.dll")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libsass_so,
        "bin\\libsass-1.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
ENV["CUDA_VISIBLE_DEVICES"] = "0,1,2,3,4,5,6,7"
using Revise
using ADSeismic
using ADCME
using MAT
using PyPlot
using Random
using DelimitedFiles
using Optim
using Dates
matplotlib.use("Agg")
close("all")
if has_gpu()
  gpu = true
  # use_gpu()
else
  gpu = false
end

data_dir = "data/acoustic"
if !ispath(data_dir)
  mkpath(data_dir)
end
figure_dir = "figure/NNFWI/marmousi/"
if !ispath(figure_dir)
  mkpath(figure_dir)
end
result_dir = "result/NNFWI/marmousi/"
if !ispath(result_dir)
  mkpath(result_dir)
end
model_dir = "NN_model/NNFWI/marmousi/"
if !ispath(model_dir)
  mkpath(model_dir)
end
loss_file = joinpath(result_dir, "loss_$(Dates.now()).txt")

################### Inversion using Automatic Differentiation #####################
reset_default_graph()
model_name = "models/marmousi2-model-smooth.mat"
# model_name = "models/BP-model-smooth.mat"

## load model setting
params = load_params(model_name, vp_ref=3e3, PropagatorKernel=1)
# params.NSTEP = 2
src = load_acoustic_source(model_name)
rcv = load_acoustic_receiver(model_name)
vp0 = matread(model_name)["vp"]
mean_vp0 = mean(vp0)
std_vp0 = std(vp0)
vp0 = constant(vp0[:,:])

## original fwi
# vp_ = Variable(vp0 / mean_vp0)
# vp = vp_ * mean_vp0
# vars = vp_

## fwi with NN
# nx = params.NX
# ny = params.NY
# dx = params.DELTAX
# dy = params.DELTAY
# z = zeros((nx+2)*(ny+2),2)
# for i = 1:nx+2
#  for j = 1:ny+2
#     z[(i-1)*(ny+2)+j, :] = [(i-1)*dx (j-1)*dy]
#   end
# end
# vp = squeeze(fc(z, [10,1])) * mean_vp0
# vp = reshape(vp, size(vp0))
Random.seed!(0);
z = constant(rand(Float32, 1, 8))
x = Generator(z, ratio=size(vp0)[1]/size(vp0)[2], base=4)
x = tf.cast(x, tf.float64)
x = x * std_vp0
vp = add_initial_model(x, vp0)

## assemble acoustic propagator model
model = x->AcousticPropagatorSolver(params, x, vp)
vars = get_collection()

## load data
std_noise = 0
Random.seed!(1234);
Rs = Array{Array{Float64,2}}(undef, length(src))
for i = 1:length(src)
  Rs[i] = readdlm(joinpath(data_dir, "marmousi-r$i.txt"))
    # Rs[i] = readdlm(joinpath(data_dir, "BP-r$i.txt"))
    Rs[i] .+= randn(size(Rs[i])) .* std(Rs[i]) .* std_noise
end

## calculate loss
if gpu
  losses, grads = compute_loss_and_grads_GPU(model, src, rcv, Rs, vars)
  grad = sum(grads)
  loss = sum(losses)
else
  [SimulatedObservation!(model(src[i]), rcv[i]) for i = 1:length(src)]
  loss = sum([sum((rcv[i].rcvv-Rs[i])^2) for i = 1:length(rcv)])
end

global_step = tf.Variable(0, trainable=false)
max_iter = 50000
lr_decayed = tf.train.cosine_decay(0.001, global_step, max_iter)
opt = AdamOptimizer(lr_decayed).minimize(loss, global_step=global_step, colocate_gradients_with_ops=true)

sess = Session(); init(sess)
loss0 = run(sess, loss)
@info "Initial loss: ", loss0

## run inversion
fp = open(loss_file, "w")
write(fp, "0,$loss0\n")
function callback(vs, iter, loss)
  if iter%50==1
    x = run(sess, vp)'
    clf()
    pcolormesh([0:params.NX+1;]*params.DELTAX/1e3,[0:params.NY+1;]*params.DELTAY/1e3,  x)
    axis("scaled")
    colorbar(shrink=0.4)
    xlabel("x (km)")
    ylabel("z (km)")
    gca().invert_yaxis()
    title("Iteration = $iter")
    savefig(joinpath(figure_dir, "inv_$(lpad(iter,5,"0")).png"), bbox_inches="tight")
    writedlm(joinpath(result_dir, "inv_$(lpad(iter,5,"0")).txt"), x)
  end
  if iter%1000 == 1
    ADCME.save(sess, joinpath(model_dir, "NNFWI_$(lpad(iter,5,"0")).mat"))
  end
  write(fp, "$iter,$loss\n")
  flush(fp)
end

## optimization using Adam
time = 0
for iter = 1:max_iter
  global time += @elapsed  _, ls, lr = run(sess, [opt, loss, lr_decayed])
  callback(nothing, iter, ls)
  println("   $iter\t$ls\t$lr")
  println(" * time: $time")
end

## optimization using BFGS
# Optimize!(sess, loss, vars=vars, grads=grad,  callback=callback)
close(fp)


text:
using Turing

@model simplegaussmodel(x) = begin
  s ~ InverseGamma(2,3)
  m ~ Normal(0,sqrt(s))
  for i in 1:length(x)
    x[i] ~ Normal(m, sqrt(s))
  end
  s, m
end


text:
using Images, ImageView
using AprilTags
using FixedPointNumbers
# Package Video4Linux is currently unregistered, and should be installed with:
# Pkg.clone("https://github.com/Affie/Video4Linux.jl.git")
# Pkg.build("Video4Linux")
using Video4Linux
using Dates

function showImage!(image, tags, imageCol)
    # Convert image to RGB
    imageCol[:,:] = RGB{N0f8}.(view(img,:,:), view(img,:,:), view(img,:,:))
    # draw color box on tag corners
    foreach(tag->drawTagBox!(imageCol, tag, width=3, drawReticle=false), tags)
    nothing
end

#Use this colour mapping when using the PS3eye.
#NOTE on Video4Linux:
#No automatic video device setup is currently performed by Video4Linux
#There are utilities to make this easy eg: qv4l2 (install with: sudo apt-get install qv4l2)
#qv4l2 provides a gui of the ioctl of the device.
#Settings such as frame size and capture format can easily be changed.
yonly = Video4Linux.YUYVonlyY(640,480)
#create a video channel that produces image frames
vidchan = Channel((c::Channel) -> videoproducer(c, yonly, devicename = "/dev/video1",
                                         iomethod = Video4Linux.IO_METHOD_MMAP, N=0 ))
##
# capture one frame
A = take!(vidchan)
# keep pointers to the same memory, preallocate variables
img = normedview(A)
imC = RGB.(view(img,:,:), view(img,:,:), view(img,:,:))
#create a new canvas for displaying the image
canvas = imshow(imC)

# Create default detector
detector = AprilTagDetector()
# settings that influence detector speed/quality
# see april tag documentation for more information on these settings
detector.nThreads = 4
detector.quad_decimate =  1.0
detector.quad_sigma = 0.0
detector.refine_edges = 0
detector.decode_sharpening = 0.25

# start asyncronious task to read video channel, detect tags and display detections
taskvid = @async begin
    i = 0
    starttime = Dates.value(now())
    while isopen(vidchan)
        i += 1
        A = take!(vidchan)
        img[:,:] = normedview(A)[:,:]

        tags = detector(img[:,:])
        showImage!(img, tags, imC)
        ImageView.imshow!(canvas["gui"]["canvas"], imC, canvas["annotations"])

        if i % 10 == 0
            @show length(tags)
            framerate = round(10/(Dates.value(now())-starttime)*1000,digits=1)
            println("frames $(framerate) per second")
            starttime = Dates.value(now())
        end
    end
end

## to finish and close the video device
stopVideoProducer()


text:
@testset "Vlasov-Poisson Landau 1D1V" begin

dev = CPU()                  # device
nx, nv = 64, 64              # grid resolution
stepper = StrangSplitting()  # timestepper
dt = 0.01                    # timestep
nsteps = 10                  # total number of time-steps

xmin, xmax = 0, 4π           # X Domain length (m)
vmin, vmax = -6, 6           # V Domain length (m)
α  = 0.5                     # Perturbation amplitude
kx = 0.5                     # Wave number of perturbation

xgrid = OneDGrid(dev, nx, xmin, xmax)
vgrid = OneDGrid(dev, nv, vmin, vmax)

f = DistributionFunction( xgrid, vgrid )

landau!(f, α, kx)

prob = VlasovProblem(f, BSLSpline(5), dev)

sol = solve!(prob, stepper, dt, nsteps )

@test true

end


text:
module VectorizedRNG

using VectorizationBase, Random, UnPack
using VectorizationBase: simd_integer_register_size, gep, _Vec, ifelse, VecUnroll, AbstractSIMD,
    rotate_right, vadd, vsub, zero_offsets, vfmadd, vfmsub, vfnmadd, shufflevector, register_size,
    cache_linesize, StaticInt, pick_vector_width, data

using Distributed: myid

export local_rng, rand!, randn!#, randexp, randexp!

abstract type AbstractVRNG{N} <: Random.AbstractRNG end
abstract type AbstractState{N,W} end

@inline vloadu(ptr::Ptr) = VectorizationBase.__vload(ptr, VectorizationBase.False(), register_size())
@inline vloadu(ptr::Ptr, i) = VectorizationBase.__vload(ptr, i, VectorizationBase.False(), register_size())
@inline vloada(ptr::Ptr) = VectorizationBase.__vload(ptr, VectorizationBase.True(), register_size())
@inline vloada(ptr::Ptr, i) = VectorizationBase.__vload(ptr, i, VectorizationBase.True(), register_size())
@inline vloada(ptr::Ptr, i, m) = VectorizationBase.__vload(ptr, i, m, VectorizationBase.True(), register_size())
@inline vstorea!(ptr::Ptr, v) = VectorizationBase.__vstore!(ptr, v, VectorizationBase.True(), VectorizationBase.False(), VectorizationBase.False(), register_size())
@inline vstorea!(ptr::Ptr, v, i) = VectorizationBase.__vstore!(ptr, v, i, VectorizationBase.True(), VectorizationBase.False(), VectorizationBase.False(), register_size())
@inline vstorea!(ptr::Ptr, v, i, m) = VectorizationBase.__vstore!(ptr, v, i, m, VectorizationBase.True(), VectorizationBase.False(), VectorizationBase.False(), register_size())
@inline vstoreu!(ptr::Ptr, v, i) = VectorizationBase.__vstore!(ptr, v, i, VectorizationBase.False(), VectorizationBase.False(), VectorizationBase.False(), register_size())
@inline vstoreu!(ptr::Ptr, v) = VectorizationBase.__vstore!(ptr, v, VectorizationBase.False(), VectorizationBase.False(), VectorizationBase.False(), register_size())

include("masks.jl")
include("api.jl")
include("special_approximations.jl")
include("xoshiro.jl")
# const GLOBAL_vPCGs = Ref{Ptr{UInt64}}()

const GLOBAL_vRNGs = Ref{Ptr{UInt64}}()

const RNG_MEM_SIZE = (5(simd_integer_register_size()*XREGISTERS + 2048*3))
local_rng(i) = Xoshift{XREGISTERS}(i*(RNG_MEM_SIZE) + GLOBAL_vRNGs[])
local_rng() = local_rng(Base.Threads.threadid() - 1)

# include("precompile.jl")
# _precompile_()

function __init__()
    # ccall(:jl_generating_output, Cint, ()) == 1 && return
  nthreads = Base.Threads.nthreads()
  GLOBAL_vRNGs[] = ptr = VectorizationBase.valloc((RNG_MEM_SIZE ÷ 8)*nthreads, UInt64)
  nstreams = XREGISTERS * nthreads * simd_integer_register_size()
  initXoshift!(ptr, nstreams)
  for tid ∈ 0:nthreads-1
    rng = local_rng(tid)
    setrandu64counter!(rng, 0x00)
    # setrandn32counter!(rng, 0x00)
    setrand64counter!(rng, 0x00)
    setrandn64counter!(rng, 0x00)
  end
end

end # module


text:
"""
    mapmodel(f, model::M, [constructor])

Applies `f` to every layer, resulting in a new model.

!!! warning

    The model layer structs are copied in the process, but the
    fields are not, hence the resulting model will have the
    same parameters!

!!! warning

    The copying is done by first recursively applying `mapmodel` to
    all children fields and then attempting to create a new object
    of type `M` with the default outer constructor `M(fields...)`.

If your custom type does not have a constructor of that form,
like e.g. `Flux.Chain`, you will have to extend this function
with a signature like

    `mapmodel(f, model::MyType, (fields) -> ...)`

where the last argument results in a `MyType(fields...)`
"""
function mapmodel(f, model::M, constructor = M) where M
    fields = []
    for (name, T) in zip(fieldnames(M), fieldtypes(M))
        value::T = getfield(model, name)
        if attrtype(model, Val(name)) == ModelUtils.ChildAttr
            push!(fields, mapmodel(f, value))
        else
            push!(fields, value)
        end
    end
    return f(constructor(fields...))
end

mapmodel(f, t::T) where T<:Tuple = (x -> mapmodel(f, x)).(t)
mapmodel(f, g::Function) = g
mapmodel(f, model::C) where C<:Chain = mapmodel(f, model, (layers) -> Chain(layers...))
mapmodel(f, model::M) where M<:Model = Model(mapmodel(f, model.layer))


text:
function legend!(fg::FigureGrid; position=:right,
                 orientation=default_orientation(position), kwargs...)

    guide_pos = guides_position(fg.figure, position)
    return legend!(guide_pos, fg; orientation, kwargs...)
end

"""
    legend!(figpos, grid; kwargs...)

Compute legend for `grid` (which should be the output of [`draw!`](@ref)) and draw it in
position `figpos`. Attributes allowed in `kwargs` are the same as `MakieLayout.Legend`.
"""
function legend!(figpos, grid; kwargs...)
    legend = compute_legend(grid)
    return isnothing(legend) ? nothing : Legend(figpos, legend...; kwargs...)
end

"""
    plottypes_attributes(entries)

Return plottypes and relative attributes, as two vectors of the same length,
for the given `entries`.
"""
function plottypes_attributes(entries)
    plottypes = PlotFunc[]
    attributes = Vector{Symbol}[]
    for entry in entries
        # FIXME: this should probably use the rescaled values
        plottype = Makie.plottype(entry.plottype, entry.positional...)
        n = findfirst(==(plottype), plottypes)
        attrs = (keys(entry.primary)..., keys(entry.named)...)
        if isnothing(n)
            push!(plottypes, plottype)
            push!(attributes, collect(Symbol, attrs))
        else
            union!(attributes[n], attrs)
        end
    end
    return plottypes, attributes
end

hassymbolkey((k, v)::Pair) = k isa Symbol

compute_legend(fg::FigureGrid) = compute_legend(fg.grid)
function compute_legend(grid::Matrix{AxisEntries})
    # assume all subplots have same scales, to be changed to support free scales
    scales = filter(hassymbolkey, first(grid).scales)

    # remove keywords that don't support legends
    for key in [:row, :col, :layout, :stack, :dodge, :group]
        pop!(scales, key, nothing)
    end

    # if no legend-worthy keyword remains return nothing
    isempty(scales) && return nothing

    titles = unique!(collect(String, (scale.label for scale in values(scales))))
    # empty strings create difficulties with the layout
    nonemptytitles = map(t -> isempty(t) ? " " : t, titles)

    plottypes, attributes = plottypes_attributes(entries(grid))

    labels = Vector{String}[]
    elements_list = Vector{Vector{LegendElement}}[]

    for title in titles
        label_attrs = [key for (key, val) in scales if val.label == title]
        uniquevalues = mapreduce(k -> datavalues(scales[k]), assert_equal, label_attrs)
        elements = map(eachindex(uniquevalues)) do idx
            local elements = LegendElement[]
            for (P, attrs) in zip(plottypes, attributes)
                shared_attrs = attrs ∩ label_attrs
                isempty(shared_attrs) && continue
                options = [attr => plotvalues(scales[attr])[idx] for attr in shared_attrs]
                append!(elements, legend_elements(P; options...))
            end
            return elements
        end
        push!(labels, map(string, uniquevalues))
        push!(elements_list, elements)
    end
    return elements_list, labels, nonemptytitles
end

# Notes

# TODO: correctly handle composite plot types (now fall back to poly)
# TODO: make legend updateable?
# TODO: allow custom attributes in legend elements?


text:
# 各点の密度を求める
function calc_density(input_point, sim_umbrella_datas, sigma_matrix)
    density = 0.0
    for i = 1:size(sim_umbrella_datas, 1)
        density += calc_gaussian(input_point, sim_umbrella_datas[i, :], sigma_matrix)
    end

    density /= size(sim_umbrella_datas, 1)
    return density
end

# 多変量ガウス
function calc_gaussian(input, mean, sigma_matrix)
    return pdf(MvNormal(mean, sigma_matrix), input)
end


text:
# ---
# title: 75. Sort Colors
# id: problem75
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Array, Two Pointers, Sort
# link: <https://leetcode.com/problems/sort-colors/description/>
# hidden: true
# ---
#
# Given an array `nums` with `n` objects colored red, white, or blue, sort them
# **[in-place](https://en.wikipedia.org/wiki/In-place_algorithm) **so that
# objects of the same color are adjacent, with the colors in the order red,
# white, and blue.
#
# Here, we will use the integers `0`, `1`, and `2` to represent the color red,
# white, and blue respectively.
#
# **Follow up:**
#
#   * Could you solve this problem without using the library's sort function?
#   * Could you come up with a one-pass algorithm using only `O(1)` constant space?
#
#
#
# **Example 1:**
#
#
#
#     Input: nums = [2,0,2,1,1,0]
#     Output: [0,0,1,1,2,2]
#
#
# **Example 2:**
#
#
#
#     Input: nums = [2,0,1]
#     Output: [0,1,2]
#
#
# **Example 3:**
#
#
#
#     Input: nums = [0]
#     Output: [0]
#
#
# **Example 4:**
#
#
#
#     Input: nums = [1]
#     Output: [1]
#
#
#
#
# **Constraints:**
#
#   * `n == nums.length`
#   * `1 <= n <= 300`
#   * `nums[i]` is `0`, `1`, or `2`.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
export ggplot_axes

split_ggargs(args...) = (),args
split_ggargs(ax::AxisId,args...) = (ax,),args
split_ggargs(ax1::AxisId,ax2::AxisId,args...) = (ax1,ax2),args
split_ggargs(ax1::AxisId,ax2::AxisId,ax3::AxisId,args...) = (ax1,ax2,ax3),args
split_ggargs(ax1::AxisId,ax2::AxisId,ax3::AxisId,ax4::AxisId,args...) =
  (ax1,ax2,ax3,ax4),args
split_ggargs(ax1::AxisId,ax2::AxisId,ax3::AxisId,ax4::AxisId,ax5::AxisId,args...) =
  (ax1,ax2,ax3,ax4,ax5),args
split_ggargs(ax1::AxisId,ax2::AxisId,ax3::AxisId,ax4::AxisId,ax5::AxisId,
  ax6::AxisId,args...) = (ax1,ax2,ax3,ax4,ax5,ax6),args
function split_ggargs(ax1::AxisId,ax2::AxisId,ax3::AxisId,ax4::AxisId,
                      ax5::AxisId,ax6::AxisId,ax7::AxisId,args...)
  error("Plotting data using 7 or more axes is not supported by ggplot backend.")
end

function ggplot_axes(data,args...;kwds...)
  ax, args = split_ggargs(args...)
  df, axes = asplotable(data,ax...;kwds...)
  if eltype(df.value) <: Complex
    df.value = abs.(df.value)
    @warn("Ignoring phase of complex value")
  end

  ggplot_axes_(df,axes,names(df)[2:end]...;args=args)
end

function ggplot_axes_(df,axes,x;args)
  R"""
  library(ggplot2)
  ggplot($df,aes_string(x=$(string(x)),y="value")) + geom_line()
  """
end

function ggplot_axes_(df,axes,x,y;args)
  R"""
  library(ggplot2)
  ggplot($df,aes_string(x=$(string(x)),y=$(string(y)))) +
      geom_raster(aes(fill=value))
  """
end

function ggplot_axes_(df,axes,x,y,z;args)
  R"""
  library(ggplot2)
  ggplot($df,aes_string(x=$(string(x)),y=$(string(y)))) +
      geom_raster(aes(fill=value)) +
      facet_wrap(as.formula(paste("~", $(string(z)))),
        labeller="label_both")
  """
end

function ggplot_axes_(df,axes,x,y,z,w;args)
  R"""
  library(ggplot2)
  ggplot($df,aes_string(x=$(string(x)),y=$(string(y)))) +
      geom_raster(aes(fill=value)) +
      facet_grid(as.formula(paste($(string(w)), "~",
        $(string(z)))),labeller="label_both")
  """
end

function ggplot_axes_(df,axes,x,y,z,w,v;args)
  R"""
  library(ggplot2)
  ggplot($df,aes_string(x=$(string(x)),y=$(string(y)))) +
      geom_raster(aes(fill=value)) +
      facet_grid(as.formula(paste($(string(w)), "~",
        $(string(z)), "+", $(string(v)))),labeller="label_both")
  """
end

function ggplot_axes_(df,axes,x,y,z,w,v,u;args)
  R"""
  library(ggplot2)
  ggplot($df,aes_string(x=$(string(x)),y=$(string(y)))) +
      geom_raster(aes(fill=value)) +
      facet_grid(as.formula(paste($(string(w)), "+", $(string(u)), "~",
        $(string(z)), "+", $(string(v)))),labeller="label_both")
  """
end

function ggplot_axes_(df,axes,args...;kwds...)
  error("Plotting data with $(length(axes)) dims along
        $(length(args)) axes is not supported.")
end


text:
"""
    check_MD5( file_path, checksum )
returns a MD5 hash from a file location.
Note: this converts Int8 representations to comma delimitted strings.
"""
get_MD5( file_path ) =  join( string.( open(md5, file_path) ), "," )

"""
    check_MD5( file_path, checksum )
Checks the result of an MD5 hash vs a stored checksum.
Note: this converts Int8 representations to comma delimitted strings.
"""
check_MD5( file_path, check_sum ) = get_MD5( file_path ) == check_sum


text:
using Documenter, Feynman

makedocs(
    modules = [Feynman],
    sitename = "Feynman.jl",
    pages = Any[
        "Home" => "index.md",
        "API" => "api.md",
    ],
)

deploydocs(
    repo = "github.com/bhgomes/Feynman.jl.git",
    target = "build",
)


text:
using Klimakoffer

world = Klimakoffer.read_geography(joinpath(@__DIR__,"..","input","world","The_World128x65.dat"),128,65)

cleared_map = Klimakoffer.clear_map(world)

clear = Klimakoffer.read_geography(joinpath(@__DIR__,"..","examples","reference_data","clear_map128x65.dat"),128,65)

diff = clear - cleared_map

res = sum(abs.(diff))


text:
mutable struct Friction{T} <: Contact{T}
    Nx::Adjoint{T,SVector{3,T}}
    D::SMatrix{2,6,T,12}
    cf::T
    offset::SVector{6,T}
    b::SVector{2,T}

    function Friction(body::Body{T}, normal::AbstractVector, cf::Real; offset::AbstractVector = zeros(3)) where T
        @assert cf>0

        # Derived from plane equation a*v1 + b*v2 + distance*v3 = p - offset
        V1, V2, V3 = orthogonalcols(normal) # gives two plane vectors and the original normal axis
        A = [V1 V2 V3]
        Ainv = inv(A)
        ainv3 = Ainv[3,SA[1; 2; 3]]
        Nx = ainv3'
        D = [[V1 V2];szeros(3,2)]'
        offset = [offset;0.0;0.0;0.0]

        new{T}(Nx, D, cf, offset, zeros(2)), body.id
    end
end

function Base.show(io::IO, mime::MIME{Symbol("text/plain")}, constraint::Friction{T}) where {T}
    summary(io, constraint)
    println(io,"")
    println(io, " Nx:     "*string(constraint.Nx))
    println(io, " D:      "*string(constraint.D))
    println(io, " cf:     "*string(constraint.cf))
    println(io, " offset: "*string(constraint.offset))
end

@inline additionalforce(friction::Friction) = friction.D'*friction.b
@inline function calcFrictionForce!(mechanism, ineqc, friction::Friction, i, body::Body)
    calcFrictionForce!(mechanism, friction, body, ineqc.γsol[2][i])
    return
end


text:
#=
P14 (*) Duplicate the elements of a list.
Example:
* (dupli '(a b c c d))
(A A B B C C C C D D)
=#

function dupli(coll)
    vcat(coll, coll)
end

@assert dupli([:a, :b, :c, :c, :d]) == [:a, :b, :c, :c, :d, :a, :b, :c, :c, :d]
@assert dupli([:a]) == [:a, :a]
@assert dupli([]) == []

println("Tests passed: JL-14.jl")


text:
using Test
using CSV
using DataFrames
using StatsBase
using BenchmarkTools
using PrettyTables

using ScottishTaxBenefitModel
using ScottishTaxBenefitModel.RunSettings: Settings, MT_Routing
using .STBOutput: make_poverty_line, summarise_inc_frame,
    dump_frames, summarise_frames, make_gain_lose

@testset "basic gain lose" begin
    pre = DataFrame( income=[1,2,3], weighted_people=[1,1,1]);
    post = DataFrame( income=[3,2,1], weighted_people=[1,1,1]);
    gl = make_gain_lose( pre, post, :income )
    @test gl.gainers == 1
    @test gl.losers == 1
    @test gl.nc == 1
    @test gl.popn == 3
end

@testset "povery line" begin

end

max_stars_repo_name: grahamstark/ScottishTaxBenefitModel

max_issues_repo_name: grahamstark/ScottishTaxBenefitModel

max_forks_repo_name: grahamstark/ScottishTaxBenefitModel


text:
# Chapter 18. Subtyping
# https://benlauwens.github.io/ThinkJulia.jl/latest/book.html#chap18

include("AllnixThinkJulia.jl")
import .AllnixThinkJulia
ns = AllnixThinkJulia

# Cards
# https://benlauwens.github.io/ThinkJulia.jl/latest/book.html#_cards
queen_of_diamonds = ns.Card(2, 12)

# Global Variables
ns.Card(3,11)

# Comparing Cards
# https://benlauwens.github.io/ThinkJulia.jl/latest/book.html#_comparing_cards

# Exercise 18-1

# Unit Testing
# https://benlauwens.github.io/ThinkJulia.jl/latest/book.html#_unit_testing
using Test

@test ns.Card(1,4) < ns.Card(2,4)
@test ns.Card(1,3) < ns.Card(1,4)

# Decks
# https://benlauwens.github.io/ThinkJulia.jl/latest/book.html#_decks
# See ch18module.jl

# Add, Remove, Shuffle and Sort
import Random.shuffle!
deck = ns.Deck()
pop!(deck)
deck
shuffle!(deck)

# Exercise 18-2
sort!(deck)

# Abstract Types and Subtyping
# https://benlauwens.github.io/ThinkJulia.jl/latest/book.html#_abstract_types_and_subtyping

# Subtyping: Group related concrete types <- Abstract Type
deck = ns.Deck()
deck isa ns.CardSet

hand = ns.Hand("new hand")

# Abstract Types and Functions
# https://benlauwens.github.io/ThinkJulia.jl/latest/book.html#_abstract_types_and_functions
deck = ns.Deck()
shuffle!(deck)
card = pop!(deck)
push!(hand, card)

ns.move!(deck, hand, 4)
@which ns.move!(deck, hand, 4)
supertype(ns.Deck)


text:
export mass
export Engine, thrust, Isp, thrust
export Rocket, update!

"A rocket produce thrust by consuming propellant"
abstract type Rocket end

"Provides propulsion for a rocket"
abstract type Engine end

"""
    mass(r::Rocket)
Mass of given object including its subcomponents.
"""
function mass end

"""
    thrust(r::Rocket)
    thrust(e::Engine)
An engine produces thrust. A rocket will have the same thrust as the sum of the engines
attached to the rocket.
"""
function thrust end

"""
    update!(r::Rocket)
Consume propellant and update mass accordingly
"""
function update! end

"""
    Isp(e::Engine)
Specific Impulse for a rocket engine. Says something about the equivalent of the gas milage
of a rocket. An engine with higher Isp will be able to change the velocity of a rocket more
per unit of fuel, than an engine with lower Isp
"""
function Isp end


text:
# test bearing

using RoME
using Test
using TensorCast
using DistributedFactorGraphs
using Manifolds: hat

@testset "Testing Bearing2D factor" begin
##
M = SpecialEuclidean(2)
ϵ = identity_element(M)
ps = [exp(M, ϵ,  hat(M, ϵ, [0.,0,0]))]
push!(ps, exp(M, ϵ,  hat(M, ϵ, [5.,0,0])))
push!(ps, exp(M, ϵ,  hat(M, ϵ, [10.,0,0])))
push!(ps, exp(M, ϵ,  hat(M, ϵ, [10.,5,0])))
push!(ps, exp(M, ϵ,  hat(M, ϵ, [10.,10,0])))
push!(ps, exp(M, ϵ,  hat(M, ϵ, [5.,10,0])))
push!(ps, exp(M, ϵ,  hat(M, ϵ, [0.,10,0])))
push!(ps, exp(M, ϵ,  hat(M, ϵ, [0.,5,0])))

push!(ps, exp(M, ϵ,  hat(M, ϵ, [0.,0,pi/4])))
push!(ps, exp(M, ϵ,  hat(M, ϵ, [0.,0,-pi/4])))

push!(ps, exp(M, ϵ,  hat(M, ϵ, [1.,2,0]))) # [4,3,0]
# push!(ps, exp(M, ϵ,  hat(M, ϵ, [1.,2,pi/2]))) # [3,-4,-pi/2]
# push!(ps, exp(M, ϵ,  hat(M, ϵ, [1.,2,pi]))) # [-4,-3,-pi]
# push!(ps, exp(M, ϵ,  hat(M, ϵ, [1.,2,-pi]))) # [-3,4,pi]

rs = [0, -pi/4, -pi/2, -3pi/4, -pi, 3pi/4, pi/2, pi/4, pi/4, -pi/4]

push!(rs, pi/4 - atan(3,4))

q = [5., 5]
m = [pi/4]

f = Pose2Point2Bearing(Normal(pi/4,0.05))

for (i,p) in enumerate(ps)
    # @show p
    r = calcFactorResidualTemporary(f, (Pose2, Point2), m, (p, q))
    # @warn "comp: $(isapprox(r, rs[i], atol=0.001))" a=r b=rs[i]
    @test isapprox(r, rs[i], atol=0.001)
end

@warn "Bearing2D, must still test factor gradients, which will also verify the sign on residual function calculations"
@test_broken false

# test sign
f = Pose2Point2Bearing(Normal(pi/2,0.001))
xi = ProductRepr([0.,0], [1. 0; 0 1])
xj = [1.,1]
res = calcFactorResidualTemporary(f, (Pose2, Point2), [], (xi, xj))
@test isapprox(res, pi/4, atol=0.1)

# test -pi +pi case
f = Pose2Point2Bearing(Normal(pi,0.001))
xi = ProductRepr([0.,0], [1. 0; 0 1])
xj = [-1, -0.001]
res = calcFactorResidualTemporary(f, (Pose2, Point2), [pi], (xi, xj))
@test isapprox(res, -0.001, atol=1e-3)

##
end

@testset "Simple Bearing2D test to give narrow vs broad belief posteriors" begin

# new factor graph
fg = initfg()

addVariable!(fg, :x1, Pose2)
addVariable!(fg, :x2, Pose2)
addFactor!(fg, [:x1], PriorPose2(MvNormal([10.0; 0; 0], diagm([0.01;0.01;0.001].^2))))
addFactor!(fg, [:x2], PriorPose2(MvNormal([0.0;10.0; 0], diagm([0.01;0.01;0.001].^2))))

# limit the scope nearby origin, but test was written because it was found with Bearing
# that many points would suddenly land very far away from the origin, even with the prior.
addVariable!(fg, :l1, Point2)
addFactor!(fg, [:l1], PriorPoint2(MvNormal([0.,0], diagm([10.,10].^2))))

initAll!(fg)

addFactor!(fg, [:x1;:l1], Pose2Point2Bearing(Normal(pi, 0.001)))
addFactor!(fg, [:x2;:l1], Pose2Point2Bearing(Normal(-pi/2,0.001)))

points = approxConv(fg, :x1l1f1, :l1)
@cast pts[j,i] := points[i][j]
#x1l1 "partial" constrian on x-axis
#FIXME very pessimistic tests, should be way better
@test sum(abs.(pts[2,:]) .< 100) > 30

L1 = approxConvBelief(fg, :x2l1f1, :l1)
points = getPoints(L1, false)
@cast pts[j,i] := points[i][j]
@test sum(abs.(pts[1,:]) .< 100) > 30

solveTree!(fg)

points = getPoints(getBelief(fg, :l1))
@cast pts[j,i] := points[i][j]
pts = collect(pts)
#FIXME check test after Bearing2D is fixed
@test_broken all(sum(abs.(pts) .< [10,10],dims=2) .> [80,80])

end

@testset "Triangulation test in 2D, 3 beacons" begin

# noise models
lmp_noise = Matrix(Diagonal([0.01;0.01].^2))

# new factor graph
fg = initfg()

# landmarks
addVariable!(fg, :l1, Point2)
addVariable!(fg, :l2, Point2)
addVariable!(fg, :l3, Point2)

addFactor!(fg, [:l1], PriorPoint2(MvNormal([10.0;1.0],lmp_noise)), graphinit=false)
addFactor!(fg, [:l2], PriorPoint2(MvNormal([10.0+sqrt(3)/2;-0.5],lmp_noise)), graphinit=false)
addFactor!(fg, [:l3], PriorPoint2(MvNormal([10.0-sqrt(3)/2;-0.5],lmp_noise)), graphinit=false)

# pose
addVariable!(fg, :x1, Pose2)
addFactor!(fg, [:x1;:l1], Pose2Point2Bearing(Normal(pi/2,0.05)), graphinit=false)
addFactor!(fg, [:x1;:l2], Pose2Point2Bearing(Normal(-pi/6,0.05)), graphinit=false)
addFactor!(fg, [:x1;:l3], Pose2Point2Bearing(Normal(-pi+pi/6,0.05)), graphinit=false)

# #TODO remove easier just for tests
# # landmarks
# addFactor!(fg, [:l1], PriorPoint2(MvNormal([10.0;0.0],lmp_noise)))
# addFactor!(fg, [:l2], PriorPoint2(MvNormal([0.0;10.0],lmp_noise)))
# addFactor!(fg, [:l3], PriorPoint2(MvNormal([10.0,10],lmp_noise)))
# # pose
# addVariable!(fg, :x1, Pose2)
# addFactor!(fg, [:x1;:l1], Pose2Point2Bearing(Normal(0.0,0.05)), graphinit=false)
# addFactor!(fg, [:x1;:l2], Pose2Point2Bearing(Normal(pi/2,0.05)), graphinit=false)
# addFactor!(fg, [:x1;:l3], Pose2Point2Bearing(Normal(pi/4,0.05)), graphinit=false)

# initManual!(fg, :x1, [30.0*randn(2,100);randn(1,100)])

## Look at graph

# drawGraph(fg)

## solve

# getSolverParams(fg).drawtree = true
# getSolverParams(fg).showtree = true

tree = solveTree!(fg)

## Look at results
#
# using RoMEPlotting
# Gadfly.set_default_plot_size(35cm,25cm)
# #
# drawPosesLandms(fg, spscale=0.2) # |> PNG("/home/dehann/Downloads/triangulation.png",20cm,15cm)
# #
# plotPose(fg, :x1, scale=0.1) # |> PNG("/home/dehann/Downloads/triangulationPose.png",20cm,15cm)
# #
# plotFactor(fg,  lsf(fg, Pose2Point2Bearing)[1])
# #
# plotLocalProduct(fg, :x1)

## complete the unit test

points = getPoints(getBelief(fg, :x1))
mean(SpecialEuclidean(2), points)

@cast pts[j,i] := getCoordinates.(Pose2, points)[i][j]
pts = collect(pts)
pts[1,:] .-= 10.0

N = size(pts,2)

@test_broken 0.8*N < sum(sqrt.(sum(pts[1:2,:].^2,dims=1)) .< 0.3)

@test_broken 0.8*N < sum(abs.(pts[3,:]) .< 0.1)

#

end

@testset "Triangulation test in 2D (opposite), 3 beacons" begin

# noise models
lmp_noise = Matrix(Diagonal([0.01;0.01].^2))

# new factor graph
fg = initfg()

# landmarks
addVariable!(fg, :l1, Point2)
addVariable!(fg, :l2, Point2)
addVariable!(fg, :l3, Point2)

addFactor!(fg, [:l1], PriorPoint2(MvNormal([-10.0;1.0-10.0],lmp_noise)), graphinit=false)
addFactor!(fg, [:l2], PriorPoint2(MvNormal([-10.0+sqrt(3)/2;-0.5-10.0],lmp_noise)), graphinit=false)
addFactor!(fg, [:l3], PriorPoint2(MvNormal([-10.0-sqrt(3)/2;-0.5-10.0],lmp_noise)), graphinit=false)

# pose
addVariable!(fg, :x1, Pose2)
addFactor!(fg, [:x1;:l1], Pose2Point2Bearing(Normal(pi/2,0.05)), graphinit=false)
addFactor!(fg, [:x1;:l2], Pose2Point2Bearing(Normal(-pi/6,0.05)), graphinit=false)
addFactor!(fg, [:x1;:l3], Pose2Point2Bearing(Normal(-pi+pi/6,0.05)), graphinit=false)

# initManual!(fg, :x1, [0.01*randn(2,100);-randn(1,100)])
# initManual!(fg, :x1, [30.0*randn(2,100);randn(1,100)])

# drawGraph(fg)

tree = solveTree!(fg)

points = getPoints(getBelief(fg, :x1))
mean(SpecialEuclidean(2), points)

@cast pts[j,i] := getCoordinates.(Pose2, points)[i][j]
pts = collect(pts)

pts[1,:] .+= 10.0
pts[2,:] .+= 10.0

N = size(pts,2)

@test_broken 0.8*N < sum(sqrt.(sum(pts[1:2,:].^2,dims=1)) .< 0.3)

@test_broken 0.8*N < sum(abs.(pts[3,:]) .< 0.1)

#
# using RoMEPlotting
# Gadfly.set_default_plot_size(35cm,25cm)
# drawPosesLandms(fg, spscale=0.2) # |> PNG("/home/dehann/Downloads/triangulation.png",20cm,15cm)
# plotPose(fg, :x1, scale=0.1) # |> PNG("/home/dehann/Downloads/triangulationPose.png",20cm,15cm)
#
# plotLocalProduct(fg, :x1)
#
# plotFactor(fg,  lsf(fg, Pose2Point2Bearing)[1])

end

#


text:
using LanguageServer
using CSTParser
using Base.Test
Range = LanguageServer.Range

@testset "LanguageServer" begin

include("test_document.jl")
include("test_communication.jl")
include("test_hover.jl")
include("test_lint.jl")
include("test_scope.jl")

end


text:
# ---
# title: 788. Rotated Digits
# id: problem788
# author: Indigo
# date: 2021-01-29
# difficulty: Easy
# categories: String
# link: <https://leetcode.com/problems/rotated-digits/description/>
# hidden: true
# ---
#
# X is a good number if after rotating each digit individually by 180 degrees,
# we get a valid number that is different from X.  Each digit must be rotated -
# we cannot choose to leave it alone.
#
# A number is valid if each digit remains a digit after rotation. 0, 1, and 8
# rotate to themselves; 2 and 5 rotate to each other (on this case they are
# rotated in a different direction, in other words 2 or 5 gets mirrored); 6 and
# 9 rotate to each other, and the rest of the numbers do not rotate to any other
# number and become invalid.
#
# Now given a positive number `N`, how many numbers X from `1` to `N` are good?
#
#
#
#     **Example:**
#     Input: 10
#     Output: 4
#     Explanation:
#     There are four good numbers in the range [1, 10] : 2, 5, 6, 9.
#     Note that 1 and 10 are not good numbers, since they remain unchanged after rotating.
#
#
# **Note:**
#
#   * N  will be in range `[1, 10000]`.
#
#
## @lc code=start
using LeetCode

function rotated_digits(N::Int)
    if N <= 9
        return count(==(1), [0, 1, -1, -1, 1, 1, -1, 0, 1][1:N])
    end
    ans, dp = 0, vcat([0, 0, 1, -1, -1, 1, 1, -1, 0, 1], fill(0, N - 9))
    for i in 1:N-1
        dp[i + 1] = (-1 in [dp[i ÷ 10 + 1], dp[i % 10 + 1]]) ? -1 : (dp[i ÷ 10 + 1] | dp[i % 10 + 1])
        ans += (dp[i + 1] == 1)
    end
    ans
end
## @lc code=end


text:
import FFAST # for mass absorption coefficienct
import CSV

"""
FFAST represents an implementation of mass-absorption coefficients and associated
edge and other elemental data.  It is a thin wrapper around the FFAST library.
"""
struct FFASTDB <: NeXLAlgorithm end

mac(elm::Element, energy::Float64, ::Type{FFASTDB})::Float64 =
    FFAST.mac(FFAST.PhotoElectricMAC, z(elm), energy)

function macU(elm::Element, energy::Float64, ::Type{FFASTDB})::UncertainValue
    macv = FFAST.mac(FFAST.PhotoElectricMAC, z(elm), energy)
    return uv(
        macv,
        min(FFAST.fractionaluncertainty(FFAST.SolidLiquid, z(elm), energy)[1], 0.9) * macv,
    )
end

let superset_edge_energies_data, subshellsindexes_data
    function load_supersetedgeenergies()
        @info "Loading William's edge energies."
        res = Dict{Tuple{Int,Int}, Float64}()
        for (z, row) in enumerate(CSV.File(joinpath(dirname(pathof(@__MODULE__)), "..", "data", "WilliamsBinding.csv"), header=true))
            for i in 1:29
                if(!ismissing(row[i]))
                    res[(z, i)]=row[i]
                end
            end
        end
        @info "Replacing with FFAST where available..."
        for (idx, _) in res
            if (idx[1] in FFAST.eachelement()) && FFAST.hasedge(idx...)
                res[idx] = FFAST.edgeenergy(idx...)
            end
        end
        return res
    end
    superset_edge_energies_data = load_supersetedgeenergies()
    subshellsindexes_data = Dict(z=>filter(ss->haskey(superset_edge_energies_data, (z,ss)), 1:29) for z in 1:99)

    global superset_edge_energies() = superset_edge_energies_data
    global subshellindices(z::Int, ::Type{FFASTDB}) = subshellsindexes_data[z]
end

edgeenergy(z::Int, ss::Int, ::Type{FFASTDB})::Float64 = superset_edge_energies()[(z,ss)]

hasedge(z::Int, ss::Int, ::Type{FFASTDB})::Bool = haskey(superset_edge_energies(), (z,ss))

# subshellindices(z::Int, ::Type{FFASTDB}) = filter(ss->haskey(superset_edge_energies(), (z,ss)), 1:29)

"""
    energy(z::Int, inner::Int, outer::Int)::Float64

Return energy (in eV) of the transition by specified inner and outer sub-shell index.
"""
energy(z::Int, inner::Int, outer::Int, ::Type{FFASTDB})::Float64 =
    edgeenergy(z, inner, FFASTDB) - edgeenergy(z, outer, FFASTDB)

"""
    jumpratio(z::Int, ss::Int, ::Type{FFASTDB}) =

Returns the jump ratio for the specified shell by index K=>1, L1=>2, L2=>3, ..., M5=>9.
"""
jumpratio(z::Int, ss::Int, ::Type{FFASTDB}) = FFAST.jumpratio(z, ss)


text:
"""
A Tiler is an iterator that, for each iteration, returns a tuple of:

- the `x`/`y` point of the center of each tile in a set of tiles that divide up a rectangular space such as a page into rows and columns (relative to current 0/0)

- the number of the tile

    tiles = Tiler(areawidth, areaheight, nrows, ncols, margin=20)

where `width`, `height` is the dimensions of the area to be tiled, `nrows`/`ncols`
is the number of rows and columns required, and `margin` is applied to all four
edges of the area before the function calculates the tile sizes required.

    tiles = Tiler(1000, 800, 4, 5, margin=20)
    for (pos, n) in tiles
    # the point pos is the center of the tile
    end

You can access the calculated tile width and height like this:

    tiles = Tiler(1000, 800, 4, 5, margin=20)
    for (pos, n) in tiles
      ellipse(pos.x, pos.y, tiles.tilewidth, tiles.tileheight, :fill)
    end
"""
type Tiler
  width::Real
  height::Real
  tilewidth::Real
  tileheight::Real
  nrows::Int
  ncols::Int
  margin::Real
  function Tiler(areawidth, areaheight, nrows::Int, ncols::Int; margin=10)
      tilewidth  = (areawidth - 2margin)/ncols
      tileheight = (areaheight - 2margin)/nrows
      new(areawidth, areaheight, tilewidth, tileheight, nrows, ncols, margin)
  end
end

function Base.start(pt::Tiler)
# return the initial state
  x = -(pt.width/2)  + pt.margin + (pt.tilewidth/2)
  y = -(pt.height/2) + pt.margin + (pt.tileheight/2)
  return (Point(x, y), 1)
end

function Base.next(pt::Tiler, state)
  # Returns the item and the next state
  # state[1] is the Point
  x = state[1].x
  y = state[1].y
  # state[2] is the tilenumber
  tilenumber = state[2]
  x1 = x + pt.tilewidth
  y1 = y
  if x1 > (pt.width/2) - pt.margin
    y1 += pt.tileheight
    x1 = -(pt.width/2) + pt.margin + (pt.tilewidth/2)
  end
  return ((Point(x, y), tilenumber), (Point(x1, y1), tilenumber + 1))
end

function Base.done(pt::Tiler, state)
  # Tests if there are any items remaining
  state[2] > (pt.nrows * pt.ncols)
end


text:
# residual calculations

function Rex(mp::modelpars_2state, Δω)
    kex = mp.kex
    pB = mp.pB
    @. real(0.5*(kex+1im*Δω-sqrt((kex+1im*Δω)^2-4im*kex*pB*Δω)))
end

"""
simulate R2 hahn-echo relaxation rates, for 1H coherence order pH, and list of field strengths B0
"""
function sim_R2(mp::modelpars_2state, ΔδH, ΔδC, R20, pH, B0)
    Δω = 1e-6 * (ΔδC*γC + pH*ΔδH*γH) .* B0
    return R20 .+ Rex(mp, Δω)
end

function sim_R2s(mp, ΔδH, ΔδC, R2_ZQ, R2_DQ, R2_QQ⁻, R2_QQ⁺, B0s)
    dq = sim_R2(mp, ΔδH, ΔδC, R2_DQ, +1, B0s)
    zq = sim_R2(mp, ΔδH, ΔδC, R2_ZQ, -1, B0s)
    dqq = sim_R2(mp, ΔδH, ΔδC, R2_QQ⁻, -3, B0s)
    qq = sim_R2(mp, ΔδH, ΔδC, R2_QQ⁺, +3, B0s)
    return R2Data(zq, dq, dqq, qq, B0s)
end

"""
simulate CPMG relaxation rates
"""
function sim_CPMG(expt::CPMGData_MQ, mp::modelpars_2state, ΔδH, ΔδC, R20)
    # basis for MQ CPMG experiment is (ZQ⁻, ZQ⁺, DQ⁻, DQ⁺) where +/- indicates 1H coherence order
    # initial magnetisation is (1/2) (ZQ⁺ + DQ⁺)
    # final magnetisation read-out is (1/2) (ZQ⁻ + DQ⁻)

    # pulse program parameters
    T = expt.Trelax
    B0 = expt.B0
    νCPMG = expt.νCPMG
    kex = mp.kex
    pB = mp.pB

    ΔωC = γC * B0 * ΔδC * 1e-6
    ΔωH = γH * B0 * ΔδH * 1e-6
    ΔωZQ = ΔωH - ΔωC
    ΔωDQ = ΔωH + ΔωC

    # initial magnetisation
    M0 = 0.5 * [0, (1-pB), 0, (1-pB), 0, pB, 0, pB]
    Mobs = [1 0 1 0 1 0 1 0]

    # set up propagators
    UπH = [
        0 0 0 1
        0 0 1 0
        0 1 0 0
        1 0 0 0
        ] |> float
    UπC = [
        0 0 1 0
        0 0 0 1
        1 0 0 0
        0 1 0 0
        ] |> float
    UπH = kron(I(2), UπH)
    UπC = kron(I(2), UπC)
    kAB = pB*kex
    kBA = (1-pB)*kex
    K = kron([-kAB kBA
               kAB -kBA
               ], I(4))
    Ω = Diagonal([0, 0, 0, 0, -ΔωZQ, ΔωZQ, -ΔωDQ, ΔωDQ])
    L = -1im*Ω + K

    R2eff = similar(νCPMG)
    for j = 1:length(νCPMG)
        ncyc = νCPMG[j] * T
        τ = T / 4ncyc

        Uτ = exp(L*τ)
        Uecho = Uτ * UπC * Uτ
        UechoN = Uecho^ncyc
        U = UechoN * UπH * UechoN

        relintensity = real(Mobs*U*M0)[1]
        R2eff[j] = R20 - 1/T*log(relintensity < 1e-6 ? 1e-6 : relintensity)
    end
    return CPMGData_MQ(νCPMG, R2eff, B0, T)
end

function sim_CPMG(expt::CPMGData_SQ1H, mp::modelpars_2state, ΔδH, ΔδC, R20)
    # basis for 1H SQ CPMG experiment is (Hx, Hy)
    # initial magnetisation is Hx
    # final magnetisation read-out is Hx

    # pulse program parameters
    T = expt.Trelax
    B0 = expt.B0
    νCPMG = expt.νCPMG
    kex = mp.kex
    pB = mp.pB

    ΔωH = γH * B0 * ΔδH * 1e-6

    # initial magnetisation
    M0 = [(1-pB), 0, pB, 0]
    Mobs = -[1 0 1 0]

    # set up propagators
    UπHx = [
        1 0
        0 -1
        ] |> float
    UπHy = [
        -1 0
        0 1
        ] |> float
    UπHx = kron(I(2), UπHx)
    UπHy = kron(I(2), UπHy)
    kAB = pB*kex
    kBA = (1-pB)*kex
    K = kron([
        -kAB kBA
        kAB -kBA
        ], I(2))
    Ω = [
        0 0 0 0
        0 0 0 0
        0 0 0 -ΔωH
        0 0 ΔωH 0
        ]
    L = Ω + K

    R2eff = similar(νCPMG)
    for j = 1:length(νCPMG)
        ncyc = νCPMG[j] * T
        τ = T / 4ncyc

        Uτ = exp(L*τ)
        Uecho = Uτ * UπHx * Uτ
        UechoN = Uecho^ncyc
        U = UechoN * UπHy * UechoN
        relintensity = real(Mobs*U*M0)[1]
        R2eff[j] = R20 - 1/T*log(relintensity < 1e-6 ? 1e-6 : relintensity)
    end
    return CPMGData_SQ1H(νCPMG, R2eff, B0, T)
end

function data(R2data::R2Data)
    return vec([R2data.ZQ; R2data.DQ; R2data.QQ⁻; R2data.QQ⁺])
end
function data(CPMGdata::T) where T<:CPMGData
    return CPMGdata.R2eff
end

function residuals_R2(expt::R2Data, mp, ΔδH, ΔδC, R2_ZQ, R2_DQ, R2_QQ⁻, R2_QQ⁺, B0s)
    sim = sim_R2s(mp, ΔδH, ΔδC, R2_ZQ, R2_DQ, R2_QQ⁻, R2_QQ⁺, B0s)
    return zscore.(data(expt) .- data(sim_R2s(mp, ΔδH, ΔδC, R2_ZQ, R2_DQ, R2_QQ⁻, R2_QQ⁺, B0s)))
end

function residuals_CPMG(expt::CPMGData, mp, ΔδH, ΔδC, R20)
    return zscore.(data(expt) .- data(sim_CPMG(expt, mp, ΔδH, ΔδC, R20)))
end

# """
# return residuals for particular coherence order pH
# """
# function residuals_R2(R2obs, mp, ΔδH, ΔδC, R20, pH, B0)
#     return zscore.(sim_R2(mp, ΔδH, ΔδC, R20, pH, B0) .- R2obs)
# end

# """
# return list of residuals for R2 measurements (all B0, all coherence orders)
# """
# function residuals_R2(R2obs::R2set, mp, ΔδH, ΔδC, R20s::R2set)
#     sim = sim_R2s(mp, ΔδH, ΔδC, R20s)
#     return vec([
#         zscore.(R2obs.ZQ - sim.ZQ)
#         zscore.(R2obs.DQ - sim.DQ)
#         zscore.(R2obs.QQ⁻ - sim.QQ⁻)
#         zscore.(R2obs.QQ⁺ - sim.QQ⁺)
#     ])
# end

max_stars_repo_name: chriswaudby/qq-hahnecho

max_issues_repo_name: chriswaudby/qq-hahnecho

max_forks_repo_name: chriswaudby/qq-hahnecho