text:
"""
Much of the code in this file is adapted from LIONSIMBA, a li-ion simulation toolbox in MATLAB.
See https://github.com/lionsimbatoolbox/LIONSIMBA for more information.
"""

function residuals_PET!(residuals, t, x, ẋ, p::AbstractModel)
    """
    First put the vector of x, ẋ, and residuals into dictionaries
    """
    states  = retrieve_states(x, p)
    ∂states = retrieve_states(ẋ, p)
    res     = retrieve_states(residuals, p)

    states[:t] = t
    states[:x] = x
    states[:ẋ] = ẋ

    """
    Calculate a few necessary auxiliary variables
    """
    build_auxiliary_states!(states, p)

    """
    Differential residuals
    """
    # Residual for the electrolyte concentration, c_e
    residuals_c_e!(res, states, ∂states, p)

    # Residual for the average surface concentration, c_s_avg
    residuals_c_s_avg!(res, states, ∂states, p)

    # If the polynomial fit for the surface concentration is used, Q
    if p.numerics.solid_diffusion == :polynomial
        residuals_Q!(res, states, ∂states, p)
    end

    # Lithium plating film thickness, film
    if p.numerics.aging == :SEI
        residuals_film!(res, states, ∂states, p)
        residuals_SOH!(res, states, ∂states, p)
    end

    # Check if the thermal dynamics are enabled.
    if p.numerics.temperature == true
        build_heat_generation_rates!(states, p)
        residuals_T!(res, states, ∂states, p)
    end

    """
    Algebraic residuals
    """
    # Residuals for ionic flux, j
    residuals_j!(res, states, p)

    # Residuals for side reaction ionic flux, j_s
    if p.numerics.aging ∈ (:SEI, :R_aging)
        residuals_j_s!(res, states, p)
    end

    # Residuals for the electrolyte potential, Φ_e
    residuals_Φ_e!(res, states, p)

    # Residuals for the solid potential, Φ_s
    residuals_Φ_s!(res, states, p)

    # Residuals for applied current density, I [not currently used]
    residuals_scalar!(res, states, p)

    """
    Compile all residuals together
    """
    build_residuals!(residuals, res, p)

    return nothing
end
function residuals_PET!(p::AbstractModel)
    θ_sym, Y_sym, YP_sym, t_sym, SOC_sym, X_applied, γ_sym, p_sym, θ_keys = get_symbolic_vars(p)
    res = similar(Y_sym)
    residuals_PET!(res, t_sym, Y_sym, YP_sym, p_sym)
end

function residuals_c_e!(res, states, ∂states, p::AbstractModel)
    """
    Calculate the electrolyte concentration residuals [mol/m³]
    """

    c_e = states[:c_e]
    T   = states[:T]
    j   = states[:j_total]

    ∂c_e = ∂states[:c_e]

    res_c_e = res[:c_e]

    Δx = Δx_values(p.N)

    # Diffusion coefficients and surface area to volume ratio
    D_eff_p, D_eff_s, D_eff_n = coeff_electrolyte_diffusion_effective(states, p)
    a_p, a_n = surface_area_to_volume_ratio(p)

    # Interpolation of the diffusion coefficients, same for electrolyte conductivities
    D_eff_p, D_eff_s, D_eff_n = interpolate_electrolyte_conductivities(D_eff_p, D_eff_s, D_eff_n, p)

    A_tot = -block_matrix_maker(p, D_eff_p, D_eff_s, D_eff_n)

    # dividing by the length and Δx
    ind = 1:p.N.p
    @views @inbounds A_tot[ind,ind] ./= (Δx.p*p.θ[:l_p])^2
    ind = (1:p.N.s) .+ ind[end]
    @views @inbounds A_tot[ind,ind] ./= (Δx.s*p.θ[:l_s])^2
    ind = (1:p.N.n) .+ ind[end]
    @views @inbounds A_tot[ind,ind] ./= (Δx.n*p.θ[:l_n])^2

    # Reset values on the lines for the interfaces conditions
    @views @inbounds A_tot[p.N.p,:]   .= 0.0
    @views @inbounds A_tot[p.N.p+1,:] .= 0.0

    # Reset values on the lines for the interfaces conditions
    @views @inbounds A_tot[p.N.p+p.N.s,:]   .= 0.0
    @views @inbounds A_tot[p.N.p+p.N.s+1,:] .= 0.0

    ## Interface between separator and cathode (last volume in the cathode)

    # Compute the common denominator at the interface
    @views @inbounds den_s = (Δx.p*p.θ[:l_p]/2 + Δx.s*p.θ[:l_s]/2)
    # Last diffusion coefficient of the cathode
    @views @inbounds last_p = D_eff_p[end-1]/(Δx.p*p.θ[:l_p])
    # Diffusion coefficient on the interface
    @views @inbounds first_s = D_eff_p[end]/den_s
    # Fix the values at the boundaries
    @views @inbounds A_tot[p.N.p,p.N.p-1:p.N.p+1] .= [last_p; -(last_p + first_s); first_s]/(Δx.p*p.θ[:l_p]*p.θ[:ϵ_p])

    ## Interface between separator and cathode (first volume in the separator)

    # First diffusion coefficient in the separator
    @views @inbounds second_s = D_eff_s[1]/(Δx.s*p.θ[:l_s])
    # Diffusion coefficient on the interface
    @views @inbounds first_s = D_eff_p[end]/den_s

    @views @inbounds A_tot[p.N.p+1,p.N.p:p.N.p+2] .= [first_s; -(first_s+second_s); second_s]/(Δx.s*p.θ[:l_s]*p.θ[:ϵ_s])

    ## Interface between separator and anode (last volume in the separator)

    # Compute the common denominator at the interface
    @views @inbounds den_s = (Δx.s*p.θ[:l_s]/2 + Δx.n*p.θ[:l_n]/2)
    # Last diffusion coefficient in the separator
    @views @inbounds last_s = D_eff_s[end-1]/(Δx.s*p.θ[:l_s])
    # Diffusion coefficient on the interface
    @views @inbounds first_n = D_eff_s[end]/den_s

    @views @inbounds A_tot[p.N.p+p.N.s,p.N.p+p.N.s-1:p.N.p+p.N.s+1] .= [last_s; -(last_s+first_n); first_n]/(Δx.s*p.θ[:l_s]*p.θ[:ϵ_s])

    ## Interface between separator and anode (first volume in the anode)

    # Compute the common denominator at the interface
    den_n = (Δx.s*p.θ[:l_s]/2 + Δx.n*p.θ[:l_n]/2)
    # First diffusion coefficient in the anode
    @views @inbounds second_n = D_eff_n[1]/(Δx.n*p.θ[:l_n])
    # Diffusion coefficient on the interface
    @views @inbounds first_n = D_eff_s[end]/den_n

    A_tot[p.N.p+p.N.s+1,p.N.p+p.N.s:p.N.p+p.N.s+2] .= [first_n; -(first_n+second_n); second_n]/(Δx.n*p.θ[:l_n]*p.θ[:ϵ_n])

    ϵ_tot = [
        ones(p.N.p).*p.θ[:ϵ_p]
        ones(p.N.s).*p.θ[:ϵ_s]
        ones(p.N.n).*p.θ[:ϵ_n]
        ]

    K = 1.0./ϵ_tot
    A_ϵ = zeros(eltype(K), (p.N.p+p.N.s+p.N.n), (p.N.p+p.N.s+p.N.n))
    ind_diagonal = diagind(A_ϵ)
    ind_neg1diag = diagind(A_ϵ, -1)
    ind_pos1diag = diagind(A_ϵ, 1)
    A_ϵ[ind_diagonal] .= K

    # Build porosities matrix
    @views @inbounds A_ϵ[ind_neg1diag] .= K[1:end-1]
    @views @inbounds A_ϵ[ind_pos1diag] .= K[1:end-1]

    @views @inbounds A_ϵ[p.N.p,p.N.p-1:p.N.p+1] .= 1.0
    @views @inbounds A_ϵ[p.N.p+1,p.N.p:p.N.p+2] .= 1.0

    @views @inbounds A_ϵ[p.N.p+p.N.s,p.N.p+p.N.s-1:p.N.p+p.N.s+1] .= 1.0
    @views @inbounds A_ϵ[p.N.p+p.N.s+1,p.N.p+p.N.s:p.N.p+p.N.s+2] .= 1.0

    A_tot = Matrix(A_tot).*Matrix(A_ϵ)

    # Write the RHS of the equation
    rhsCe = A_tot*c_e

    ind_p = (1:p.N.p)
    ind_n = (1:p.N.n) .+ (p.N.p+p.N.s)

    ν_p,ν_s,ν_n = p.numerics.thermodynamic_factor(c_e.p, c_e.s, c_e.n, T.p, T.s, T.n, p)

    rhsCe[ind_p] .+= K[ind_p].*(1-p.θ[:t₊]).*ν_p.*a_p.*j.p
    # nothing for the separator since a_s = 0
    rhsCe[ind_n] .+= K[ind_n].*(1-p.θ[:t₊]).*ν_n.*a_n.*j.n

    # Write the residual of the equation
    res_c_e .= rhsCe .- ∂c_e

    return nothing
end

function residuals_c_s_avg!(res, states, ∂states, p::T) where {T<:Union{AbstractModelSolidDiff{:polynomial},AbstractModelSolidDiff{:quadratic}}}
    """
    Calculate the solid particle concentrations residuals with quadratic or polynomial approximations [mol/m³]
    """
    j = states[:j]

    ∂c_s_avg = ∂states[:c_s_avg]

    res_c_s_avg = res[:c_s_avg]

    # Cathode
    rhsCs_p = -3j.p/p.θ[:Rp_p]

    # Anode
    rhsCs_n = -3j.n/p.θ[:Rp_n]

    res_c_s_avg .= [rhsCs_p; rhsCs_n] .- ∂c_s_avg

    return nothing
end
function residuals_c_s_avg!(res, states, ∂states, p::T) where {jac,temp,T<:AbstractModel{jac,temp,:Fickian,:finite_difference}}
    """
    Calculate the volume-averaged solid particle concentration residuals using a 9th order accurate finite difference method (FDM) [mol/m³]
    """
    j = states[:j]
    c_s_avg = states[:c_s_avg]

    ∂c_s_avg = ∂states[:c_s_avg]

    res_c_s_avg = res[:c_s_avg]

    # First, retreive the diffusion coefficients
    D_sp_eff, D_sn_eff = coeff_solid_diffusion_effective(states, p)

    function rhs_func(c_s_avg,j,Rp,D_s_eff,N,N_r)
        # Matrices needed for first and second order derivatives
        deriv = derivative_matrices_first_and_second_order(N_r)

        rhsCs = zeros(eltype(j), N*N_r)
        @inbounds for i in 1:N
            ind = (1:N_r) .+ N_r*(i-1)
            c_s = @inbounds @views c_s_avg[ind]

            # First order derivatives matrix multiplication
            ∂ᵣc_s = deriv[1](c_s)

            # Boundary condition at r = 1
            ∂ᵣc_s[N_r] = -j[i]/D_s_eff[i]*Rp

            # Boundary condition at r = 0
            ∂ᵣc_s[1] = 0

            # Second order derivatives matrix multiplication
            ∂ᵣᵣc_s = deriv[2](c_s)

            # Neumann BC at r = 1
            ∂ᵣᵣc_s[end] += 50deriv[2].Δx*∂ᵣc_s[N_r]*(deriv[2].coeff)

            # Make the RHS vector for this particle
            @inbounds rhsCs[ind] .= (D_s_eff[i]./Rp^2) .* [
                3∂ᵣᵣc_s[1]
                ∂ᵣᵣc_s[2:end]+2.0./range(1/(N_r-1), 1, length=N_r-1).*∂ᵣc_s[2:end]
                ]

        end
        return rhsCs
    end

    rhsCs_p = rhs_func(c_s_avg.p,j.p,p.θ[:Rp_p],D_sp_eff,p.N.p,p.N.r_p)
    rhsCs_n = rhs_func(c_s_avg.n,j.n,p.θ[:Rp_n],D_sn_eff,p.N.n,p.N.r_n)

    res_c_s_avg .= [rhsCs_p; rhsCs_n] .- ∂c_s_avg

    return nothing
end
function residuals_c_s_avg!(res, states, ∂states, p::T) where {jac,temp,T<:AbstractModel{jac,temp,:Fickian,:spectral}}
    """
    BETA: Calculate the volume-averaged solid particle concentration residuals using a spectral method [mol/m³]
    """
    j = states[:j]
    c_s_avg = states[:c_s_avg]

    ∂c_s_avg = ∂states[:c_s_avg]

    res_c_s_avg = res[:c_s_avg]

    function cheb(N::Int64)
        x = cos.(range(0,π,length=N+1))
        c = [2.0; ones(N-1); 2.0].*(-1).^(0:N)

        X = repeat(x,1,N+1)
        dX = X .- X'
        D = (c*(1.0./c)')./(dX .+ LinearAlgebra.I(N+1)) # off-diagonal entries
        D[diagind(D)] .-= sum(D', dims=1)[:]

        return D, x
    end

    # First, retreive the diffusion coefficients
    D_sp_eff, D_sn_eff = coeff_solid_diffusion_effective(states, p)

    function rhs_func(c_s_avg,j,Rp,D_s_eff,N,N_r)
        diffusion_matrix, radial_position = cheb(N_r-1)

        rhsCs = zeros(eltype(j),N*N_r)
        @inbounds for i in 1:N
            ind = ind = (1:N_r) .+ N_r*(i-1)
            c_s = c_s_avg[ind]

            ∂ᵣc_s = diffusion_matrix*reverse(c_s)
            ∂ᵣc_s[1] = -j[i]*Rp*0.5/D_s_eff[i] # modified BC value due to cheb scheme
            ∂ᵣc_s[end] = 0 # no flux BC

            rhs_numerator = reverse(diffusion_matrix*(4*D_s_eff[i]*((radial_position .+ 1).^2).*∂ᵣc_s/(Rp^2)))

            rhs_limit_vector = (4*D_s_eff[i]/Rp^2)*3*(diffusion_matrix*∂ᵣc_s) # limit at r_tilde tends to -1 (at center)

            @inbounds rhsCs[ind] .= [
                rhs_limit_vector[end] # L'hopital's rule at the center of the particle
                rhs_numerator[2:end]./((reverse(radial_position[1:end-1]) .+ 1).^2)
                ]
        end
        return rhsCs
    end

    rhsCs_p = rhs_func(c_s_avg.p,j.p,p.θ[:Rp_p],D_sp_eff,p.N.p,p.N.r_p)
    rhsCs_n = rhs_func(c_s_avg.n,j.n,p.θ[:Rp_n],D_sn_eff,p.N.n,p.N.r_n)

    res_c_s_avg .= [rhsCs_p; rhsCs_n] .- ∂c_s_avg

    return nothing
end

function residuals_Q!(res, states, ∂states, p::AbstractModelSolidDiff{:polynomial})
    """
    residuals_Q! is used to implement the three parameters reduced sol for solid phase diffusion [mol/m⁴]
    This sol has been taken from the paper "Efficient Macro-Micro Scale Coupled
    Modeling of Batteries" by Subramanian et al.
    """

    Q = states[:Q]
    j = states[:j]

    ∂Q = ∂states[:Q]

    res_Q = res[:Q]

    # Diffusion coefficients for the solid phase
    D_sp_eff, D_sn_eff = coeff_solid_diffusion_effective(states, p)

    rhsQ_p = @. (-30D_sp_eff*Q.p - 45/2*j.p)/p.θ[:Rp_p]^2
    rhsQ_n = @. (-30D_sn_eff*Q.n - 45/2*j.n)/p.θ[:Rp_n]^2

    res_Q .= [rhsQ_p; rhsQ_n] .- ∂Q

    return nothing
end

function residuals_film!(res, states, ∂states, p::AbstractModel)
    """
    residuals_film! describes the dynamics of the solid-electrolyte layer at the anode side [m]
    """

    j_s = states[:j_s]

    ∂film = ∂states[:film]

    res_film = res[:film]

    rhs_film = -j_s.*p.θ[:M_n]./p.θ[:ρ_n]

    res_film .= rhs_film .- ∂film

    return nothing
end

function residuals_SOH!(res, states, ∂states, p::AbstractModel)
    """
    residuals_SOH! integrates the SOH when aging is enabled and there are losses
    """

    j_s = states[:j_s]

    ∂SOH = ∂states[:SOH][1]

    res_SOH = res[:SOH]

    j_s_int = -trapz(extrapolate_section(j_s, p, :n)...)
    j_s_int *= const_Faradays*surface_area_to_volume_ratio(p)[2]/(3600*calc_I1C(p.θ))

    rhs_SOH = -j_s_int

    res_SOH .= rhs_SOH - ∂SOH

    return nothing
end

function residuals_T!(res, states, ∂states, p)
    """
    Calculate the 1D temperature residuals [K]
    """
    T         = states[:T]
    I_density = states[:I][1]
    Q_rev     = states[:Q_rev]
    Q_rxn     = states[:Q_rxn]
    Q_ohm     = states[:Q_ohm]

    ∂T = ∂states[:T]

    res_T = res[:T]

    Δx = Δx_values(p.N)

    T_BC_sx =  p.θ[:h_cell]*(p.θ[:T_amb]-T[1])/(Δx.a*p.θ[:l_a])
    T_BC_dx = -p.θ[:h_cell]*(T[end]-p.θ[:T_amb])/(Δx.z*p.θ[:l_z])

    block_tridiag(N) = spdiagm(
        -1 => ones(eltype(I_density),N-1),
        0 => -[1;2ones(eltype(I_density),N-2);1],
        +1 => ones(eltype(I_density),N-1),
        )

    # Positive current collector
    A_a = p.θ[:λ_a].*block_tridiag(p.N.a)

    # Cathode
    A_p = p.θ[:λ_p].*block_tridiag(p.N.p)

    # Separator
    A_s = p.θ[:λ_s].*block_tridiag(p.N.s)

    # Anode
    A_n = p.θ[:λ_n].*block_tridiag(p.N.n)

    # Negative current collector
    A_z = p.θ[:λ_z].*block_tridiag(p.N.z)

    A_z[end, end-1:end] = [p.θ[:λ_z] -p.θ[:λ_z]]

    # Divide the matrices by (Δx*l)^2

    A_a /= (Δx.a*p.θ[:l_a])^2
    A_p /= (Δx.p*p.θ[:l_p])^2
    A_s /= (Δx.s*p.θ[:l_s])^2
    A_n /= (Δx.n*p.θ[:l_n])^2
    A_z /= (Δx.z*p.θ[:l_z])^2

    A_tot = zeros(eltype(A_a), (p.N.p+p.N.s+p.N.n+p.N.a+p.N.z), (p.N.p+p.N.s+p.N.n+p.N.a+p.N.z))

    ind = 1:p.N.a
    A_tot[ind,ind] = A_a

    ind = (1:p.N.p) .+ ind[end]
    A_tot[ind,ind] = A_p

    ind = (1:p.N.s) .+ ind[end]
    A_tot[ind,ind] = A_s

    ind = (1:p.N.n) .+ ind[end]
    A_tot[ind,ind] = A_n

    ind = (1:p.N.z) .+ ind[end]
    A_tot[ind,ind] = A_z

    ## Interfaces

    # Interface between aluminium current collector & cathode. We
    # are in the last volume of the current collector
    β_a_p   = (Δx.a*p.θ[:l_a]/2)/(Δx.a*p.θ[:l_a]/2+Δx.p*p.θ[:l_p]/2)
    λ_a_p   = harmonic_mean(β_a_p,p.θ[:λ_a],p.θ[:λ_p])

    den_a_p = Δx.p*p.θ[:l_p]/2 +Δx.a*p.θ[:l_a]/2
    last_a  = p.θ[:λ_a] / (Δx.a*p.θ[:l_a])
    first_p = λ_a_p/den_a_p

    A_tot[p.N.a,p.N.a-1:p.N.a+1] .= [last_a; -(last_a+first_p); first_p]/(Δx.a*p.θ[:l_a])

    # Interface between aluminium current collector & cathode. We
    # are in the first volume of the cathode

    den_a_p  = Δx.p*p.θ[:l_p]/2 +Δx.a*p.θ[:l_a]/2
    second_p = p.θ[:λ_p] / (Δx.p*p.θ[:l_p])
    first_p  = λ_a_p/den_a_p

    A_tot[p.N.a+1,p.N.a:p.N.a+2] .= [first_p; -(second_p+first_p); second_p]/(Δx.p*p.θ[:l_p])

    # Interface between cathode & separator. We
    # are in the last volume of the cathode
    β_p_s = (Δx.p*p.θ[:l_p]/2)/(Δx.s*p.θ[:l_s]/2+Δx.p*p.θ[:l_p]/2)
    λ_p_s = harmonic_mean(β_p_s,p.θ[:λ_p],p.θ[:λ_s])

    den_p_s = Δx.p*p.θ[:l_p]/2 +Δx.s*p.θ[:l_s]/2
    last_p  = p.θ[:λ_p] / (Δx.p*p.θ[:l_p])
    first_s = λ_p_s/den_p_s

    A_tot[p.N.a+p.N.p,p.N.a+p.N.p-1:p.N.a+p.N.p+1] .= [last_p; -(last_p+first_s); first_s]/(Δx.p*p.θ[:l_p])

    # Interface between cathode & separator. We
    # are in the first volume of the separator
    den_p_s  = Δx.p*p.θ[:l_p]/2 +Δx.s*p.θ[:l_s]/2
    second_s = p.θ[:λ_s] / (Δx.s*p.θ[:l_s])
    first_s  = λ_p_s/den_p_s

    A_tot[p.N.a+p.N.p+1,p.N.a+p.N.p:p.N.a+p.N.p+2] .= [first_s; -(second_s+first_s); second_s]/(Δx.s*p.θ[:l_s])

    # Interface between separator anode. We
    # are in the last volume of the separator
    β_s_n = (Δx.s*p.θ[:l_s]/2)/(Δx.s*p.θ[:l_s]/2+Δx.n*p.θ[:l_n]/2)
    λ_s_n = harmonic_mean(β_s_n,p.θ[:λ_s],p.θ[:λ_n])

    den_s_n = Δx.n*p.θ[:l_n]/2 +Δx.s*p.θ[:l_s]/2
    last_s  = p.θ[:λ_s] / (Δx.s*p.θ[:l_s])
    first_n = λ_s_n/den_s_n

    A_tot[p.N.a+p.N.p+p.N.s,p.N.a+p.N.p+p.N.s-1:p.N.a+p.N.p+p.N.s+1] .= [last_s; -(last_s+first_n); first_n]/(Δx.s*p.θ[:l_s])

    # Interface between separator anode. We
    # are in the first volume of the anode

    den_s_n  = Δx.n*p.θ[:l_n]/2 +Δx.s*p.θ[:l_s]/2
    second_n = p.θ[:λ_n] / (Δx.n*p.θ[:l_n])
    first_n  = λ_s_n/den_s_n

    A_tot[p.N.a+p.N.p+p.N.s+1,p.N.a+p.N.p+p.N.s:p.N.a+p.N.p+p.N.s+2] .= [first_n; -(first_n+second_n); second_n]/(Δx.n*p.θ[:l_n])

    # Interface between anode & negative current collector. We
    # are in the last volume of the anode
    β_n_z = (Δx.n*p.θ[:l_n]/2)/(Δx.z*p.θ[:l_z]/2+Δx.n*p.θ[:l_n]/2)
    λ_n_z = harmonic_mean(β_n_z,p.θ[:λ_n],p.θ[:λ_z])

    den_n_co = Δx.n*p.θ[:l_n]/2 +Δx.z*p.θ[:l_z]/2
    last_n   = p.θ[:λ_n] / (Δx.n*p.θ[:l_n])
    first_co = λ_n_z/den_n_co

    A_tot[p.N.a+p.N.p+p.N.s+p.N.n,p.N.a+p.N.p+p.N.s+p.N.n-1:p.N.a+p.N.p+p.N.s+p.N.n+1] .= [last_n; -(last_n+first_co); first_co]/(Δx.n*p.θ[:l_n])

    # Interface between anode & negative current collector. We
    # are in the first volume of the negative current collector

    den_n_co  = Δx.n*p.θ[:l_n]/2 +Δx.z*p.θ[:l_z]/2
    second_co = p.θ[:λ_z] / (Δx.z*p.θ[:l_z])
    first_co  = λ_n_z/den_n_co

    A_tot[p.N.a+p.N.p+p.N.s+p.N.n+1,p.N.a+p.N.p+p.N.s+p.N.n:p.N.a+p.N.p+p.N.s+p.N.n+2] .= [first_co; -(second_co+first_co); second_co]/(Δx.z*p.θ[:l_z])

    Q_rev_tot = [
        zeros(p.N.a)
        Q_rev.p
        zeros(p.N.s)
        Q_rev.n
        zeros(p.N.z)
        ]

    Q_rxn_tot = [
        zeros(p.N.a)
        Q_rxn.p
        zeros(p.N.s)
        Q_rxn.n
        zeros(p.N.z)
        ]

    Q_ohm_tot = [
        (I_density^2)./repeat([p.θ[:σ_a]], p.N.a)
        Q_ohm.p
        Q_ohm.s
        Q_ohm.n
        (I_density^2)./repeat([p.θ[:σ_z]], p.N.z)
        ]

    BC = [
        T_BC_sx
        zeros(length(T)-2)
        T_BC_dx
        ]

    ρ_Cp = [
        repeat([p.θ[:ρ_a]*p.θ[:Cp_a]], p.N.a)
        repeat([p.θ[:ρ_p]*p.θ[:Cp_p]], p.N.p)
        repeat([p.θ[:ρ_s]*p.θ[:Cp_s]], p.N.s)
        repeat([p.θ[:ρ_n]*p.θ[:Cp_n]], p.N.n)
        repeat([p.θ[:ρ_z]*p.θ[:Cp_z]], p.N.z)
        ]

    rhsT   = A_tot*T
    rhsT .+= Q_rev_tot
    rhsT .+= Q_rxn_tot
    rhsT .+= Q_ohm_tot
    rhsT .+= BC
    rhsT ./= ρ_Cp

    res_T .= rhsT .- ∂T

    return nothing
end

function residuals_j!(res, states, p::AbstractModel)
    """
    Calculate the molar flux density of Li-ions residuals at the electrode-electrolyte interface [mol/(m²•s)]
    """
    c_s_star = states[:c_s_star]
    c_e = states[:c_e]
    T = states[:T]
    j = states[:j]

    η = states[:η]

    res_j = res[:j]

    # Calculate the reaction rates
    k_p_eff, k_n_eff = coeff_reaction_rate(states, p)

    # only use for MHC reaction rate
    λ_MHC_p = haskey(p.θ, :λ_MHC_p) ? p.θ[:λ_MHC_p] : 0
    λ_MHC_n = haskey(p.θ, :λ_MHC_n) ? p.θ[:λ_MHC_n] : 0

    j_p_calc = p.numerics.rxn_p(c_s_star.p, c_e.p, T.p, η.p, k_p_eff, λ_MHC_p, p.θ[:c_max_p], p)
    j_n_calc = p.numerics.rxn_n(c_s_star.n, c_e.n, T.n, η.n, k_n_eff, λ_MHC_n, p.θ[:c_max_n], p)

    res_j .= [j_p_calc; j_n_calc] .- j

    return nothing
end

function residuals_j_s!(res, states, p::AbstractModel)
    """
    Calculate the molar flux density side reaction residuals due to SEI resistance [mol/(m²•s)]
    """
    j_s = states[:j_s]
    Φ_s = states[:Φ_s]
    Φ_e = states[:Φ_e]
    film = states[:film]
    T = states[:T]
    I_density = states[:I][1]

    res_j_s = res[:j_s]

    F = const_Faradays
    R = const_Ideal_Gas

    I1C = calc_I1C(p)

    R_film = (p.θ[:R_SEI] .+ film./p.θ[:k_n_aging])

    η_s = Φ_s.n .- Φ_e.n .- p.θ[:Uref_s] .- F.*j_s.*R_film
    α = 0.5

    j_s_calc = -abs.((p.θ[:i_0_jside].*(I_density/I1C)^p.θ[:w]./F).*(-exp.(-α.*F./(R.*T.n).*η_s)))

    # Only activate the side reaction during charge
    j_s_calc .= [IfElse.ifelse(I_density > 0, x, 0) for x in j_s_calc]

    # side reaction residuals
    res_j_s .= j_s .- j_s_calc

    return nothing
end

function residuals_Φ_e!(res, states, p::AbstractModel)
    """
    residuals_Φ_e! evaluates residuals for the electrolyte potential equation discretized using method of lines, [V]
    """

    j   = states[:j_total]
    Φ_e = states[:Φ_e]
    c_e = states[:c_e]
    T   = states[:T]

    K_eff = states[:K_eff]

    res_Φ_e = res[:Φ_e]

    Δx = Δx_values(p.N)
    a_p, a_n = surface_area_to_volume_ratio(p)

    R = const_Ideal_Gas
    F = const_Faradays

    # Interpolate the K_eff values to the edge of the control volume
    K̂_eff_p, K̂_eff_s, K̂_eff_n = interpolate_electrolyte_conductivities(K_eff.p, K_eff.s, K_eff.n, p)

    A_tot = block_matrix_maker(p, K̂_eff_p, K̂_eff_s, K̂_eff_n)

    ind = 1:p.N.p
    A_tot[ind,ind] ./= (Δx.p*p.θ[:l_p])
    ind = (1:p.N.s) .+ ind[end]
    A_tot[ind,ind] ./= (Δx.s*p.θ[:l_s])
    ind = (1:p.N.n) .+ ind[end]
    A_tot[ind,ind] ./= (Δx.n*p.θ[:l_n])

    # Φ_e(x = L) = 0
    A_tot[end, end-1:end] .= [0.0, 1.0]

    ## Interfaces Cathode [last volume of the positive]

    # Here we are in the last volume of the positive
    den = (Δx.p*p.θ[:l_p]/2+Δx.s*p.θ[:l_s]/2)
    last_p = K̂_eff_p[end-1]/(Δx.p*p.θ[:l_p])
    A_tot[p.N.p,p.N.p-1:p.N.p+1] .= [-last_p, (last_p+K̂_eff_p[end]/den), -K̂_eff_p[end]/den]

    ## Interfaces Cathode [first volume of the separator]

    # Here we are in the first volume of the separator
    first_s = K̂_eff_s[1]/(Δx.s*p.θ[:l_s])
    A_tot[p.N.p+1,p.N.p:p.N.p+2] .= [-K̂_eff_p[end]/den, (first_s+K̂_eff_p[end]/den), -first_s]

    ## Interfaces Cathode [last volume of the separator]
    # Here we are in the last volume of the separator
    den = (Δx.n*p.θ[:l_n]/2+Δx.s*p.θ[:l_s]/2)
    last_s = K̂_eff_s[end-1]/(Δx.s*p.θ[:l_s])
    A_tot[p.N.p+p.N.s,p.N.p+p.N.s-1:p.N.p+p.N.s+1] .= [-last_s, (last_s+K̂_eff_s[end]/den), -K̂_eff_s[end]/den]

    ## Interfaces Cathode [first volume of the negative]
    # Here we are inside the first volume of the anode
    first_n = K̂_eff_n[1]/(Δx.n*p.θ[:l_n])
    A_tot[p.N.p+p.N.s+1,p.N.p+p.N.s:p.N.p+p.N.s+2] .= [-K̂_eff_s[end]/den, (first_n+K̂_eff_s[end]/den), -first_n]

    ## Electrolyte concentration interpolation
    # Evaluate the interpolation of the electrolyte concentration values at the
    # edges of the control volumes.
    c̄_e_p, c̄_e_s, c̄_e_n = PETLION.interpolate_electrolyte_concentration(c_e, p)
    ## Temperature interpolation
    # Evaluate the temperature value at the edges of the control volumes
    T̄_p, T̄_s, T̄_n = PETLION.interpolate_temperature(T, p)
    ## Electrolyte fluxes
    # Evaluate the interpolation of the electrolyte concentration fluxes at the
    # edges of the control volumes.
    ∂ₓc_e_p, ∂ₓc_e_s, ∂ₓc_e_n = PETLION.interpolate_electrolyte_concetration_fluxes(c_e, p)

    ## RHS arrays
    ν_p,ν_s,ν_n = p.numerics.thermodynamic_factor(c_e.p, c_e.s, c_e.n, T.p, T.s, T.n, p)
    ν = [ν_p;ν_s;ν_n]

    K = 2R.*(1-p.θ[:t₊]).*ν[1:end-1]/F

    prod_tot = [
        K̂_eff_p.*T̄_p.*∂ₓc_e_p./c̄_e_p # p
        K̂_eff_s.*T̄_s.*∂ₓc_e_s./c̄_e_s # s
        K̂_eff_n[1:end-1].*T̄_n.*∂ₓc_e_n./c̄_e_n # n
    ]

    prod_tot[2:end] .-= prod_tot[1:end-1]

    f = -K.*prod_tot

    ind_p = (1:p.N.p)
    ind_n = (1:p.N.n) .+ (p.N.p+p.N.s)

    f[ind_p] .+= (Δx.p*p.θ[:l_p]*F*a_p)*j.p
    f[ind_n[1:end-1]] .+= (Δx.n*p.θ[:l_n]*F*a_n)*j.n[1:end-1]

    # The boundary condition enforces that Φ_e(x=L) = 0
    append!(f, 0.0)

    # Return the residual value for the electrolyte potential
    res_Φ_e .= A_tot*Φ_e .- f

    return nothing
end

function residuals_Φ_s!(res, states, p::AbstractModel)
    """
    residuals_Φ_s! evaluates the residuals of the solid potential equation [V]
    """

    j = states[:j_total]
    Φ_s = states[:Φ_s]
    I_density = states[:I][1]

    res_Φ_s = res[:Φ_s]

    F = const_Faradays
    Δx = Δx_values(p.N)

    # Effective conductivities and surface area to volume ratio
    σ_eff_p, σ_eff_n = conductivity_effective(p)
    a_p, a_n = surface_area_to_volume_ratio(p)

    ## Cathode

    # RHS for the solid potential in the cathode and anode
    f_p = @. p.θ[:l_p]^2*Δx.p^2*a_p*F*j.p/σ_eff_p
    f_n = @. p.θ[:l_n]^2*Δx.n^2*a_n*F*j.n/σ_eff_n

    # Additional term at the electrode-current collector interface
    f_p[1]   += -I_density*(Δx.p*p.θ[:l_p]/σ_eff_p)
    f_n[end] += +I_density*(Δx.n*p.θ[:l_n]/σ_eff_n)

    block_tridiag(N) = spdiagm(
        -1 => ones(eltype(j.p),N-1),
        0 => -[1;2ones(eltype(j.p),N-2);1],
        +1 => ones(eltype(j.p),N-1),
        )

    A_p = block_tridiag(p.N.p)
    A_n = block_tridiag(p.N.n)

    A_n *= Φ_s.n
    A_n .-= f_n

    A_p *= Φ_s.p
    A_p .-= f_p

    res_Φ_s .= [A_p; A_n]

    return nothing
end

function residuals_scalar!(res, states, p::AbstractModel)
    """
    *** THE RESIDUALS ARE HANDLED IN `scalar_residual.jl` ***

    This function sets an arbitrary value for the current residuals [C-rate]
    """

    res_I = res[:I]

    res_I .= 0.0

    return nothing
end


text:
#RSDeltaSigma: designLCBP() algorithm
#-------------------------------------------------------------------------------

#==Types
===============================================================================#
"""`LCBPParameters`

A struct containing the (`n`, `OSR`, `Hinf`, `f0`, `t` and `form`) arguments
for the LCBP filter, plus the following fields:
 - `t`: A two-element vector containing the start and end times of the
   feedback waveform. Default `= [0 1]`.
 - `l`: The inductances in each tank.
 - `c`: The capacitances in each tank.
 - `gu`: The transconductance from the u input to each of the tanks.
   The final term is the voltage gain from u to the comparator input.
   Default `= [1 0 .. ]`. (1 by n+1)
 - `gv`: The transconductance from the v output to each of the tanks. (1 by n+1)
 - `gw`: The inter-stage transconductances. Default `= [1 1 .. ]`. (1 by n-1)
 - `gx`: The gains from the top of each tank resistor to the comparator.
   Default `= [0 0 .. 1]`. (1 by n)
 - `rx`: The resistances inserted between the output of each interstage
   transconductance and top of each tank. Default `= [0 ..]`. (1 by n)
   Note that rx[1] is not used.
 - `rl`: The series resistance of each inductor. Default `= [0 ..]`. (1 by n)
 - `gc`: The parallel conductances of each tank. Default `= [0 ..]`. (1 by n)
"""
mutable struct LCBPParameters
	n::Int
	OSR::Int
	Hinf::Float64
	f0::Float64
	t
	form::Symbol
	l::Array{Float64}; c::Array{Float64}
	gu::Array{Float64}; gv::Array{Float64}; gw::Array{Float64}; gx::Array{Float64}
	rx::Array{Float64}; rl::Array{Float64}; gc::Array{Float64}
end
LCBPParameters(n=3; OSR=64, Hinf=1.6, f0=0.25, t=[0 1], form::Symbol=:FB) =
	LCBPParameters(n, OSR, Hinf, f0, t, form, [], [], [], [], [], [], [], [], [])

#==
===============================================================================#
"""`(H,L0,ABCD,k)=LCparam2tf(param,k=1)`

Compute the transfer functions (and the Q. gain, k) of an LC modulator.

# Arguments
 - `param`: A LCBPParameters struct.
 - `k`: The value to assume for the quantizer gain.

# Note:
 - If `k` is omitted or if `k` is defaulted, `k` is taken to be 1.
 - If `k` is `0`, a value is computed by the formula
   `k = mean(abs(y))/mean(y.^2)`, where y is the quantizer input
   sequence when the modulator is fed with a small input at f0.

# Output
 - `H`: The closed-loop noise transfer function (in z).
 - `L0`: The open-loop tf (in s) from u to v.  G  may be evaluated
   using `evalContSTF(L0,H,f)`.
 - `k`: The value of the quantizer gain used to compute the above tfs.
 - `ABCD`: The state-space description of the discrete-time system

For the conversion from continuous to discrete time, corrections to the
standard MATLAB formulae need to be applied in order for the sampled-data
value to be taken at the END of the time interval.

This effectively makes t1>0 and requires corrections if t2>=1 and gv(n+1)~=0.
More complex formulae are needed for t2>=1, so for the moment this code
requires t2<1.
"""
function LCparam2tf(param::LCBPParameters, k::Float64=1.0)
	l = param.l; c = param.c
	f0 = mean(1 / sqrt.(l .* c)) / (2π *pi)
	n = length(l)

	gu = param.gu
		isempty(gu) && (gu = vcat(1, zeros(n))) #Defaults, if needed
	gu = padr(gu, n+1, 0) #Pad out with zeros, if needed.
	gv = padr(param.gv, n+1, 0) #Pad out with zeros, if needed.
	gw = param.gw
		isempty(gw) && (gw = ones(n-1)) #Defaults, if needed
	gx = param.gx
		isempty(gx) && (gx = vcat(zeros(n-1), 1)) #Defaults, if needed
	rx = param.rx
	if isempty(rx)
		rx = zeros(n)
	elseif 1==length(rx)
		throw("Not sure: Original code did not seem to make sense.")
		rx = ones(n)
	end
	t = param.t
		isempty(t) && (t = [0 1]) #Defaults, if needed
	rl = param.rl
	if isempty(rl)
		rl = zeros(n)
	elseif 1==length(rl)
		throw("Not sure: Original code did not seem to make sense.")
		rl = ones(n)
	end
	gc = param.gc
	if isempty(gc)
		gc = zeros(n)
	elseif 1==length(gc)
		throw("Not sure: Original code did not seem to make sense.")
		gc = ones(n)
	end

	#Form the state-space description of the modulator
	#The states are ordered v1, i1, v2, i2 ...
	n2 = 2*n
	local y
	ABc = zeros(n2,n2+2); CDc = zeros(1,n2+2)
	for i = 1:n
		i2 = 2*i-1
		#y represents the voltage at the top of the resistor stacked on the tank.
		if i == 1
			y = vcat(1, zeros(n2+1))
		else
			ABc[i2,:] = gw[i-1]*y
			y = rx[i]*gw[i-1]*y; y[i2] = 1
		end
		ABc[i2,n2+1:n2+2] = [gu[i] -gv[i]]
		ABc[i2:i2+1,i2:i2+1] = [-gc[i] -1; 1 -rl[i]]
		ABc[i2,:] = ABc[i2,:]/c[i]
		ABc[i2+1,:] = ABc[i2+1,:]/l[i]
		CDc = CDc + (gx[i]*y)'
	end
	CDc = CDc + [zeros(1,n2) gu[n+1] -gv[n+1]]

	#THIS IS THE NEW CODE !!
	(Ac,Bc,Cc,Dc) = partitionABCD([ABc;CDc], 2)
	sys_c = _ss( Ac, Bc[:,2], Cc, Dc[2] )
	(sys,) = mapCtoD(sys_c, t=t, f0=f0)
	#Augment the input matrix.
	A=sys.A
	B=[padb(Bc[:,1],size(sys.B,1)) sys.B]
	C=sys.C
	D=[Dc[1] sys.D]

	#Compute L0; use the LC parameters to compute the poles and thereby
	#ensure exact agreement with the zeros in H.
	#!! THIS IS GOING TO NEED FIXING
	s_poles = zeros(2*n) .+ 0j #Make complex
	for i in 1:n
		_v = [l[i]*c[i], rl[i]*c[i] + gc[i]*l[i], 1+rl[i]*gc[i]]
		s_poles[2*i-1:2*i] = _roots(_v)
	end
	LF0 = _ss(Ac, Bc[:,1], Cc, Dc[1])
	L0 = _zpk( LF0 )
	L0.p = s_poles

	#Compute H. Use the poles in L0 to compute the zeros in H.
	ABCD =[A B; C D]
	if k==0 #Estimate k by simulatiing the modulator with a small input
		w0 = mean(1/sqrt.(l .* c))
		H = calculateTF(ABCD)
		H.z = exp.(s_poles)
		stf0 = abs(evalTFP(L0, H, w0/(2π)))
		u = 0.1/stf0*sin.(w0*(0:10000))
		y = simulateDSM(u, [A B; C D]).y
		k = mean(abs.(y))/mean(y .^ 2)
	end
	(H, STF) = calculateTF(ABCD, [k])
	H.z = exp.(s_poles)

	#Correct L0k to include the quantizer gain
	L0.k = L0.k*k
	return (H,L0,ABCD,k)
end

#Dead code???
function XLCoptparam2tf(x, param)
	n = param.n
	#Uncomment the following line to take into account the quantizer gain
	#before returning results and doing the plots
	#(H, L0, ABCD, k) = LCparam2tf(param, 0)
	#Uncomment the following lines to yield parameters corresponding to k=1
	#param.gu = k*param.gu; param.gv = k*param.gv
	#ABCD[2*n+1, :] = ABCD[2*n+1, :] / k
	#ABCD[:, 2*n+[1 2]] = ABCD[:, 2*n+[1 2]] * k

	#Now fix up the inputs for 0dB gain at f0.
	gain_f0 = abs(evalTFP(L0, H, f0))
	param.gu = param.gu/gain_f0; L0.k = L0.k/gain_f0
	ABCD[:,2*n+1] = ABCD[:,2*n+1]/gain_f0
	if (:FF==form) || (:FFR==form)
		param.gu[n+1] = 1 #For the flattest STF
	end
	if dbg
		@warn("LCplotTF: not yet implemented")
		#LCplotTF(H,L0,param)
	end
	return (H, L0, ABCD, param)
end

"""`LCoptparam2tf(...)`

Convert optimization parameters to transfer functions.
"""
function LCoptparam2tf(x,param)
	n = param.n
	param.gw = ones(n-1)

	if :FB == param.form
		param.gu = vcat(1, zeros(n))
		param.gv = vcat(x)
		param.gx = vcat(zeros(n-1), 1)
		param.rx = zeros(n)
	elseif :FF == param.form
		param.gu = vcat(1, zeros(n-1), 1)
		param.gv = vcat(1, zeros(n))
		param.gx = vcat(x)
		param.rx = zeros(n)
	elseif :FFR == param.form
		param.gu = vcat(1, zeros(n-1), 1)
		param.gv = vcat(x[1], zeros(n))
		param.gx = vcat(zeros(n-1), 1)
		param.rx = vcat(0, x[2:n])
	elseif :GEN == param.form
		param.gv = x[1:n]
		param.rx = vcat(0, x[n+1:2*n-1])
	else
		form = param.form
		throw("form=:$form is not supported.")
	end

	(H, L0, ABCD) = LCparam2tf(param)
	return (H,L0,ABCD,param)
end

"""`f = LCObj1a(x, param, max_radius, dbg)`

The objective function for the initial optimization process used to put
the roots of the denominator of the LCBP NTF inside the unit circle.
"""
function LCObj1a(x, param, max_radius, dbg)
	return f=1 #No actual objective
end

"""`(C, Ceq) = LCObj1b(x, param, max_radius, dbg)`

The constraints function for the initial optimization process used to put
the roots of the denominator of the LCBP NTF inside the unit circle.
"""
function LCObj1b(x, param, max_radius, dbg)
	(H,) = LCoptparam2tf(x,param)
	objective = 1 #No actual objective

	rmax = maximum(abs.(H.p[:]))
	C = rmax - max_radius
	Ceq = []
	if dbg
		@show(x)
		@show(rmax)
	end
	return (C, Ceq)
end

#The objective function
function LCObj2a(x, param, rmax, dbg)
	if dbg
		@show(x)
	end
	(H, L0, ABCD) = LCoptparam2tf(x, param)
	#The objective is the in-band noise, in dB
	f1 = param.f0 - 0.25/param.OSR
	f2 = param.f0 + 0.25/param.OSR
	objective = dbv(rmsGain(H,f1,f2)/sqrt(3*param.OSR))
	println(@sprintf("N0^2 = %.1fdB", objective))
	return objective
end

"""`(C, Ceq) = LCObj2b(x,param,rmax,dbg)`
% The constraints: r-rmax, Hinf-Hinfdesired, quantizer input,
% and for the FB form, |stf|-1.1.
"""
function LCObj2b(x, param, rmax, dbg)
	(H, L0, ABCD) = LCoptparam2tf(x, param)
	Ceq = []
	local C

	#The constraints are on the maximum radius of the poles of the NTF
	#the infinity-norm of the NTF, the peaking of the STF (for the 'FB' form),
	#and the maximum input to the quantizer.
	max_radius = maximum(abs(H.p[:]))
	H_inf = infnorm(H)
	stf0 = abs( evalTFP(L0, H, param.f0) )
	y = simulateDSM(0.5/stf0*sin(2π*param.f0*(0:1000)), ABCD).y
	ymax = maximum(abs.(y))
	if :FB == param.form
		f1 = param.f0 - 0.25/param.OSR
		f2 = param.f0 + 0.25/param.OSR
		stf1 = abs( evalTFP(L0, H, f1) )/stf0
		stf2 = abs( evalTFP(L0, H, f2) )/stf0
		C = vcat(100*(max_radius .- rmax), H_inf-param.Hinf, stf1-1.01, stf2-1.01, (ymax-5)/10)
		#C = vcat(100*(max_radius .- rmax), H_inf-param.Hinf, stf1-1.01, stf2-1.01, log(ymax/5))
	else
		C = vcat(100*(max_radius .- rmax), H_inf-param.Hinf, (ymax-5)/10)
	end

	if dbg
		println(@sprintf("constraints = %8.5f", C))
				println(@sprintf("rmax = %5.3f, Hinf = %4.2f, ymax = %5.3f\n",
					max_radius, H_inf, ymax)
				)
	end
	return (C, Ceq)
end

#==Main algorithm
===============================================================================#
"""`(param,H,L0,ABCD,x)=designLCBP(n=3,OSR=64,opt=2,Hinf=1.6,f0=1/4,t=[0 1],form='FB',x0,dbg=0)`

Design an LC modulator of order `2*n` with center frequency `f0`.
The feedback waveform is a rectangular pulse from `t[1]` to `t[2]`.

# Arguments
 - `n`: The number of LC tanks.
 - `OSR`: The oversampling ratio.
 - `opt`: A flag for NTF zero optimization. See synthesizeNTF.
 - `Hinf`: The target out-of-band gain of the NTF.
 - `f0`: The center frequency of the modulator, normalized to fs.
 - `t`: A 2-element vector containing the (normalized) feedback pulse edge times.
 - `form`: The modulator topology. One of `:FB`, `:FF`, `:FFR` or `:GEN`.
 - `x0`: An initial guess at the parameter vector for the optimizer.
 - `dbg`: A flag for enabling debugging/progress report output.

# Output
 - `param`: A LCBPParameters struct.
 - `H`: The noise transfer function of the modulator, taking into account
   the post-design quantizer gain.
 - `L0`: A description of the open-loop TF from u to v, for use with evalTFP().
 - `ABCD`: The state-space representation of the discrete-time equivalent.
"""
function designLCBP(n=3; OSR=64, opt=2, Hinf=1.6, f0=0.25, t=[0 1],
		form::Symbol=:FB, x0=NaN, dbg::Bool=false)
	param = LCBPParameters(n, OSR=OSR, Hinf=Hinf, f0=f0, t=t, form=form)
	#Compute the resonant frequencies for the LC tanks
	#and the L and C values.
	local w
	if opt>0
		w = 2π*(f0 .+ 0.25/OSR*ds_optzeros(n,opt))
	else
		w = 2π*f0*ones(n)
	end
	param.l = (1 ./ w); param.c = (1 ./ w) #Impedance scaling is possible.
	#Find values for the parameters that yield a stable system.
	#?? Could some control theory help here ??
	x = x0
	if isnan(x0)
		x = collect(2:n+1)*1.0 #Needs to be Float64
	elseif length(x0) != n
		throw("Initial guess (x0) has the wrong number of parameters.")
	end
	if dbg
		println("Iterations for initial pole placement:")
	end
	max_radius = 0.97
	opt = Optim.Options(x_tol=0.001, f_tol=1, iterations=1000)
	lx = Float64[0,0,0]; ux = Float64[10,10,10]
	lc = [-Inf]; uc = [max_radius^2]
	fun1!(x) = LCObj1a(x, param, max_radius, dbg)
	con_c1!(c, x) = ((C, Ceq)=LCObj1b(x, param, max_radius, dbg); C)
#	dfc = TwiceDifferentiableConstraints(con_c!, con_jacobian!, con_h!,lx,ux,lc,uc)
	dfc = TwiceDifferentiableConstraints(con_c1!,lx,ux,lc,uc)#,:forward)
#	x = optimize(fun1!, dfc, x, IPNewton()).minimum
	x = optimize(fun1!, con_c1!, x, LBFGS()).minimum
	#x = fmincon(@LCObj1a,x,[],[],[],[],[],[],@LCObj1b,options,param,0.97,dbg);
	(H,) = LCoptparam2tf(x,param)
	rmax = maximum(abs.(H.p[:]))
	if rmax>max_radius #Failure to converge!
		throw("Unable to find parameter values which stabilize the system.")
	end
	#Do the main optimization for the feedback parameters.
	if dbg
		println("\nParameter Optimization:")
	end
	max_radius = 0.97
	opt = Optim.Options(x_tol=0.001, f_tol=1, iterations=1000,
		show_trace=dbg, #Extra debugging info
		g_tol=0.01 #Solver-specific
	)
	lx = Float64[]; ux = Float64[]
	lc = [-Inf]; uc = [max_radius^2]
	fun2!(x) = LCObj2a(x, param, max_radius, dbg)
	con_c2!(c, x) = ((C, Ceq)=LCObj2b(x, param, max_radius, dbg); C)
	dfc = TwiceDifferentiableConstraints(con_c2!,lx,ux,lc,uc,:forward)
	x = optimize(fun2!, dfc, x, IPNewton()).minimum
	#x = fmincon(@LCObj2a,x,[],[],[],[],[],[],@LCObj2b,options,param,0.97,dbg);
	return (param,H,L0,ABCD,x)
end

#Last line

ext: jl


text:
using DbgMacro

function foo()
    m = [1 2; 3 4]
    @dbg 1+2 "Hello" m
    42
end

expected_foo_err = """
Main.TestDbg:test_dbg.jl:5  1 + 2 = 3
Main.TestDbg:test_dbg.jl:5  "Hello" = "Hello"
Main.TestDbg:test_dbg.jl:5  m = [1 2; 3 4]
"""

using Suppressor: @capture_err

err = @capture_err global result = foo()
@test result == 42
@test err == expected_foo_err


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "SPIRV_Cross"
version = v"2021.01.15"

# Collection of sources required to complete build
sources = [
    GitSource("https://github.com/KhronosGroup/SPIRV-Cross.git", "9acb9ec31f5a8ef80ea6b994bb77be787b08d3d1")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/SPIRV-Cross/

install_license LICENSE

CMAKE_FLAGS=()
CMAKE_FLAGS+=(-DCMAKE_BUILD_TYPE=Release)
CMAKE_FLAGS+=(-DCMAKE_INSTALL_PREFIX=${prefix})
CMAKE_FLAGS+=(-DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN})
CMAKE_FLAGS+=(-DCMAKE_CROSSCOMPILING:BOOL=ON)
CMAKE_FLAGS+=(-DBUILD_SHARED_LIBS=ON)
CMAKE_FLAGS+=(-DSPIRV_CROSS_SHARED=ON)

mkdir build && cd build

cmake  ${CMAKE_FLAGS[@]} ..
make -j${nproc}
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = expand_cxxstring_abis(supported_platforms()) # build on all supported platforms

# The products that we will ensure are always built
products = [
    LibraryProduct("libspirv-cross-c-shared", :libspirv_cross),
    ExecutableProduct("spirv-cross", :spirv_cross)
]

# Dependencies that must be installed before this package can be built
dependencies = Dependency[
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; julia_compat = "1.6")

max_stars_repo_name: jakebolewski/Yggdrasil

max_issues_repo_name: jakebolewski/Yggdrasil

max_forks_repo_name: jakebolewski/Yggdrasil


text:
@cdef :PySet_New PyPtr (PyPtr,)
@cdef :PyFrozenSet_New PyPtr (PyPtr,)
@cdef :PySet_Add Cint (PyPtr, PyPtr)

const PySet_Type__ref = Ref(PyNULL)
PySet_Type() = pyglobal(PySet_Type__ref, :PySet_Type)

PySet_Check(o) = Py_TypeCheckFast(o, PySet_Type())

PySet_CheckExact(o) = Py_TypeCheckExact(o, PySet_Type())

const PyFrozenSet_Type__ref = Ref(PyNULL)
PyFrozenSet_Type() = pyglobal(PyFrozenSet_Type__ref, :PyFrozenSet_Type)

PyFrozenSet_Check(o) = Py_TypeCheckFast(o, PyFrozenSet_Type())

PyFrozenSet_CheckExact(o) = Py_TypeCheckExact(o, PyFrozenSet_Type())

PyAnySet_Check(o) = PySet_Check(o) || PyFrozenSet_Check(o)

PyAnySet_CheckExact(o) = PySet_CheckExact(o) || PyFrozenSet_CheckExact(o)

_PySet_FromIter(r::PyPtr, xs) = begin
    try
        for x in xs
            xo = PyObject_From(x)
            isnull(xo) && (Py_DecRef(r); return PyNULL)
            err = PySet_Add(r, xo)
            Py_DecRef(xo)
            ism1(err) && (Py_DecRef(r); return PyNULL)
        end
        return r
    catch err
        Py_DecRef(r)
        PyErr_SetString(PyExc_Exception(), "Julia error: $err")
        return PyNULL
    end
end

PySet_FromIter(xs) = begin
    r = PySet_New(C_NULL)
    isnull(r) && return PyNULL
    _PySet_FromIter(r, xs)
end

PyFrozenSet_FromIter(xs) = begin
    r = PyFrozenSet_New(C_NULL)
    isnull(r) && return PyNULL
    _PySet_FromIter(r, xs)
end

PySet_From(x::Union{AbstractSet,AbstractVector,Tuple}) = PySet_FromIter(x)

PyFrozenSet_From(x::Union{AbstractSet,AbstractVector,Tuple}) = PyFrozenSet_FromIter(x)


text:
module Observables

export Observable, on, off, onany, connect!, obsid, async_latest, throttle

import Base.Iterators.filter

const addhandler_callbacks = []
const removehandler_callbacks = []

abstract type AbstractObservable{T} end

# Internal function that doesn't need invokelatest!
abstract type InternalFunction <: Function end

function observe(::S) where {S<:AbstractObservable}
    error("observe not defined for AbstractObservable $S")
end

"""
Like a `Ref` but updates can be watched by adding a handler using `on`.
"""
mutable struct Observable{T} <: AbstractObservable{T}
    listeners::Vector{Any}
    val::T
    Observable{T}() where {T} = new{T}([])
    Observable{T}(val) where {T} = new{T}([], val)
end

function Base.copy(observable::Observable{T}) where T
    result = Observable{T}(observable[])
    on(observable) do value
        result[] = value
    end
    return result
end

Observable(val::T) where {T} = Observable{T}(val)

observe(x::Observable) = x

Base.convert(::Type{Observable}, x::AbstractObservable) = x
Base.convert(::Type{Observable{T}}, x::AbstractObservable{T}) where {T} = x

function Base.convert(::Type{Observable{T}}, x::AbstractObservable) where {T}
    result = Observable{T}(convert(T, x[]))
    on(x) do value
        result[] = convert(T, value)
    end
    return result
end

Base.convert(::Type{T}, x) where {T<:Observable} = T(x)

function Base.getproperty(obs::Observable, field::Symbol)
    if field === :val
        return getfield(obs, field)
    elseif field === :listeners
        return getfield(obs, field)
    elseif field === :id
        return obsid(obs)
    else
        error("Field $(field) not found")
    end
end

"""
    notify!(observable::AbstractObservable)

Pushes an updates to all listeners of `observable`
"""
function notify!(observable::AbstractObservable)
    observable[] = observable[]
    return
end

function Base.show(io::IO, x::Observable{T}) where T
    println(io, "Observable{$T} with $(length(x.listeners)) listeners. Value:")
    if isdefined(x, :val)
        show(io, x.val)
    else
        println(io, "not assigned yet!")
    end
end

Base.show(io::IO, ::MIME"application/prs.juno.inline", x::Observable) = x

"""
    on(f, observable::AbstractObservable)

Adds function `f` as listener to `observable`. Whenever `observable`'s value
is set via `observable[] = val` `f` is called with `val`.
"""
function on(f, observable::AbstractObservable)
    push!(listeners(observable), f)
    for g in addhandler_callbacks
        g(f, observable)
    end
    return f
end

"""
    off(observable::AbstractObservable, f)

Removes `f` from listeners of `observable`.
"""
function off(observable::AbstractObservable, f)
    callbacks = listeners(observable)
    for (i, f2) in enumerate(callbacks)
        if f === f2
            deleteat!(callbacks, i)
            for g in removehandler_callbacks
                g(observable, f)
            end
            return
        end
    end
    throw(KeyError(f))
end

"""
    observable[] = val

Updates the value of an `Observable` to `val` and call its listeners.
"""
function Base.setindex!(observable::Observable, val; notify=(x)->true)
    observable.val = val
    for f in listeners(observable)
        if notify(f)
            if f isa InternalFunction
                f(val)
            else
                Base.invokelatest(f, val)
            end
        end
    end
end

####################################################
# tasks & channel api

Observable(val::Channel{T}) where {T} = Observable{T}(val)
Observable(val::Task) = Observable{Any}(val)
function Observable{T}(val::Union{Task, Channel}) where {T}
    observable = Observable{T}()
    observable[] = val
    return observable
end

function Base.setindex!(observable::Observable, val_async::Task; notify=x->true)
    return @async begin
        try
            val = fetch(val_async)
            setindex!(observable, val, notify=notify)
        catch e
            Base.showerror(stderr, e)
            Base.show_backtrace(stderr, catch_backtrace())
        end
    end
end

function Base.setindex!(observable::Observable, channel::Channel; notify=x->true)
    return @async begin
        try
            for val in channel
                setindex!(observable, val, notify=notify)
                yield()
            end
        catch e
            Base.showerror(stderr, e)
            Base.show_backtrace(stderr, catch_backtrace())
        end
    end
end

function Base.setindex!(observable::AbstractObservable, val; notify=x->true)
    Base.setindex!(observe(observable), val; notify=notify)
end

function setexcludinghandlers(observable::AbstractObservable, val, pred=x->true)
    setindex!(observable, val; notify=pred)
end

"""
    observable[]

Returns the current value of `observable`.
"""
Base.getindex(observable::Observable) = observable.val

Base.getindex(observable::AbstractObservable) = getindex(observe(observable))

### Utilities

"""
    to_value(x::Union{Any, AbstractObservable})
Extracts the value of an observable, and returns the object if it's not an observable!
"""
to_value(observable::AbstractObservable) = observable[]
to_value(x) = x

"""
    obsid(observable::Observable)

Gets a unique id for an observable!
"""
obsid(observable::Observable) = string(objectid(observable))
obsid(observable::AbstractObservable) = obsid(observe(observable))

listeners(observable::Observable) = observable.listeners
listeners(observable::AbstractObservable) = listeners(observe(observable))

struct OnUpdate{F, Args} <: InternalFunction
    f::F
    args::Args
end

(ou::OnUpdate)(_) = ou.f(map(to_value, ou.args)...)

"""
    onany(f, args...)

Calls `f` on updates to any oservable refs in `args`.
`args` may contain any number of `Observable` objects.
`f` will be passed the values contained in the refs as the respective argument.
All other objects in `args` are passed as-is.
"""
function onany(f, args...)
    callback = OnUpdate(f, args)
    for observable in args
        (observable isa AbstractObservable) && on(callback, observable)
    end
end

struct MapUpdater{F, T} <: InternalFunction
    f::F
    observable::Observable{T}
end

function (mu::MapUpdater)(args...)
    mu.observable[] = mu.f(args...)
end

"""
    map!(f, observable::AbstractObservable, args...)

Updates `observable` with the result of calling `f` with values extracted from args.
`args` may contain any number of `Observable` objects.
`f` will be passed the values contained in the refs as the respective argument.
All other objects in `args` are passed as-is.
"""
function Base.map!(f, observable::AbstractObservable, os...)
    onany(MapUpdater(f, observable), os...)
    return observable
end

"""
    connect!(o1::AbstractObservable, o2::AbstractObservable)

Forwards all updates from `o2` to `o1`
"""
connect!(o1::AbstractObservable, o2::AbstractObservable) = map!(identity, o2, o1)

"""
    map(f, observable::AbstractObservable, args...)

Creates a new observable ref which contains the result of `f` applied to
values extracted from args. The second argument `observable` must be an observable ref for
dispatch reasons. `args` may contain any number of `Observable` objects.
`f` will be passed the values contained in the refs as the respective argument.
All other objects in `args` are passed as-is.
"""
function Base.map(f, observable::AbstractObservable, os...;
                  init=f(observable[], map(to_value, os)...))
    map!(f, Observable{Any}(init), observable, os...)
end

Base.eltype(::AbstractObservable{T}) where {T} = T

"""
    async_latest(observable::AbstractObservable, n=1)

Returns an `Observable` which drops all but
the last `n` updates to `observable` if processing the updates
takes longer than the interval between updates.

This is useful if you want to pass the updates from,
say, a slider to a plotting function that takes a while to
compute. The plot will directly compute the last frame
skipping the intermediate ones.

# Example:
```
observable = Observable(0)
function compute_something(x)
    for i=1:10^8 rand() end # simulate something expensive
    println("updated with \$x")
end
o_latest = async_latest(observable, 1)
on(compute_something, o_latest) # compute something on the latest update

for i=1:5
    observable[] = i
end
```
"""
function async_latest(input::AbstractObservable{T}, n=1) where T
    buffer = T[]
    cond = Condition()
    lck  = ReentrantLock() # advisory lock for access to buffer
    output = Observable{T}(input[]) # output

    @async while true
        while true # while !isempty(buffer) but with a lock
            # transact a pop
            lock(lck)
            if isempty(buffer)
                unlock(lck)
                break
            end
            upd = pop!(buffer)
            unlock(lck)

            output[] = upd
        end
        wait(cond)
    end

    on(input) do val
        lock(lck)
        if length(buffer) < n
            push!(buffer, val)
        else
            while length(buffer) >= n
                pop!(buffer)
            end
            pushfirst!(buffer, val)
        end
        unlock(lck)
        notify(cond)
    end

    output
end

# TODO: overload broadcast on v0.6
include("observablepair.jl")
include("flatten.jl")
include("time.jl")
include("macros.jl")

end # module


text:
using SymbolicUtils, SpecialFunctions
using SymbolicUtils.Rewriters
using Symbolics

Incomplete_beta(z, a, b) = beta(a, b)*beta_inc(a, b, z)[1]

# These are direct computing, we need to change to the numerical methods in http://www.mymathlib.com/c_source/functions/fresnel_sin_cos_integrals/fresnel_auxiliary_sine_integral.c
AuxiliaryFresnelSin(x) = sqrt(2/pi)*quadgk(t->exp(-2x*t)*sin(t^2), 0, Inf)[1]
AuxiliaryFresnelCos(x) = sqrt(2/pi)*quadgk(t->exp(-2x*t)*cos(t^2), 0, Inf)[1]

@register_symbolic Incomplete_beta(z::Number, a::Number, b::Number)
@register_symbolic AuxiliaryFresnelCos(x::Number)
@register_symbolic AuxiliaryFresnelSin(x::Number)
@register_symbolic Hypergeometric1F1(a::Number, b::Number, c::Number, z::Number)
@register_symbolic Struve(p::Number, x::Number)# Here p is the order of Struve function
@register_symbolic MStruve(p::Number, x::Number)# Here p is the order of modified Struve function
@register_symbolic Legendre(p::Number, x::Number)# Here p is the order of Legendre function
@register_symbolic SinIntegral(x::Number)
@register_symbolic HyperSinIntegral(x::Number)

SEMIDIFFRULES = [
        # CONSTANTS AND POWERS
        @acrule sqrt(~x) => sqrt(pi)/2
        #@acrule((~x) => 2*sqrt(~x/pi))
        @acrule (~x)^(~α) => gamma(~α+1)/gamma(~α+1/2)*(~x)^(~α-1/2)

        # BINOMIALS
        @acrule sqrt(1+~x) => 1/sqrt(pi*~x) + atan(sqrt(~x))/sqrt(pi)
        @acrule sqrt(1-~x) => 1/sqrt(pi*~x) - atanh(sqrt(~x))/sqrt(pi)
        @acrule 1/sqrt(1+~x) => 1/(sqrt(pi*~x)*(1+~x))
        @acrule 1/sqrt(1-~x) => 1/(sqrt(pi*~x)*(1-~x))
        @acrule 1/(1+~x) => (sqrt(1+~x)-sqrt(~x)*asinh(sqrt(~x)))/(sqrt(pi*~x)*(1+~x)^(3/2))
        @acrule 1/(1-~x) => (sqrt(1-~x)+sqrt(~x)*asinh(sqrt(~x)))/(sqrt(pi*~x)*(1-~x)^(3/2))
        @acrule 1/(1+~x)^(3/2) => (1-~x)/(sqrt(pi*~x)*(1+~x)^2)
        @acrule 1/(1-~x)^(3/2) => (1+~x)/(sqrt(pi*~x)*(1-~x)^2)
        @acrule 1/(1-~x)^(~p) => -Incomplete_beta(~x, -1/2, ~p+1/2)/(2*sqrt(pi)*(1-~x)^(~p+1/2))
        @acrule 1/sqrt(~x*(1-~x)) => (ellipe(~x)-(1-~x)*ellipk(~x))/(~x*(1-~x)*sqrt(pi))
        @acrule 1/sqrt(~x*(1+~x)) => (ellipe(~x/(1+~x))-ellipk(~x/(1+~x)))/(~x*sqrt(pi*(1+~x)))
        @acrule sqrt(~x/(1-~x)) => ellipe(~x)/((1-~x)*sqrt(pi))
        # Doesn't support Legendre
        @acrule 1/(sqrt(~x)*(1+~x)) => -sqrt(pi)/(2*(1+~x)^(3/2))
        @acrule 1/(sqrt(~x)*(1-~x)) => sqrt(pi)/(2*(1-~x)^(3/2))
        @acrule sqrt(~x)/~x => sqrt(pi)/(2*(1+~x)^(3/2))
        @acrule sqrt(~x)/(-~x) => sqrt(pi)/(2*(1-~x)^(3/2))
        @acrule sqrt(~x)/(1-~x)^2 => sqrt(pi)*(2+~x)/(4*(1-~x)^(5/2))
        #@acrule((~x)^(~p)/(1-(~x))^((~p)+3/2) => (((~p)+1/2+1/2*(~x))*gamma(~p+1)*(~x)^((~p)-1/2))/(gamma((~p+3/2)*(1-(~x))^(2+(~p))))
        @acrule sqrt((1-~x)/~x) => (ellipe(~x)-ellipk(~x))/(~x*sqrt(pi))
        @acrule sqrt((1+~x)/~x) => ((1+~x)*ellipe(~x/(1+~x)) - ellipk(~x/(1+~x)))/(~x*sqrt(pi*(1+~x)))
        @acrule  (~x)^(~p)/(1+~x)^~r => gamma(~p+1)/gamma(~p+1/2)*(~x)^(~p-1/2)* Hypergeometric1F1(~r, ~p+1, p+1/2, ~x)
        @acrule  (~x)^(~p)/(1-~x)^~r => gamma(~p+1)/gamma(~p+1/2)*(~x)^(~p-1/2)* Hypergeometric1F1(~r, ~p+1, p+1/2, -~x)

        # EXPONENTIAL AND RELATED FUNCTIONS
        @acrule exp(~x) => 1/sqrt(pi*~x) + exp(~x)*erf(sqrt(~x))
        @acrule exp(-~x) => 1/sqrt(pi*~x) - 2/sqrt(pi)*dawson(sqrt(-~x))
        @acrule exp(~x)*erf(sqrt(~x)) => exp(~x)
        @acrule dawson(sqrt(~x)) => 1/2*sqrt(pi)*exp(-~x)
        @acrule exp(~x)*erfc(sqrt(~x)) => 1/sqrt(pi*~x) - exp(~x)*erfc(sqrt(~x))
        @acrule exp(~x)*erfc(-sqrt(~x)) => 1/sqrt(pi*~x) + exp(~x)*erf(-sqrt(~x))
        @acrule erf(~x) => exp(-1/2*~x)*besseli(0, 1/2*~x)
        @acrule exp(~x)/sqrt(~x) => 1/2*sqrt(pi)*exp(1/2*~x)*(besseli(1, 1/2*~x)+besseli(0, 1/2*~x))
        @acrule exp(-~x)/sqrt(~x) => 1/2*sqrt(pi)*exp(-1/2*~x)*(besseli(1, 1/2*~x)-besseli(0, 1/2*~x))
        @acrule cosh(sqrt(~x))/sqrt(~x) => sqrt(pi/~x)*besseli(1, sqrt(~x))/2
        @acrule  (1-cos(sqrt(~x)))/~x => sqrt(pi)/2*Struve0(sqrt(~x))
        @acrule  (1-cosh(sqrt(~x)))/~x => sqrt(pi)/2*MStruve0(sqrt(~x))

        @acrule sin(~x) => sin(~x+pi/4)-sqrt(2)*(AuxiliaryFresnelCos(sqrt(2*~x/π)))
        @acrule cos(~x) => 1/(sqrt(π*~x))+cos(~x+1/4*π)-sqrt(2)*AuxiliaryFresnelSin(sqrt(2*~x/π))
        @acrule sinh(~x) => dawson(sqrt(~x))/sqrt(pi)-(exp(~x)*erf(sqrt(~x)))/2
        @acrule cosh(~x) => 1/sqrt(pi*~x)+exp(~x)*erf(sqrt(~x))/2-dawson(sqrt(~x))/sqrt(pi)
        @acrule asin(sqrt(~x))/sqrt(1-~x) => sqrt(pi)/(2*(1-~x))
        @acrule atan(sqrt(~x)) => 1/2*sqrt(pi/(1+~x))
        @acrule asinh(sqrt(~x))/sqrt(1+~x) => sqrt(pi)/(2*(1+~x))
        @acrule atanh(sqrt(~x)) => 1/2*sqrt(pi/(1-~x))

        # BESSEL AND STRUVE FUNCTIONS
        @acrule besselj(0, sqrt(~x)) => cos(sqrt(~x))/sqrt(pi*(~x))
        @acrule besselj(1, sqrt(~x))/sqrt(~x) => (cos(sqrt(~x))+sqrt(~x)*sin(sqrt(~x))-1)/(sqrt(~x)*(~x)^(3/2))
        @acrule  besselj(~v, sqrt(~x))/(~x)^(~v/2) => 1/(2^(~v)*gamma(~v+1)*sqrt(pi*~x))- Struve(~v+1/2, sqrt(~x))/(sqrt(2)*(~x)^(~v/2+1/4))
        @acrule sqrt(~x)*besselj(1, sqrt(~x)) => sin(sqrt(~x))/sqrt(pi)
        @acrule (~x)^(~v/2)*besselj(1, sqrt(~x)) => (~x)^((~v)/2-1/4)*besselj(~v-1/2, sqrt(~x))/sqrt(2)
        @acrule besseli(0, sqrt(~x)) => cosh(sqrt(~x))/sqrt(pi*~x)
        @acrule besseli(1, sqrt(~x))/sqrt(~x) => (cosh(sqrt(~x))-sqrt(~x)*sinh(sqrt(~x))-1)/(sqrt(pi)*(~x)^(3/2))
        @acrule  besseli(~v, sqrt(~x))/((~x)^(~v/2)) => 1/(2^~v*gamma(~v+1)*sqrt(pi*(~x))) + Legendre(~v+1/2, sqrt(~x))/(sqrt(2)*(~x)^(~v/2+1/4))
        @acrule sqrt(~x)*besseli(1, sqrt(~x)) => sinh(sqrt(~x))/sqrt(pi)
        @acrule  (~x)^(~v/2)*besseli(~v, sqrt(~x)) => (~x)^((~v)/2-1/4)*besseli(~v-1/2, sqrt(~x))/sqrt(2)
        @acrule exp(-~x)*besseli(0, ~x) => exp(-2*~x)/sqrt(pi*(~x))
        @acrule  Struve(0, sqrt(~x)) => sin(sqrt(~x))/sqrt(pi*(~x))
        @acrule  Struve(0, sqrt(~x))/(~x) => (sin(sqrt(~x)) - SinIntegral(sqrt(~x)))/(sqrt(pi)*(~x)^(3/2))
        @acrule sqrt(~x)*Struve(1, sqrt(~x)) =>(1-cos(sqrt(~x)))/sqrt(pi)
        @acrule (~x)^(~v/2)*Struve(~v, sqrt(~x)) => (~x)^((~v)/2-1/4)*Struve(~v-1/2, sqrt(~x))/sqrt(2)

        @acrule MStruve(0, sqrt(~x)) => sinh(sqrt(~x))/sqrt(pi*(~x))
        @acrule MStruve(0, sqrt(~x))/(~x) => (sinh(sqrt(~x)) - HyperSinIntegral(sqrt(~x)))/(sqrt(pi)*(~x)^(3/2))
        @acrule sqrt(~x)*MStruve(1, sqrt(~x)) => (cosh(sqrt(~x)-1))/sqrt(pi)
        @acrule  (~x)^(~v/2)*MStruve(~v, sqrt(~x)) => (~x)^((~v)/2-1/4)*MStruve(~v-1/2, sqrt(~x))/sqrt(2)

        # MISCELLANEOUS FUNCTIONS
        @acrule log(~x) => log(4*~x)/sqrt(pi*~x)
        @acrule sqrt(~x)*log(~x) => sqrt(pi)/2*(log(~x/4)+2)
        @acrule log(~x)/sqrt(~x) => sqrt(pi)/~x

        @acrule ellipk(~x) => sqrt(pi)/(2*sqrt(~x*(1-~x)))
        @acrule ellipe(~x) => 1/2*sqrt(pi*(1-~x)/(~x))

]

DIFFRULES = Chain(SEMIDIFFRULES)

"""
# Symbolic fractional differentiation

    semidiff(fun)

```semidiff``` uses SymbolicUtils.jl and Symbolics.jl to compute symbolic fractional differentiation.

### Example

```julia-repl
julia> using SymbolicUtils
julia> @syms x
julia> semidiff(log(x))
log(4x) / sqrt(πx)
```
"""
semidiff(x) = DIFFRULES(x)


text:
module AWSMetadataUtilities

using ..AWSExceptions
using Base64
using GitHub
using OrderedCollections: OrderedDict, LittleDict
using JSON
using Mocking

"""
    _get_aws_sdk_js_files(repo_name::String, auth::GitHub.OAuth2)

Get a list of all AWS service API definition files from the `awsk-sdk-js` GitHub repository.

# Arguments
- `repo_name::String`: aws-sdk-js repository
- `auth::GitHub.OAuth2`: OAuth2 credentials for authenticated requests

# Returns
- `Array{OrderedDict}`: Array of Dictionaries which contains information about each AWS Service
"""
function _get_aws_sdk_js_files(repo_name::String, auth::GitHub.OAuth2)
    master_tree = @mock tree(repo_name, "master"; auth=auth)
    apis_sha = [t for t in master_tree.tree if t["path"]=="apis"][1]["sha"]
    files = @mock tree(repo_name, apis_sha)
    files = files.tree

    filter!(f -> endswith(f["path"], ".normal.json"), files)

    return _filter_latest_service_version(files)
end

"""
    _filter_latest_service_version(services::AbstractArray)

Return a list of all AWS Services and their latest version.

# Arguments
- `services::AbstractArray`: Array of all AWS APIs and their versions

# Returns
- `Vector{OrderedDict}`: Array of the latest AWS Service definitions
"""
function _filter_latest_service_version(services::AbstractArray)
    seen_services = Set{String}()
    latest_versions = OrderedDict[]

    for service in reverse(services)
        service_name, _ = _get_service_and_version(service["path"])

        if !(service_name in seen_services)
            push!(seen_services, service_name)
            push!(latest_versions, service)
        end
    end

    return latest_versions
end

"""
    _get_service_and_version(filename::String)

Get the `service` and `version` from a filename.

# Arguments
- `filename`: Name of the file, ex. `{Service}-{Version}.normal.json`

# Returns
- `Tuple{String, String}`: `service` and `version` for an AWS service
"""
function _get_service_and_version(filename::String)
    try
        # Remove ".normal.json" suffix
        service_and_version = join(split(filename, '.')[1:end-2],'.')

        service_and_version = split(service_and_version, '-')
        service = join(service_and_version[1:end-3], '-')
        version = join(service_and_version[end-2:end], '-')

        return (service, version)
    catch e
        if e isa BoundsError
            throw(InvalidFileName("$filename is an invalid AWS JSON filename."))
        else
            rethrow()
        end
    end
end

"""
    _generate_low_level_definitions(services::AbstractArray{<:AbstractDict}, repo_name::String, auth::GitHub.OAuth2)

Get the low-level definitions for all AWS Services.

# Arguments
- `services::Abstract{<:AbstractDict}`: List of AWS Services to generate low-level definitions
- `repo_name::String`: aws-sdk-js repository
- `auth::GitHub.OAuth2`: OAuth2 credentials for authenticated requests

# Returns
- `Vector{String}`: Array of low-level service code to be written into `AWSServices.jl`
"""
function _generate_low_level_definitions(services::AbstractArray, repo_name::String, auth::GitHub.OAuth2)
    service_definitions = String[]

    for service in services
        service_name = service["path"]
        @info "Generating low level wrapper for $service_name"

        # Get the contents of the ${Service}.normal.json file
        service_blob = blob(repo_name, service["sha"]; auth=auth)
        service = JSON.parse(String(base64decode(service_blob.content)))
        service_metadata = service["metadata"]

        definition = _generate_low_level_definition(service_metadata)
        push!(service_definitions, definition)
    end

    return service_definitions
end

"""
    _generate_low_level_definition(service::AbstractDict)

Get the low-level definition for an AWS Service.

# Arguments
- `service::AbstractDict`: AWS Service to generate the low-level code definition

# Returns
- `String`: Low-level definition for the AWS Service

# Throws
- `ProtocolNotFound`: If `service` is using a protocol that's not either `[rest-xml, rest-json, json, query]`
"""
function _generate_low_level_definition(service::AbstractDict)
    protocol = service["protocol"]
    service_name = service["endpointPrefix"]
    service_id = replace(lowercase(service["serviceId"]), ' ' => '_')
    api_version = service["apiVersion"]

    service_specifics = LittleDict{String, String}()

    if service_id == "glacier"
        service_specifics[service_id] = "LittleDict(\"x-amz-glacier-version\" => \"$(service["apiVersion"])\")"
    end

    if protocol == "rest-xml"
        return "const $service_id = AWS.RestXMLService(\"$service_name\", \"$api_version\")"
    elseif protocol in ("ec2", "query")
        return "const $service_id = AWS.QueryService(\"$service_name\", \"$api_version\")"
    elseif protocol == "rest-json" && haskey(service_specifics, service_id)
        return "const $service_id = AWS.RestJSONService(\"$service_name\", \"$api_version\", $(service_specifics[service_id]))"
    elseif protocol == "rest-json"
        return "const $service_id = AWS.RestJSONService(\"$service_name\", \"$api_version\")"
    elseif protocol == "json"
        json_version = service["jsonVersion"]
        target = service["targetPrefix"]
        return "const $service_id = AWS.JSONService(\"$service_name\", \"$api_version\", \"$json_version\", \"$target\")"
    else
        throw(ProtocolNotDefined("$service_id is using a new protocol; $protocol which is not supported."))
    end
end

"""
    _clean_documentation(documentation::String)

Clean up the documentation to make it Julia compiler and human-readable.

- Remove any HTML tags
- Remove any dollar signs
- Remove any backslashes
- Escape any double-quotes

# Arguments
- `documentation::String`: Documentation to be cleaned

# Returns
- `String`: Cleaned documentation string
"""
function _clean_documentation(documentation::String)
    documentation = replace(documentation, r"\<.*?\>" => "")
    documentation = replace(documentation, '$' => "")
    documentation = replace(documentation, '\\' => "")
    documentation = replace(documentation, '"' => "\\\"")

    return documentation
end

"""
    _clean_uri(uri::String)

Replace URI parameters with the appropriate syntax for Julia interpolation.

Find all URI parameters, and apply the following replacements:
* { => \$(
* } => )
* - => _
* + => empty string

Example:
"/v1/configurations/{configuration-id}" => "/v1/configurations/\$(configuration_id)"

# Arguments
- `uri::String`: URI to be cleaned

# Returns
- `String`: Cleaned URI
"""
function _clean_uri(uri::String)
    uri_parameters = eachmatch(r"{.*?}", uri) # Match anything surrounded in "{ }"

    for param in uri_parameters
        match = param.match
        original_match = match

        match = replace(match, '{' => "\$(")  # Replace { with $(
        match = replace(match, '}' => ')')  # Replace } with )
        match = replace(match, '-' => '_')  # Replace hyphens with underscores
        match = replace(match, '+' => "")  # Remove +

        uri = replace(uri, original_match => match)
    end

    return uri
end

"""
    _format_function_name(function_name::String)

Convert a function name from CamelCase to snake_case

# Arguments
- `function_name::String`: Function name to be converted to snake_case

# Returns:
- `String`: snake_case version of the function_name
"""
function _format_function_name(function_name::String)
    # Replace a string of uppercase characters with themselves prefaced by an underscore
    # [A-Z](?![A-Z]) => Match a single uppercase character that is not followed by another uppercase character
    # |(A-Z]{1,})    => Match 1-Infinite amounts of uppercase characters
    function_name = replace(function_name, r"[A-Z](?![A-Z])|([A-Z]{1,})" => s"_\g<0>")

    # Lowercase the entire string
    function_name = lowercase(function_name)

    # Chop off the leading underscore
    function_name = chop(function_name, head=1, tail=0)

    return function_name
end

"""
    _get_function_parameters(input, shapes)

Get the required and optional parameters for a given operation.

# Arguments
- `input::String`: The input object for the given operation
- `shapes::AbstractDict{String, Any}`: Dictionary of all shapes for this AWS Service

# Returns
- `Tuple(Dict, Dict)`: (required_parameters, optional_parameters)
"""
function _get_function_parameters(input::String, shapes::AbstractDict{String})
    """
        _get_parameter_name(parameter::String, input_shape::AbstractDict{String, <:Any})

    Find the correct parameter name for making requests. Certain ones have a specific locationName, for Batch requests this locationName is nested one shape deeper.

    # Arguments
    - `parameter::String`: Name of the original parameter
    - `input_shape::AbstractDict{String, <:Any}`: The parameter shape

    # Returns
    - `String`: Either the original parameter name, the locationName for the parameter, or the locationName nested one shape deeper
    """
    function _get_parameter_name(parameter::String, input_shape::AbstractDict{String, <:Any})
        # If the parameter has a locationName, return it
        if haskey(input_shape["members"][parameter], "locationName")
            return input_shape["members"][parameter]["locationName"]
        end

        # Check to see if the nested shape has a locationName
        nested_shape = input_shape["members"][parameter]["shape"]
        nested_member = get(shapes[nested_shape], "member", Dict())

        # If nested_shape[member] exists return locationName (if exists), otherwise return the original parameter name
        return get(nested_member, "locationName", parameter)
    end

    required_parameters = LittleDict{String, Any}()
    optional_parameters = LittleDict{String, Any}()
    input_shape = shapes[input]

    if haskey(input_shape, "required")
        for parameter in input_shape["required"]
            parameter_name = _get_parameter_name(parameter, input_shape)

            # Check if the parameter needs to be in a certain place
            parameter_location = get(input_shape["members"][parameter], "location", "")

            documentation = _clean_documentation(get(input_shape["members"][parameter], "documentation", ""))

            required_parameters[parameter_name] = LittleDict{String, String}(
                "location" => parameter_location,
                "documentation" => documentation,
            )
        end
    end

    if haskey(input_shape, "members")
        for (member_key, member_value) in input_shape["members"]
            parameter_name = get(input_shape["members"][member_key], "locationName", member_key)

            if !haskey(required_parameters, parameter_name)
                documentation = _clean_documentation(get(member_value, "documentation", ""))
                idempotent = get(member_value, "idempotencyToken", false)

                optional_parameters[parameter_name] = LittleDict{String, Union{String, Bool}}(
                    "documentation" => documentation,
                    "idempotent" => idempotent,
                )
            end
        end
    end

    return (sort(required_parameters), sort(optional_parameters))
end

"""
    _generate_high_level_definitions(
        service_name::String,
        protocol::String,
        operations::AbstractDict{String, Any},
        shapes::AbstractDict{String, Any}
    )

Generate high-level definitions for the `service`.

# Arguments
- `service_name::String`: AWS `service`
- `protocol::String`: API protocol type
- `operations::AbstractDict{String, Any}`: Service API operations
- `shapes::AbstractDict{String, Any}`: All input shapes for this service

# Return
- `Array`: Array of all high-level definitions and documentation to be written into `services/{Service}.jl`
"""
function _generate_high_level_definitions(
    service_name::String,
    protocol::String,
    operations::AbstractDict,
    shapes::AbstractDict
)
    operation_definitions = String[]

    for (_, operation) in operations
        name = operation["name"]
        method = operation["http"]["method"]
        request_uri = operation["http"]["requestUri"]

        documentation = ""

        if haskey(operation, "documentation")
            documentation = _clean_documentation(operation["documentation"])
        end

        required_parameters = Dict{String, Any}()
        optional_parameters = Dict{String, Any}()

        if haskey(operation, "input")
            required_parameters, optional_parameters = _get_function_parameters(operation["input"]["shape"], shapes)
        end

        operation_definition = _generate_high_level_definition(
            service_name,
            protocol,
            name,
            method,
            request_uri,
            required_parameters,
            optional_parameters,
            documentation
        )

        push!(operation_definitions, operation_definition)
    end

    return operation_definitions
end

"""
    _generate_high_level_definition(
        service_name::String,
        protocol::String,
        name::String,
        method::String,
        request_uri::String,
        required_parameters::Dict,
        optional_parameters::Dict,
        documentation::String
    )

Generate the high-level definition for a services function.

# Arguments
- `name::String`: Name of this services function
- `service_name::String`: Name of the AWS Service
- `protocol::String`: Request protocol for the AWS Service
- `method::String`: Request method for this services function
- `request_uri::String`: Endpoint for this services function
- `required_parameters::AbstractDict{Any, Any}`: Required parameters for this services function
- `optional_parameters::AbstractDict{Any, Any}`: Optional parameters for this services function
- `documentation::String`: Documentation for this services function

# Return
- `String`: High-level definition in Julia code for the services function
"""
function _generate_high_level_definition(
    service_name::String,
    protocol::String,
    name::String,
    method::String,
    request_uri::String,
    required_parameters::AbstractDict,
    optional_parameters::AbstractDict,
    documentation::String
)
    """
        _generate_rest_operation_defintion(required_params::AbstractDict{String, Any}, optional_params::AbstractDict{String, Any}, function_name::String, service_name::String, method::String, request_uri::String)

    Generate function definition for a service request given required, header and idempotent parameters.

    # Arguments
    - `required_params::AbstractDict{String, Any}`: All required parameters which the function needs
    - `optional_params::AbstractDict{String, Any}`: All optional parameters which can be passed in
    - `function_name::String`: Name of the function which we are generating a definition for
    - `service_name::String`: Name of the AWS service which this function is being defined for
    - `method::String`: The REST method that this function requires
    - `request_uri::String`: The endpoint of the function

    # Returns
    - `String`: Function definition as Julia code
    """
    function _generate_rest_operation_defintion(required_params::AbstractDict, optional_params::AbstractDict, function_name::String, service_name::String, method::String, request_uri::String)
        request_uri = replace(request_uri, '{' => "\$(")  # Replace { with $(
        request_uri = replace(request_uri, '}' => ')')  # Replace } with )
        request_uri = replace(request_uri, '+' => "")  # Remove + from the request URI

        # Pre Julia-1.3 workaround
        req_keys = [replace(key, "-"=>"_") for key in collect(keys(required_params))]

        required_params = filter(p -> (p[2]["location"] != "uri"), required_params)
        header_params = filter(p -> (p[2]["location"] == "header"), required_params)
        required_params = setdiff(required_params, header_params)

        idempotent_params = filter(p -> (p[2]["idempotent"]), optional_params)

        req_kv = ["\"$(p[1])\"=>$(replace(p[1], "-" => "_"))" for p in required_params]
        header_kv = ["\"$(p[1])\"=>$(replace(p[1], "-" => "_"))" for p in header_params]
        idempotent_kv = ["\"$(p[1])\"=>string(uuid4())" for p in idempotent_params]

        required_keys = !isempty(req_keys)
        headers = !isempty(header_params)
        idempotent = !isempty(idempotent_params)

        req_str = !isempty(req_kv) ? "Dict{String, Any}($(join(req_kv, ", "))" : ""
        params_str = (!isempty(req_kv) || idempotent) ? "$(join(vcat(req_kv, idempotent_kv), ", "))" : ""
        headers_str = headers ? "\"headers\"=>Dict{String, Any}($(join(header_kv, ", ")))" : ""
        params_headers_str = "Dict{String, Any}($(join([s for s in (params_str, headers_str) if !isempty(s)], ", ")))"

        formatted_function_name = _format_function_name(function_name)

        if required_keys && (idempotent || headers)
            return """\n
                $formatted_function_name($(join(req_keys, ", ")); aws_config::AWSConfig=global_aws_config()) = $service_name(\"$method\", \"$request_uri\", $params_headers_str; aws_config=aws_config)
                $formatted_function_name($(join(req_keys, ", ")), args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = $service_name(\"$method\", \"$request_uri\", Dict{String, Any}(mergewith(_merge, $params_headers_str, args)); aws_config=aws_config)
                """
        elseif !required_keys && (idempotent || headers)
            return """\n
                $formatted_function_name(; aws_config::AWSConfig=global_aws_config()) = $service_name(\"$method\", \"$request_uri\", $params_headers_str; aws_config=aws_config)
                $formatted_function_name(args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = $service_name(\"$method\", \"$request_uri\", Dict{String, Any}(mergewith(_merge, $params_headers_str, args)); aws_config=aws_config)
                """
        elseif required_keys && !isempty(req_kv)
            return """\n
                $formatted_function_name($(join(req_keys, ", ")); aws_config::AWSConfig=global_aws_config()) = $service_name(\"$method\", \"$request_uri\", $req_str); aws_config=aws_config)
                $formatted_function_name($(join(req_keys, ", ")), args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = $service_name(\"$method\", \"$request_uri\", Dict{String, Any}(mergewith(_merge, $req_str), args)); aws_config=aws_config)
                """
        elseif required_keys
            return """\n
                $formatted_function_name($(join(req_keys, ", ")); aws_config::AWSConfig=global_aws_config()) = $service_name(\"$method\", \"$request_uri\"; aws_config=aws_config)
                $formatted_function_name($(join(req_keys, ", ")), args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = $service_name(\"$method\", \"$request_uri\", args; aws_config=aws_config)
                """
        else
            return """\n
                $formatted_function_name(; aws_config::AWSConfig=global_aws_config()) = $service_name(\"$method\", \"$request_uri\"; aws_config=aws_config)
                $formatted_function_name(args::AbstractDict{String, Any}; aws_config::AWSConfig=global_aws_config()) = $service_name(\"$method\", \"$request_uri\", args; aws_config=aws_config)
                """
        end
    end

    """
        _generate_json_query_operation_definition(required_params::AbstractDict, optional_params::AbstractDict, function_name::String, service_name::String)

    Generate a JSON/Query high level definition.

    # Arguments
    - `required_params::AbstractDict`: Required parameters for the function
    - `optional_params::AbstractDict`: Optional parmaeters for the function
    - `function_name::String`: Name of the function to generate code
    - `service_name`: Name of the service associated with the function

    # Returns
    - `String`: Function definition as Julia code
    """
    function _generate_json_query_opeation_definition(required_params::AbstractDict, optional_params::AbstractDict, function_name::String, service_name::String)
        req_keys = [replace(key, '-' => '_') for key in collect(keys(required_params))]

        idempotent_params = filter(p -> (p[2]["idempotent"]), optional_params)

        req_kv = ["\"$(p[1])\"=>$(replace(p[1], "-" => "_"))" for p in required_params]
        idempotent_kv = ["\"$(p[1])\"=>string(uuid4())" for p in idempotent_params]

        required = !isempty(req_kv)
        idempotent = !isempty(idempotent_kv)

        formatted_function_name = _format_function_name(function_name)

        if required && idempotent
            return """\n
                $formatted_function_name($(join(req_keys, ", ")); aws_config::AWSConfig=global_aws_config()) = $service_name(\"$function_name\", Dict{String, Any}($(join(req_kv, ", ")), $(join(idempotent_kv, ", "))); aws_config=aws_config)
                $formatted_function_name($(join(req_keys, ", ")), args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = $service_name(\"$function_name\", Dict{String, Any}(mergewith(_merge, Dict{String, Any}($(join(req_kv, ", ")), $(join(idempotent_kv, ", "))), args)); aws_config=aws_config)
                """
        elseif required
            return """\n
                $formatted_function_name($(join(req_keys, ", ")); aws_config::AWSConfig=global_aws_config()) = $service_name(\"$function_name\", Dict{String, Any}($(join(req_kv, ", "))); aws_config=aws_config)
                $formatted_function_name($(join(req_keys, ", ")), args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = $service_name(\"$function_name\", Dict{String, Any}(mergewith(_merge, Dict{String, Any}($(join(req_kv, ", "))), args)); aws_config=aws_config)
                """
        elseif idempotent
            return """\n
                $formatted_function_name(; aws_config::AWSConfig=global_aws_config()) = $service_name(\"$function_name\", Dict{String, Any}($(join(idempotent_kv, ", "))); aws_config=aws_config)
                $formatted_function_name(args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = $service_name(\"$function_name\", Dict{String, Any}(mergewith(_merge, Dict{String, Any}($(join(idempotent_kv, ", "))), args)); aws_config=aws_config)
                """
        else
            return """\n
                $formatted_function_name(; aws_config::AWSConfig=global_aws_config()) = $service_name(\"$function_name\"; aws_config=aws_config)
                $formatted_function_name(args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = $service_name(\"$function_name\", args; aws_config=aws_config)
                """
        end
    end

    """
        _generate_docstring(function_name, documentation, required_parameters, optional_parameters)

    Generate the docstring for the `function_name`.

    # Arguments
    - `function_name`: Name of the function
    - `documentation`: Documentation associated with the function
    - `required_parameters`: Required parameters for the function
    - `optional_parameters`: Optional parameters for the function

    # Returns
    - `String`: Docstring for the function
    """
    function _generate_docstring(function_name, documentation, required_parameters, optional_parameters)
        operation_definition = """
            $(repeat('"', 3))
                $function_name()

            $documentation\n
            """

        # Add in the required parameters if applicable
        if !isempty(required_parameters)
            operation_definition *= "# Required Parameters\n"

            for (required_key, required_value) in required_parameters
                operation_definition *= "- `$required_key`: $(required_value["documentation"])\n"
            end

            operation_definition *= "\n"
        end

        # Add in the optional parameters if applicable
        if !isempty(optional_parameters)
            operation_definition *= "# Optional Parameters\n"

            for (optional_key, optional_value) in optional_parameters
                operation_definition *= "- `$optional_key`: $(optional_value["documentation"])\n"
            end
        end

        operation_definition *= repeat('"', 3)
    end

    operation_definition = _generate_docstring(name, documentation, required_parameters, optional_parameters)

    if protocol in ("json", "query", "ec2")
        operation_definition *= _generate_json_query_opeation_definition(required_parameters, optional_parameters, name, service_name)
    elseif protocol in ("rest-json", "rest-xml")
        operation_definition *= _generate_rest_operation_defintion(required_parameters, optional_parameters, name, service_name, method, request_uri)
    else
        throw(ProtocolNotDefined("$function_name is using a new protocol; $protocol which is not supported."))
    end

    return operation_definition
end
end  # Module


text:
let
    frame = Frame([1,2,3])
    @test length(frame) == 3
    @test capacity(frame) == 3
    for i in 1 : 3
        @test frame[i] == i
    end

    frame = Frame([1,2,3], capacity=5)
    @test length(frame) == 3
    @test capacity(frame) == 5
    for i in 1 : 3
        @test frame[i] == i
    end

    @test_throws ErrorException Frame([1,2,3], capacity=2)

    frame = Frame(Int)
    @test length(frame) == 0
    @test capacity(frame) > 0

    frame = Frame(Int, 2)
    @test length(frame) == 0
    @test capacity(frame) == 2

    frame[1] = 999
    frame[2] = 888

    @test frame[1] == 999
    @test frame[2] == 888
    @test length(frame) == 0 # NOTE: length does not change
    @test capacity(frame) == 2

    empty!(frame)
    @test length(frame) == 0
    @test capacity(frame) == 2

    push!(frame, 999)
    push!(frame, 888)
    @test length(frame) == 2
    @test capacity(frame) == 2

    @test_throws BoundsError push!(frame, 777)

    frame = Frame([999,888])
    deleteat!(frame, 1)
    @test length(frame) == 1
    @test capacity(frame) == 2
    @test frame[1] == 888

    deleteat!(frame, 1)
    @test length(frame) == 0
    @test capacity(frame) == 2

    frame = Frame([1,2,3])
    frame2 = copy(frame)
    for i in 1 : 3
        @test frame[i] == frame2[i]
    end
    frame[1] = 999
    @test frame2[1] == 1
end

let
    frame = EntityFrame(Int, Float64, String)
    @test eltype(frame) == Entity{Int,Float64,String}

    frame = Frame([Entity(1,1,"A"),Entity(2,2,"B"),Entity(3,3,"C")])
    @test  in("A", frame)
    @test  in("B", frame)
    @test  in("C", frame)
    @test !in("D", frame)
    @test findfirst("A", frame) == 1
    @test findfirst("B", frame) == 2
    @test findfirst("C", frame) == 3
    @test findfirst("D", frame) == nothing
    @test id2index(frame, "A") == 1
    @test id2index(frame, "B") == 2
    @test id2index(frame, "C") == 3
    @test_throws BoundsError id2index(frame, "D")

    frame = Frame([Entity(1,1,"A"),Entity(2,2,"B"),Entity(3,3,"C")])
    @test get_by_id(frame, "A") == frame[1]
    @test get_by_id(frame, "B") == frame[2]
    @test get_by_id(frame, "C") == frame[3]

    delete!(frame, Entity(2,2,"B"))
    @test frame[1] == Entity(1,1,"A")
    @test frame[2] == Entity(3,3,"C")
    @test length(frame) == 2

    delete!(frame, "A")
    @test frame[1] == Entity(3,3,"C")
    @test length(frame) == 1

    frame = Frame([Entity(1,1,1),Entity(2,2,2)], capacity=3)
    @test get_first_available_id(frame) == 3
end

max_stars_repo_name: tawheeler/Records

max_issues_repo_name: tawheeler/Records

max_forks_repo_name: tawheeler/Records


text:
function get_data_paths(config)
    data_paths = []
    for i = 1:nrow(config)
        row = config[i, :]
        file = string(
            "n_factors_",
            row.n_factors,
            "_n_items_",
            row.n_items,
            "_meanstructure_",
            row.meanstructure,
            ".csv"
            )
        push!(data_paths, file)
    end
    return data_paths
end

function read_files(dir, data_paths)
    data = Vector{Any}()
    for i in 1:length(data_paths)
        push!(data, DataFrame(CSV.File(dir*"/"*data_paths[i])))
    end
    return data
end

function gen_CFA_RAM(nfact, nitem)
    nfact = Int64(nfact)
    nitem = Int64(nitem)

    ## Model definition
    nobs = nfact*nitem
    nnod = nfact+nobs
    n_latcov = Int64(nfact*(nfact-1)/2)
    npar = 2nobs + n_latcov
    Symbolics.@variables x[1:npar]

    #F
    Ind = collect(1:nobs)
    Jnd = collect(1:nobs)
    V = fill(1,nobs)
    F = sparse(Ind, Jnd, V, nobs, nnod)

    #A
    Ind = collect(1:nobs)
    Jnd = vcat([fill(nobs+i, nitem) for i in 1:nfact]...)
    V = [x...][1:nobs]
    A = sparse(Ind, Jnd, V, nnod, nnod)

    #S
    Ind = collect(1:nnod)
    Jnd = collect(1:nnod)
    V = [[x...][nobs+1:2nobs]; fill(1.0, nfact)]
    S = sparse(Ind, Jnd, V, nnod, nnod)
    xind = 2nobs+1
    for i in nobs+1:nnod
        for j in i+1:nnod
            S[i,j] = x[xind]
            S[j,i] = x[xind]
            xind = xind+1
        end
    end

    return RAMMatrices(;A = A, S = S, F = F, parameters = x)
end

function gen_model(nfact, nitem, data, estimator, backend)

    ram_matrices = gen_CFA_RAM(nfact, nitem)

    if backend == "Optim.jl"
        semdiff = SemDiffOptim(
            LBFGS(
                ;linesearch = BackTracking(order=3),
                alphaguess = InitialHagerZhang()
                ),
            Optim.Options(
                ;f_tol = 1e-10,
                x_tol = 1.5e-8)
            )
    elseif backend =="NLopt.jl"
        semdiff = SemDiffNLopt(
            :LD_LBFGS,
            nothing
            )
    end

    # imply and loss
    if estimator == "ML"
        model = Sem(
            ram_matrices = ram_matrices,
            data = data,
            imply = RAM,
            diff = semdiff
        )
    elseif estimator == "GLS"
        model = Sem(
            ram_matrices = ram_matrices,
            data = data,
            imply = RAM,
            loss = (SemWLS, ),
            diff = semdiff
        )
    else
        error("unknown estimator")
    end

    return model
end

function gen_models(config, data_vec)
    models = []
    for i = 1:nrow(config)
        row = config[i, :]
        model = gen_model(row.n_factors, row.n_items, Matrix(data_vec[i]), row.Estimator, row.backend)
        push!(models, model)
    end
    return models
end

function benchmark_models(models)
    benchmarks = []
    for model in models
        bm = @benchmark sem_fit($model)
        push!(benchmarks, bm)
    end
    return benchmarks
end

function get_fits(models)
    fits = []
    for model in models
        fit = sem_fit(model)
        push!(fits, fit)
    end
    return fits
end

function compare_estimates(fits, par_vec, config)
    correct =
        [compare_estimate(fit, estimate, n_factors, n_items) for
                (fit, estimate, n_factors, n_items) in zip(fits, par_vec, config.n_factors, config.n_items)]
    return correct
end

function compare_estimate(fit, estimate, n_factors, n_items)

    nfact = Int64(n_factors)
    nitem = Int64(n_items)

    nobs = nfact*nitem
    nnod = nfact+nobs

    nrows = size(estimate, 1)

    solution = fit.minimizer
    cov_ind = Int(nrows-nfact*(nfact-1)/2+1):nrows
    par_ind = [1:2nobs..., cov_ind...]

    return StructuralEquationModels.compare_estimates(estimate.est[par_ind], solution, 0.01)

end


text:
module SSLTestLandCover
using SimpleSDMLayers
using Test

lc1 = SimpleSDMPredictor(EarthEnv, LandCover, 1; left=-5.0, right=7.0, bottom=30.0, top=45.0)
@test typeof(lc1) <: SimpleSDMPredictor

lc1 = SimpleSDMPredictor(EarthEnv, LandCover, 1; full=true, left=-5.0, right=7.0, bottom=30.0, top=45.0)
@test typeof(lc1) <: SimpleSDMPredictor

lc1 = SimpleSDMPredictor(EarthEnv, LandCover, 1:3; left=-5.0, right=7.0, bottom=30.0, top=45.0)
@test eltype(lc1) <: SimpleSDMPredictor

he1 = SimpleSDMPredictor(EarthEnv, HabitatHeterogeneity, 1; left=-5.0, right=7.0, bottom=30.0, top=45.0)
@test typeof(he1) <: SimpleSDMPredictor

he1 = SimpleSDMPredictor(EarthEnv, HabitatHeterogeneity, 1:3; left=-5.0, right=7.0, bottom=30.0, top=45.0)
@test eltype(he1) <: SimpleSDMPredictor

end


text:
using Printf

function rd(num::Float64, digits::Int64=2) # round
    if digits == 2
      return pd(@printf("%.2e", num))
    elseif digits == 3
      return pd(@printf("%.3e", num))
    elseif digits == 4
      return pd(@printf("%.4e", num))
    elseif digits == 5
      return pd(@printf("%.5e", num))
    else
      error("digits $digits not supported for rd")
    end
end

function pd(input)
    return pd(string(input))
end

function pd(str::String)
    return rpad(str, 11 ," ")
end

function pd(input, pad_len::Int64)
    return pd(string(input), pad_len)
end

function pd(str::String, pad_len::Int64)
    return rpad(str, pad_len ," ")
end

function solve_quadratic(a::Float64, b::Float64, c::Float64)
    val = -b/(2 * a)
    pm = sqrt(b^2 - 4 * a * c) / (2 * a)

    return (val + pm, val - pm)
end

function my_warn(str::String)
    println("WARNING: ", str)
end

function densest_col(mat) #::SparseMatrixCSC{Float64,Int32})
    density = 0
    for i = 1:size(mat,2)
      density = max(density, nnz(mat[:,i]))
    end

    return density
end

function densest_row(mat) #::SparseMatrixCSC{Float64,Int32})
    #=density = 0
    for i = 1:size(mat,1)
      density = max(density, nnz(mat[i,:]))
    end

    return density=#
    return densest_col(mat')
end


text:
using DelayDiffEq, DiffEqBase
using Test

@test DiffEqBase.undefined_exports(DelayDiffEq) == []

max_stars_repo_name: JuliaDiffEq/DelayDiffEq

max_issues_repo_name: JuliaDiffEq/DelayDiffEq

max_forks_repo_name: JuliaDiffEq/DelayDiffEq


text:
## Analyzing memory allocation

immutable MallocInfo
    bytes::Int
    filename::Compat.UTF8String
    linenumber::Int
end

sortbybytes(a::MallocInfo, b::MallocInfo) = a.bytes < b.bytes

"""
    analyze_malloc_files(files) -> Vector{MallocInfo}

Iterates through the given list of filenames and return a `Vector` of
`MallocInfo`s with allocation information.
"""
function analyze_malloc_files(files)
    bc = MallocInfo[]
    for filename in files
        open(filename) do file
            for (i,ln) in enumerate(eachline(file))
                tln = strip(ln)
                if !isempty(tln) && isdigit(tln[1])
                    s = split(tln)
                    b = parse(Int, s[1])
                    push!(bc, MallocInfo(b, filename, i))
                end
            end
        end
    end
    sort(bc, lt=sortbybytes)
end

function find_malloc_files(dirs)
    files = Compat.String[]
    for dir in dirs
        filelist = readdir(dir)
        for file in filelist
            file = joinpath(dir, file)
            if isdir(file)
                append!(files, find_malloc_files(file))
            elseif endswith(file, "jl.mem")
                push!(files, file)
            end
        end
    end
    files
end
find_malloc_files(file::Compat.String) = find_malloc_files([file])

analyze_malloc(dirs) = analyze_malloc_files(find_malloc_files(dirs))
analyze_malloc(dir::Compat.String) = analyze_malloc([dir])

isfuncexpr(ex::Expr) =
    ex.head == :function || (ex.head == :(=) && typeof(ex.args[1]) == Expr && ex.args[1].head == :call)
isfuncexpr(arg) = false

# Support Unix command line usage like `julia Coverage.jl $(find ~/.julia/v0.3 -name "*.jl.mem")`
if !isinteractive()
    bc = analyze_malloc_files(ARGS)
    println(bc)
end


text:
mutable struct Label{T <: Union{String,Vector{String}}}
    text::T
    scale::Float64
    position::Luxor.Point
    halign::Symbol
    valign::Symbol
    angle::Float64
    leading::Float64
end

function Label( ;
    text,
    scale::Real,
    position,
    halign,
    valign,
    angle::Real,
    leading::Real,
)
    return Label(
        text,
        convert(Float64, scale),
        position,
        halign,
        valign,
        convert(Float64, angle),
        convert(Float64, leading),
    )
end

# labstr = Label( ; text="hieeeee", scale=1.1, position=Luxor.Point(0,10), halign=:center, valign=:middle, angle=0.0, leading=0)
# labvec = Label( ; text=["hieeeee","leoooooo"], scale=1.1, position=Luxor.Point(0,10), halign=:center, valign=:middle, angle=0.0, leading=0)
# labmis = Label( ; text=missing, scale=1.1, position=Luxor.Point(0,10), halign=:center, valign=:middle, angle=0.0, leading=0)


text:
################################################################################

# three-layered autoencoder
"""

    buildAutoencoder(inputLayer::T;
    nnParams) where T <: Integer

# Description
Build a three-layered autoencoder.

# Arguments
`inputLayer` number of neurons on input.

`nnParams` neural network hyperparameters.

See also: [`modelTrain!`](@ref)
"""
function buildAutoencoder(inputLayer::T; nnParams) where T <: Integer
  @info("Building three-layered autoencoder...")
  args = nnParams()
  return Chain(
    Dense(inputLayer, args.λ, args.σ),
    Dense(args.λ, inputLayer, args.σ),
  )
end

################################################################################

max_stars_repo_name: DanielRivasMD/MindReader

max_issues_repo_name: DanielRivasMD/MindReader

max_forks_repo_name: DanielRivasMD/MindReader


text:
function combining_sigma(sigma::Dict; rules="Lorentz-Berthelot")
    components = [collect(i)[1] for i in keys(sigma)]
    pairs = [Set(i) for i in collect(combinations(components, 2))]
    combined_sigmas = Dict{Set{String},Float64}()
    for pair in pairs
        if !haskey(sigma,union(pair))
            if rules == "Lorentz-Berthelot"
                combined_sigma = (sum(sigma[Set([i])] for i in pair))/2
            end
            push!(combined_sigmas, pair => combined_sigma)
        end
    end
    return combined_sigmas
end

function combining_epsilon(epsilon::Dict, sigma::Dict, k::Dict; rules_k="Lorentz-Berthelot-mod", rules_no_k="Lorentz-Berthelot")
    components = [collect(i)[1] for i in keys(epsilon)]
    pairs = [Set(i) for i in collect(combinations(components, 2))]
    combined_epsilons = Dict{Set{String},Float64}()
    for pair in pairs
        if !haskey(epsilon,union(pair))
            if haskey(k, pair)
                if rules_k == "Lorentz-Berthelot-mod"
                    combined_epsilon = (1-k[pair])*sqrt(prod(epsilon[Set([i])] for i in pair))
                end
            else
                if rules_no_k == "Lorentz-Berthelot"
                    combined_epsilon = sqrt(prod(epsilon[Set([i])] for i in pair))
                elseif rules_no_k == "Hudson-McCoubrey"
                    combined_epsilon = 2^3*sqrt(prod(sigma[Set([i])]^3 for i in pair))/sum(sigma[Set([i])] for i in pair)^3*sqrt(prod(epsilon[Set([i])] for i in pair))
                end
            end
            push!(combined_epsilons, pair => combined_epsilon)
        end
    end
    return combined_epsilons
end

function combining_lambda(lambda::Dict; rules="Lorentz-Berthelot")
    components = [collect(i)[1] for i in keys(lambda)]
    pairs = [Set(i) for i in collect(combinations(components, 2))]
    combined_lambdas = Dict{Set{String},Float64}()
    for pair in pairs
        if !haskey(lambda,union(pair))
            if rules == "Lorentz-Berthelot"
                combined_lambda = 3+sqrt(prod(lambda[Set([i])]-3 for i in pair))
            end
            push!(combined_lambdas, pair => combined_lambda)
        end
    end
    return combined_lambdas
end