text:
"""
```
shock_decompositions(m, system, histshocks)

shock_decompositions(system, forecast_horizons, histshocks, start_index,
    end_index)

shock_decompositions(m, system, histshocks, start_date, end_date)

shock_decompositions(system, forecast_horizons, histshocks, start_index,
    end_index, regime_inds)
```

### Inputs

- `system::System{S}` or `RegimeSwitchingSystem{S}`: state-space system matrices
- `forecast_horizons::Int`: number of periods ahead to forecast
- `histshocks::Matrix{S}`: matrix of size `nshocks` x `hist_periods` of
  historical smoothed shocks
- `start_index::Int`: first index from which to return computed shock
  decompositions
- `end_index::Int`: last index for which to return computed shock decompositions
- `start_date::Date`: initial date for deterministic trends
- `end_date::Date`: final date for deterministic trends
- `regime_inds::Vector{UnitRange}`: indices of the data corresponding to each regime.

where `S<:AbstractFloat`.

### Outputs

- `states::Array{S, 3}`: matrix of size `nstates` x `nperiods` x `nshocks` of
  state shock decompositions
- `obs::Array{S, 3}`: matrix of size `nobs` x `nperiods` x `nshocks` of
  observable shock decompositions
- `pseudo::Array{S, 3}`: matrix of size `npseudo` x `nperiods` x `nshocks` of
  pseudo-observable shock decompositions

where `nperiods = `end_index - start_index + 1`.
"""
function shock_decompositions(m::AbstractDSGEModel,
    system::System{S}, histshocks::Matrix{S}) where {S<:AbstractFloat}

    horizon   = forecast_horizons(m)
    start_ind = index_shockdec_start(m)
    end_ind   = index_shockdec_end(m)

    shock_decompositions(system, horizon, histshocks, start_ind, end_ind)
end

function shock_decompositions(system::System{S},
    forecast_horizons::Int, histshocks::Matrix{S},
    start_index::Int, end_index::Int) where {S<:AbstractFloat}

    # Setup
    nshocks      = size(system[:RRR], 2)
    nstates      = size(system[:TTT], 2)
    nobs         = size(system[:ZZ], 1)
    npseudo      = size(system[:ZZ_pseudo], 1)
    histperiods  = size(histshocks, 2)
    allperiods   = histperiods + forecast_horizons

    states = zeros(S, nstates, allperiods, nshocks)
    obs    = zeros(S, nobs, allperiods, nshocks)
    pseudo = zeros(S, npseudo, allperiods, nshocks)

    # Check dates
    if forecast_horizons <= 0 || start_index < 1 || end_index > allperiods
        throw(DomainError("At least one of forecast_horizons, the index for the start date of the shock decomposition, and the index for the end date of the shock decomposition is outside the permissible range."))
    end

    # Set constant system matrices to 0
    system = zero_system_constants(system)

    z0 = zeros(S, nstates)
    for i = 1:nshocks
        # Isolate single shock
        shocks = zeros(S, nshocks, allperiods)
        shocks[i, 1:histperiods] = histshocks[i, :]

        # Iterate state space forward
        states[:, :, i], obs[:, :, i], pseudo[:, :, i], _ = forecast(system, z0, shocks)
    end

    # Return shock decompositions in appropriate range
    if start_index == 1 && end_index == allperiods
        return states, obs, pseudo
    else
        range = start_index:end_index
        return states[:, range, :], obs[:, range, :], pseudo[:, range, :]
    end
end

function shock_decompositions(m::AbstractDSGEModel{S},
                              system::RegimeSwitchingSystem{S}, histshocks::Matrix{S},
                              start_date::Dates.Date = date_presample_start(m),
                              end_date::Dates.Date = prev_quarter(date_forecast_start(m))) where {S<:AbstractFloat}
    # Set up dates
    horizon     = forecast_horizons(m)
    start_index = index_shockdec_start(m)
    end_index   = index_shockdec_end(m)

    # Set up regime indices. Note that we do not need to account
    # for ZLB split b/c `histshocks` should have zeros for anticipated shocks
    # in the pre-ZLB periods.
    regime_inds = regime_indices(m, start_date, end_date)
    if regime_inds[1][1] < 1 # remove periods occuring before desired start date
        regime_inds[1] = 1:regime_inds[1][end]
    end

    shock_decompositions(system, horizon, histshocks, start_index, end_index, regime_inds)
end

function shock_decompositions(system::RegimeSwitchingSystem{S},
                              forecast_horizons::Int, histshocks::Matrix{S},
                              start_index::Int, end_index::Int,
                              regime_inds::Vector{UnitRange{Int}}) where {S<:AbstractFloat}

    # Setup
    nshocks     = size(system[1, :RRR], 2)
    nstates     = size(system[1, :TTT], 2)
    nobs        = size(system[1, :ZZ], 1)
    npseudo     = size(system[1, :ZZ_pseudo], 1)
    histperiods = size(histshocks, 2)
    allperiods  = histperiods + forecast_horizons

    states = zeros(S, nstates, allperiods, nshocks)
    obs    = zeros(S, nobs,    allperiods, nshocks)
    pseudo = zeros(S, npseudo, allperiods, nshocks)
    shocks = zeros(S, nshocks, histperiods)

    # Check dates
    if forecast_horizons <= 0 || start_index < 1 || end_index > allperiods
        throw(DomainError("At least one of forecast_horizons, the index for the start date of the shock decomposition, and the index for the end date of the shock decomposition is outside the permissible range."))
    end

    # Set constant system matrices to 0 and start at stationary steady state
    system = zero_system_constants(system)
    z0     = zeros(S, nstates)

    fcast_shocks = zeros(S, nshocks, forecast_horizons) # these are always zero
    for i = 1:nshocks
        # Isolate single shock
        shocks[i, :] = histshocks[i, :]
        # Use forecast to iterate state-space system forward without shocks or the ZLB procedure
        for (reg_num, reg_ind) in enumerate(regime_inds)
            init_state = (reg_num == 1) ? zeros(S, nstates) : states[:, reg_ind[1] - 1, i] # update initial state

            states[:, reg_ind, i], obs[:, reg_ind, i], pseudo[:, reg_ind, i], _ = forecast(system[reg_num], init_state, shocks[:, reg_ind])
        end

        fcast_inds = (histperiods + 1):allperiods
        states[:, fcast_inds, i], obs[:, fcast_inds, i], pseudo[:, fcast_inds, i], _ =
            forecast(system[length(regime_inds)], states[:, regime_inds[end][end], i], fcast_shocks)

        # zero out shocks
        shocks[i, :] .= 0.
    end

    # Return shock decompositions in appropriate range
    if start_index == 1 && end_index == allperiods
        return states, obs, pseudo
    else
        range = start_index:end_index
        return states[:, range, :], obs[:, range, :], pseudo[:, range, :]
    end
end

"""
```
deterministic_trends(m, system, z0)

deterministic_trends(system, z0, nperiods, start_index, end_index)

deterministic_trends(m, system, z0, start_date, end_date)

deterministic_trends(system, z0, nperiods, regime_inds, regimes)
```

Compute deterministic trend values of states, observables, and
pseudo-observables, given a model object and system matrices. The deterministic
trend for a single draw is simply the series that would be obtained by iterating
the state-space system forward, beginning from a state vector `z0` in the last
presample period.

### Inputs

- `m::AbstractDSGEModel`: model object
- `system::System{S}` or `RegimeSwitchingSystem`: state-space system matrices
- `z0`::Vector{S}: initial state vector
- `start_date::Date`: initial date for deterministic trends
- `end_date::Date`: final date for deterministic trends
- `regime_inds::Vector{UnitRange}`: indices of the data corresponding to each regime.
- `regimes::UnitRange`: which regimes are involved in the date range for which
    we want to compute the deterministic trends

where `S<:AbstractFloat`.

### Outputs

- `states::Matrix{S}`: matrix of size `nstates` x `nperiods` of state
  steady-state values
- `obs::Matrix{S}`: matrix of size `nobs` x `nperiods` of observable
  steady-state values
- `pseudo::Matrix{S}`: matrix of size `npseudo` x `nperiods` of
  pseudo-observable steady-state values

where `nperiods` is the number of quarters between
`date_shockdec_start(m)` and `date_shockdec_end(m)`, inclusive.
"""
function deterministic_trends(m::AbstractDSGEModel{S},
                              system::System{S}, z0::Vector{S}) where {S<:AbstractFloat}

    # Dates: We compute the deterministic trend starting from the
    # first historical period.  However, since it is only used to
    # compute shock decompositions, we truncate and only store
    # results for periods corresponding to the shockdec period.
    nperiods    = subtract_quarters(date_forecast_end(m), date_mainsample_start(m)) + 1
    start_index = index_shockdec_start(m)
    end_index   = index_shockdec_end(m)

    deterministic_trends(system, z0, nperiods, start_index, end_index)
end

function deterministic_trends(system::System{S}, z0::Vector{S}, nperiods::Int,
                              start_index::Int, end_index::Int) where {S<:AbstractFloat}

    # Set constant system matrices to 0
    system = zero_system_constants(system)

    # Construct matrix of 0 shocks for entire history and forecast horizon
    nshocks  = size(system[:RRR], 2)
    shocks   = zeros(S, nshocks, nperiods)

    # Use forecast to iterate state-space system forward without shocks or the ZLB procedure
    states, obs, pseudo, _ = forecast(system, z0, shocks)

    # Return deterministic trends in appropriate range
    if start_index == 1 && end_index == nperiods
        return states, obs, pseudo
    else
        range = start_index:end_index
        return states[:, range], obs[:, range], pseudo[:, range]
    end
end

function deterministic_trends(m::AbstractDSGEModel{S},
                              system::RegimeSwitchingSystem{S}, z0::Vector{S},
                              start_date::Dates.Date = date_presample_start(m),
                              end_date::Dates.Date = prev_quarter(date_forecast_start(m))) where {S<:AbstractFloat}

    # Dates: We compute the deterministic trend starting from the
    # first historical period.  However, since it is only used to
    # compute shock decompositions, we truncate and only store
    # results for periods corresponding to the shockdec period.
    nperiods    = subtract_quarters(date_forecast_end(m), date_mainsample_start(m)) + 1
    start_index = index_shockdec_start(m)
    end_index   = index_shockdec_end(m)

    # Now intersect these periods with regime indices
    regime_inds = regime_indices(m, start_date, end_date) # do not need to account for ZLB split b/c shocks don't matter for trends
    if regime_inds[1][1] < 1
        regime_inds[1] = 1:regime_inds[1][end]
    end
    regime_inds[end] = regime_inds[end][1]:end_index      # if the end index is in the middle of a regime or is past the regime's end

    return deterministic_trends(system, z0, nperiods, start_index, end_index, regime_inds)
end

function deterministic_trends(system::RegimeSwitchingSystem{S}, z0::Vector{S}, nperiods::Int,
                              start_index::Int, end_index::Int,
                              regime_inds::Vector{UnitRange{Int}}) where {S<:AbstractFloat}

    # Set constant system matrices to 0
    system = zero_system_constants(system)

    # Initialize storage matrices
    nshocks = size(system[1, :RRR], 2) # this is used in the forecasting loops to create shock matrices
    nstates = size(system[1, :TTT], 1)
    nobs    = size(system[1, :ZZ], 1)
    npseudo = size(system[1, :ZZ_pseudo], 1)
    states = Matrix{S}(undef, nstates, nperiods)
    obs    = Matrix{S}(undef, nobs,    nperiods)
    pseudo = Matrix{S}(undef, npseudo, nperiods)

    # Use forecast to iterate state-space system forward without shocks or the ZLB procedure
    for (reg_num, reg_ind) in enumerate(regime_inds)
        init_state = (reg_num == 1) ? z0 : states[:, reg_ind[1] - 1] # update initial state

        # Construct matrix of 0 shocks for this regime
        shocks = zeros(S, nshocks, length(reg_ind))
        states[:, reg_ind], obs[:, reg_ind], pseudo[:, reg_ind], _ = forecast(system[reg_num], init_state, shocks)
    end
    if regime_inds[end][end] < nperiods
        shocks   = zeros(S, nshocks, nperiods - regime_inds[end][end])
        last_ind = (1 + regime_inds[end][end]):nperiods
        states[:, last_ind], obs[:, last_ind], pseudo[:, last_ind], _ = forecast(system[length(regime_inds)], states[:, regime_inds[end][end]], shocks)
    end

    if start_index == 1 && end_index == nperiods
        return states, obs, pseudo
    else
        range = start_index:end_index
        return states[:, range], obs[:, range], pseudo[:, range]
    end
end

"""
```
trends(system::System{S}) where {S<:AbstractFloat}
trends(system::RegimeSwitchingSystem{S}) where {S<:AbstractFloat}
```

Compute trend (steady-state) states, observables, and pseudo-observables. The
trend is used for plotting shock decompositions.
"""
function trends(system::System{S}) where {S<:AbstractFloat}

    state_trend  = system[:CCC]
    obs_trend    = system[:ZZ]*system[:CCC] + system[:DD]
    pseudo_trend = system[:ZZ_pseudo]*system[:CCC] + system[:DD_pseudo]

    return state_trend, obs_trend, pseudo_trend
end

function trends(system::RegimeSwitchingSystem{S}) where {S<:AbstractFloat}

    nreg  = n_regimes(system)
    input = 1:nreg

    state_trends  = Matrix{S}(undef, length(system[1, :CCC]), nreg)
    obs_trends    = Matrix{S}(undef, length(system[1, :DD]), nreg)
    pseudo_trends = Matrix{S}(undef, length(system[1, :DD_pseudo]), nreg)
    for i in input
        state_trends[:, i]  = system[i, :CCC]
        obs_trends[:, i]    = system[i, :ZZ] * system[i, :CCC] + system[i, :DD]
        pseudo_trends[:, i] = system[i, :ZZ_pseudo] * system[i, :CCC] + system[i, :DD_pseudo]
    end

    return state_trends, obs_trends, pseudo_trends
end


text:
using Test
using NiLang, NiLang.AD

@testset "integrate" begin
    FT = Float64
    n = 100
    h = FT(π/n)
    dt = FT(0.01)
    α = FT(4e-2)
    @i function step!(dest::AbstractArray{T}, src::AbstractArray{T}; α, h, dt) where T
        n ← length(dest)
        @invcheckoff for i=1:n
            @routine begin
                @zeros T cum g h2 αcum
                cum += src[mod1(i+1, n)] + src[mod1(i-1, n)]
                cum -= 2*src[i]
                αcum += cum * α
                h2 += h^2
                g += αcum/h2
            end
            dest[i] += src[i]
            dest[i] += dt*g
            ~@routine
        end
        n → length(dest)
    end
    x = zeros(FT, n)
    x[n÷2] = 1
    #state = Dict{Int,Vector{FT}}()
    k = 4
    N = 100
    x_last = NiLang.direct_emulate(step!, FT.(x); N=N, α=α, h=h, dt=dt)
    log1 = NiLang.BennettLog()
    log2 = NiLang.BennettLog()
    log3 = NiLang.BennettLog()
    _, x_last_b, _ = bennett(step!, zero(FT.(x)), FT.(x); k=k, N=N, α=α, h=h, dt=dt, logger=log1)
    _, x_last_b2 = bennett!(step!, Dict(1=>FT.(x)); k=k, N=N, α=α, h=h, dt=dt, logger=log2)
    _, x_last_b3 = bennett!([step! for _=1:100], Dict(1=>FT.(x)); k=k, N=N, α=α, h=h, dt=dt, logger=log3)
    @test sum(x_last_b) ≈ 1
    @test x_last ≈ x_last_b
    @test x_last ≈ x_last_b2[N+1]
    @test x_last ≈ x_last_b3[N+1]
    @test length(log1.fcalls) > length(log2.fcalls)
    @test length(log1.fcalls) < 2*length(log2.fcalls)
    @test length(log3.fcalls) == length(log2.fcalls)

    @i function loss(out, step, y, x; kwargs...)
        bennett((@skip! step), y, x; kwargs...)
        out += y[n÷2]
    end
    @i function loss2(out, step, d; N, kwargs...)
        bennett!((@skip! step), d; N, kwargs...)
        out += d[N+1][n÷2]
    end
    _, _, _, gx = NiLang.AD.gradient(loss, (0.0, step!, zero(x), copy(x)); iloss=1, k=k, N=N, α=α, h=h, dt=dt)
    _, _, gx2 = NiLang.AD.gradient(loss2, (0.0, step!, Dict(1=>copy(x))); iloss=1, k=k, N=N, α=α, h=h, dt=dt)
    x_last_2 = NiLang.direct_emulate(step!, (x2=copy(x); x2[n÷2]+=1e-5; FT.(x2)); N=N, α=α, h=h, dt=dt)
    @test gx[n÷2] ≈ (x_last_2 - x_last)[n÷2]/1e-5
    @test gx2[1][n÷2] ≈ (x_last_2 - x_last)[n÷2]/1e-5
end

ext: jl


text:
function generate_raw_moment_eqs(drift_eqs::AbstractVector{Equation}, diff::AbstractArray,
                                 m_order::Int, vars, name, ps, iv)

    N = length(drift_eqs)
    drift = [e.rhs for e in drift_eqs]
    diff_mat = diff*transpose(diff)

    smap = Dict(Pair.(vars, 1:length(vars)))
    drift_factors, drift_powers, drift_order = polynomial_propensities(drift, iv, smap)
    diff_factors, diff_powers, diff_order = polynomial_propensities(value.(diff_mat), iv, smap)
    q_order = m_order + max(drift_order - 1, diff_order - 2)

    iter_all = construct_iter_all(N, q_order)
    iter_m = filter(x -> 1 < sum(x) <= m_order, iter_all)
    iter_q = filter(x -> m_order < sum(x) <= q_order, iter_all)
    iter_1 = filter(x -> sum(x) == 1, iter_all)

    μ = define_μ(iter_all, iv)
    eqs = Equation[]

    for i in 1:N
        poly = sum( drift_factors[i][j] * μ[Tuple(drift_powers[i][j])] for j in 1:length(drift_factors[i]) )
        push!(eqs, Differential(iv)(μ[iter_1[i]]) ~ expand(poly))
    end

    for iter in iter_m

        term1 = 0
        for i in 1:N
            iszero(iter[i]) ? continue :
            gradf = iter[i]
            gradh = iter .- iter_1[i]
            term1 += gradf*sum( μ[gradh .+ Tuple(drift_powers[i][j])] * drift_factors[i][j]
                                for j in 1:length(drift_factors[i]) )
        end

        term2 = 0
        for i in 1:N
            for j in 1:N
                ( i == j && isone(iter[i]) ) || (iszero(iter[i]) || iszero(iter[j])) ? continue :
                hji = iter .- (iter_1[i] .+ iter_1[j])
                hji_factor = i == j ? iter[i]*(iter[i]-1) : iter[i]*iter[j]
                term2 += hji_factor * sum( μ[hji .+ Tuple(diff_powers[i, j][k])] * diff_factors[i, j][k]
                                                for k in 1:length(diff_factors[i, j]))
            end
        end

        push!(eqs, Differential(iv)(μ[iter]) ~ expand(term1 + term2/2))
    end

    odename = Symbol(name, "_raw_moment_eqs_m", m_order)
    odes = ODESystem(eqs, iv, extract_variables(eqs, μ), ps; name=odename)

    RawMomentEquations(odes, μ, N, m_order, q_order, iter_all, iter_m, iter_q, iter_1)
end

generate_raw_moment_eqs(sys::SDESystem, m_order::Int) = generate_raw_moment_eqs(equations(sys), get_noiseeqs(sys), m_order,
                                                                                states(sys), nameof(sys), parameters(sys), get_iv(sys))

#=
function generate_raw_moment_eqs(drift_eqs::AbstractVector{Equation}, diff::AbstractArray{T},
                                 m_order::Int, vars, name, ps, iv) where T <: Union{Sym, Num}

    N = length(drift_eqs)
    drift = [e.rhs for e in drift_eqs]
    diff_mat = diff*transpose(diff)

    q_order = m_order + max(maximum(degree(eq, vars) for eq in drift) - 1, maximum(degree(eq, vars) for eq in diff_mat) - 2)
    iter_all = construct_iter_all(N, q_order)
    iter_m = filter(x -> 1 < sum(x) <= m_order, iter_all)
    iter_q = filter(x -> m_order < sum(x) <= q_order, iter_all)
    iter_1 = filter(x -> sum(x) == 1, iter_all)

    μ = define_μ(iter_all, iv)
    mono_to_moment = OrderedDict(prod(vars.^iter) => μ[iter] for iter in iter_all)
    mom_eqs = Equation[]
    moms = []
    for iter in vcat(iter_1, iter_m)
        mono = prod(vars.^iter)
        H = hessian(mono, vars)
        poly = sum(drift .* gradient(mono, vars)) + 1/2*sum(diff_mat[i,:]'*H[:,i] for i in 1:N )
        push!(mom_eqs, Differential(iv)(mono_to_moment[mono]) ~ poly_subs(poly, mono_to_moment, ps, true))
    end

    odename = Symbol(name, "_raw_moment_eqs_m", m_order)
    odes = ODESystem(mom_eqs, iv, moms, ps; name=odename)

    RawMomentEquations(odes, μ, N, m_order, q_order, iter_all, iter_m, iter_q, iter_1)
end

=#


text:
module TestClient

using HTTP
include(joinpath(dirname(pathof(HTTP)), "../test/resources/TestRequest.jl"))
using .TestRequest
using Sockets
using JSON
using Test
using URIs

status(r) = r.status
@testset "Custom HTTP Stack" begin
   @testset "Low-level Request" begin
        wasincluded = Ref(false)
        result = TestRequest.get("https://httpbin.org/ip"; httptestlayer=wasincluded)
        @test status(result) == 200
        @test wasincluded[]
    end
end

@testset "Client.jl - $sch" for sch in ["http", "https"]
    @testset "GET, HEAD, POST, PUT, DELETE, PATCH" begin
        @test status(HTTP.get("$sch://httpbin.org/ip")) == 200
        @test status(HTTP.head("$sch://httpbin.org/ip")) == 200
        @test status(HTTP.post("$sch://httpbin.org/ip"; status_exception=false)) == 405
        @test status(HTTP.post("$sch://httpbin.org/post")) == 200
        @test status(HTTP.put("$sch://httpbin.org/put")) == 200
        @test status(HTTP.delete("$sch://httpbin.org/delete")) == 200
        @test status(HTTP.patch("$sch://httpbin.org/patch")) == 200
    end

    @testset "ASync Client Requests" begin
        @test status(fetch(@async HTTP.get("$sch://httpbin.org/ip"))) == 200
        @test status(HTTP.get("$sch://httpbin.org/encoding/utf8")) == 200
    end

    @testset "Query to URI" begin
        r = HTTP.get(URI(HTTP.URI("$sch://httpbin.org/response-headers"); query=Dict("hey"=>"dude")))
        h = Dict(r.headers)
        @test (haskey(h, "Hey") ? h["Hey"] == "dude" : h["hey"] == "dude")
    end

    @testset "Cookie Requests" begin
        empty!(HTTP.access_threaded(Dict{String, Set{HTTP.Cookie}}, HTTP.CookieRequest.default_cookiejar))
        r = HTTP.get("$sch://httpbin.org/cookies", cookies=true)

        body = String(r.body)
        @test replace(replace(body, " "=>""), "\n"=>"")  == "{\"cookies\":{}}"

        r = HTTP.get("$sch://httpbin.org/cookies/set?hey=sailor&foo=bar", cookies=true)
        @test status(r) == 200

        body = String(r.body)
        @test replace(replace(body, " "=>""), "\n"=>"")  == "{\"cookies\":{\"foo\":\"bar\",\"hey\":\"sailor\"}}"

        r = HTTP.get("$sch://httpbin.org/cookies/delete?hey")
        cookies = JSON.parse(String(r.body))["cookies"]
        @test length(cookies) == 1
        @test cookies["foo"] == "bar"
    end

    @testset "Client Streaming Test" begin
        r = HTTP.post("$sch://httpbin.org/post"; body="hey")
        @test status(r) == 200

        # stream, but body is too small to actually stream
        r = HTTP.post("$sch://httpbin.org/post"; body="hey", stream=true)
        @test status(r) == 200

        r = HTTP.get("$sch://httpbin.org/stream/100")
        @test status(r) == 200

        bytes = r.body
        a = [JSON.parse(l) for l in split(chomp(String(bytes)), "\n")]
        totallen = length(bytes) # number of bytes to expect

        io = IOBuffer()
        r = HTTP.get("$sch://httpbin.org/stream/100"; response_stream=io)
        seekstart(io)
        @test status(r) == 200

        b = [JSON.parse(l) for l in eachline(io)]
        @test all(zip(a, b)) do (x, y)
            x["args"] == y["args"] &&
            x["id"] == y["id"] &&
            x["url"] == y["url"] &&
            x["origin"] == y["origin"] &&
            x["headers"]["Content-Length"] == y["headers"]["Content-Length"] &&
            x["headers"]["Host"] == y["headers"]["Host"] &&
            x["headers"]["User-Agent"] == y["headers"]["User-Agent"]
        end
    end

    @testset "Client Body Posting - Vector{UTF8}, String, IOStream, IOBuffer, BufferStream" begin
        @test status(HTTP.post("$sch://httpbin.org/post"; body="hey")) == 200
        @test status(HTTP.post("$sch://httpbin.org/post"; body=UInt8['h','e','y'])) == 200
        io = IOBuffer("hey"); seekstart(io)
        @test status(HTTP.post("$sch://httpbin.org/post"; body=io)) == 200
        tmp = tempname()
        open(f->write(f, "hey"), tmp, "w")
        io = open(tmp)
        @test status(HTTP.post("$sch://httpbin.org/post"; body=io, enablechunked=false)) == 200
        close(io); rm(tmp)
        f = Base.BufferStream()
        write(f, "hey")
        close(f)
        @test status(HTTP.post("$sch://httpbin.org/post"; body=f, enablechunked=false)) == 200
    end

    @testset "Chunksize" begin
        #     https://github.com/JuliaWeb/HTTP.jl/issues/60
        #     Currently httpbin.org responds with 411 status and “Length Required”
        #     message to any POST/PUT requests that are sent using chunked encoding
        #     See https://github.com/kennethreitz/httpbin/issues/340#issuecomment-330176449
        println("client transfer-encoding chunked")
        @test status(HTTP.post("$sch://httpbin.org/post"; body="hey", #=chunksize=2=#)) == 200
        @test status(HTTP.post("$sch://httpbin.org/post"; body=UInt8['h','e','y'], #=chunksize=2=#)) == 200
        io = IOBuffer("hey"); seekstart(io)
        @test status(HTTP.post("$sch://httpbin.org/post"; body=io, #=chunksize=2=#)) == 200
        tmp = tempname()
        open(f->write(f, "hey"), tmp, "w")
        io = open(tmp)
        @test status(HTTP.post("$sch://httpbin.org/post"; body=io, #=chunksize=2=#)) == 200
        close(io); rm(tmp)
        f = Base.BufferStream()
        write(f, "hey")
        close(f)
        @test status(HTTP.post("$sch://httpbin.org/post"; body=f, #=chunksize=2=#)) == 200
    end

    @testset "Incomplete response with known content length" begin
        server = Sockets.listen(ip"0.0.0.0", 8080)
        try
            task = @async HTTP.listen("0.0.0.0", 8080; server=server) do http
                HTTP.setstatus(http, 200)
                HTTP.setheader(http, "Content-Length" => "64") # Promise 64 bytes...
                HTTP.startwrite(http)
                HTTP.write(http, rand(UInt8, 63)) # ...but only send 63 bytes.
                # Close the stream so that eof(stream) is true and the client isn't
                # waiting forever for the last byte.
                HTTP.close(http.stream)
            end

            err = try
                HTTP.get("http://localhost:8080"; retry=false)
            catch err
                err
            end
            @test err isa HTTP.IOError
            @test err.e isa EOFError

        finally
            # Shutdown
            try; close(server); wait(task); catch; end
            HTTP.ConnectionPool.closeall()
        end
    end

    @testset "ASync Client Request Body" begin
        f = Base.BufferStream()
        write(f, "hey")
        t = @async HTTP.post("$sch://httpbin.org/post"; body=f, enablechunked=false)
        #fetch(f) # fetch for the async call to write it's first data
        write(f, " there ") # as we write to f, it triggers another chunk to be sent in our async request
        write(f, "sailor")
        close(f) # setting eof on f causes the async request to send a final chunk and return the response
        @test status(fetch(t)) == 200
    end

    @testset "Client Redirect Following - $read_method" for read_method in ["GET", "HEAD"]
        @test_skip status(HTTP.request(read_method, "$sch://httpbin.org/redirect/1")) ==200
        @test_skip status(HTTP.request(read_method, "$sch://httpbin.org/redirect/1", redirect=false)) == 302
        @test_skip status(HTTP.request(read_method, "$sch://httpbin.org/redirect/6")) == 302 #over max number of redirects
        @test_skip status(HTTP.request(read_method, "$sch://httpbin.org/relative-redirect/1")) == 200
        @test_skip status(HTTP.request(read_method, "$sch://httpbin.org/absolute-redirect/1")) == 200
        @test_skip status(HTTP.request(read_method, "$sch://httpbin.org/redirect-to?url=http%3A%2F%2Fgoogle.com")) == 200
    end

    @testset "Client Basic Auth" begin
        @test status(HTTP.get("$sch://user:pwd@httpbin.org/basic-auth/user/pwd")) == 200
        @test status(HTTP.get("$sch://user:pwd@httpbin.org/hidden-basic-auth/user/pwd")) == 200
        @test status(HTTP.get("$sch://test:%40test@httpbin.org/basic-auth/test/%40test")) == 200
    end

    @testset "Misc" begin
        @test status(HTTP.post("$sch://httpbin.org/post"; body="√")) == 200
        r = HTTP.request("GET", "$sch://httpbin.org/ip")
        @test status(r) == 200

        uri = HTTP.URI("$sch://httpbin.org/ip")
        r = HTTP.request("GET", uri)
        @test status(r) == 200
        r = HTTP.get(uri)
        @test status(r) == 200

        r = HTTP.request("GET", "$sch://httpbin.org/ip")
        @test status(r) == 200

        uri = HTTP.URI("$sch://httpbin.org/ip")
        r = HTTP.request("GET", uri)
        @test status(r) == 200

        r = HTTP.get("$sch://httpbin.org/image/png")
        @test status(r) == 200

        # ensure we can use AbstractString for requests
        r = HTTP.get(SubString("http://httpbin.org/ip",1))

        # canonicalizeheaders
        @test status(HTTP.get("$sch://httpbin.org/ip"; canonicalizeheaders=false)) == 200
    end
end

if !isempty(get(ENV, "PIE_SOCKET_API_KEY", "")) && get(ENV, "JULIA_VERSION", "") == "1"
    println("found pie socket api key, running websocket tests")
    pie_socket_api_key = ENV["PIE_SOCKET_API_KEY"]
    @testset "openraw client method - $socket_protocol" for socket_protocol in ["wss", "ws"]
        # WebSockets require valid headers.
        headers = Dict(
            "Upgrade" => "websocket",
            "Connection" => "Upgrade",
            "Sec-WebSocket-Key" => "dGhlIHNhbXBsZSBub25jZQ==",
            "Sec-WebSocket-Version" => "13")

        socket, response = HTTP.openraw("GET", "$socket_protocol://free3.piesocket.com/v3/http_test_channel?api_key=$pie_socket_api_key&notify_self", headers)

        @test response.status == 101

        # This is an example text frame from RFC 6455, section 5.7. It sends the text "Hello" to the
        # echo server, and so we expect "Hello" back, in an unmasked frame.
        frame = UInt8[0x81, 0x85, 0x37, 0xfa, 0x21, 0x3d, 0x7f, 0x9f, 0x4d, 0x51, 0x58]

        write(socket, frame)

        # The frame we expect back looks like `expectedframe`.
        expectedframe = UInt8[0x81, 0x05, 0x48, 0x65, 0x6c, 0x6c, 0x6f]

        # Note the spec for read says:
        #     read(s::IO, nb=typemax(Int))
        # Read at most nb bytes from s, returning a Vector{UInt8} of the bytes read.
        # ... so read will return less than 7 bytes unless we wait first:
        eof(socket)
        actualframe = read(socket, 7)
        @test expectedframe == actualframe

        close(socket)
    end
end

@testset "HTTP.open accepts method::Symbol" begin
    @test status(HTTP.open(x -> x, :GET, "http://httpbin.org/ip")) == 200
end

@testset "readtimeout" begin
    @test_throws HTTP.TimeoutRequest.ReadTimeoutError begin
        HTTP.get("http://httpbin.org/delay/5"; readtimeout=1, retry=false)
    end
    HTTP.get("http://httpbin.org/delay/1"; readtimeout=2, retry=false)
end

@testset "Retry all resolved IP addresses" begin
    # See issue https://github.com/JuliaWeb/HTTP.jl/issues/672
    # Bit tricky to test, but can at least be tested if localhost
    # resolves to both IPv4 and IPv6 by listening to the respective
    # interface
    alladdrs = getalladdrinfo("localhost")
    if ip"127.0.0.1" in alladdrs && ip"::1" in alladdrs
        for interface in (IPv4(0), IPv6(0))
            server = listen(interface, 8080)
            try
                @async HTTP.listen(string(interface), 8080; server=server) do http
                    HTTP.setstatus(http, 200)
                    HTTP.startwrite(http)
                    HTTP.write(http, "hello, world")
                end
                req = HTTP.get("http://localhost:8080")
                @test req.status == 200
                @test String(req.body) == "hello, world"
            finally
                close(server)
                HTTP.ConnectionPool.closeall()
            end
        end
    end
end

@testset "Sockets.get(sock|peer)name(::HTTP.Stream)" begin
    server = listen(IPv4(0), 8080)
    try
        @async HTTP.listen("0.0.0.0", 8080; server=server) do http
            sock = Sockets.getsockname(http)
            peer = Sockets.getpeername(http)
            str = sprint() do io
                print(io, sock[1], ":", sock[2], " - ", peer[1], ":", peer[2])
            end
            HTTP.setstatus(http, 200)
            HTTP.setheader(http, "Content-Length" => string(sizeof(str)))
            HTTP.startwrite(http)
            HTTP.write(http, str)
        end

        # Tests for Stream{TCPSocket}
        HTTP.open("GET", "http://localhost:8080") do http
            # Test server peer/sock
            reg = r"^127\.0\.0\.1:8080 - 127\.0\.0\.1:(\d+)$"
            m = match(reg, read(http, String))
            @test m !== nothing
            server_peerport = parse(Int, m[1])
            # Test client peer/sock
            sock = Sockets.getsockname(http)
            @test sock[1] == ip"127.0.0.1"
            @test sock[2] == server_peerport
            peer = Sockets.getpeername(http)
            @test peer[1] == ip"127.0.0.1"
            @test peer[2] == 8080
        end
    finally
        close(server)
        HTTP.ConnectionPool.closeall()
    end

    # Tests for Stream{SSLContext}
    HTTP.open("GET", "https://julialang.org") do http
        sock = Sockets.getsockname(http)
        if VERSION >= v"1.2.0"
            @test sock[1] in Sockets.getipaddrs()
        end
        peer = Sockets.getpeername(http)
        @test peer[1] in Sockets.getalladdrinfo("julialang.org")
        @test peer[2] == 443
    end
end

@testset "input verification of bad URLs" begin
    # HTTP.jl#527, HTTP.jl#545
    url = "julialang.org"
    @test_throws ArgumentError("missing or unsupported scheme in URL (expected http(s) or ws(s)): $(url)") HTTP.get(url)
    url = "ptth://julialang.org"
    @test_throws ArgumentError("missing or unsupported scheme in URL (expected http(s) or ws(s)): $(url)") HTTP.get(url)
    url = "http:julialang.org"
    @test_throws ArgumentError("missing host in URL: $(url)") HTTP.get(url)
end

@testset "Implicit request headers" begin
    server = listen(IPv4(0), 8080)
    try
        tsk = @async HTTP.listen("0.0.0.0", 8080; server=server) do http
            data = Dict{String,String}(http.message.headers)
            HTTP.setstatus(http, 200)
            HTTP.startwrite(http)
            HTTP.write(http, sprint(JSON.print, data))
        end
        old_user_agent = HTTP.DefaultHeadersRequest.USER_AGENT[]
        default_user_agent = "HTTP.jl/$VERSION"
        # Default values
        HTTP.setuseragent!(default_user_agent)
        d = JSON.parse(IOBuffer(HTTP.get("http://localhost:8080").body))
        @test d["Host"] == "localhost:8080"
        @test d["Accept"] == "*/*"
        @test d["User-Agent"] == default_user_agent
        # Overwriting behavior
        headers = ["Host" => "http.jl", "Accept" => "application/json"]
        HTTP.setuseragent!("HTTP.jl test")
        d = JSON.parse(IOBuffer(HTTP.get("http://localhost:8080", headers).body))
        @test d["Host"] == "http.jl"
        @test d["Accept"] == "application/json"
        @test d["User-Agent"] == "HTTP.jl test"
        # No User-Agent
        HTTP.setuseragent!(nothing)
        d = JSON.parse(IOBuffer(HTTP.get("http://localhost:8080").body))
        @test !haskey(d, "User-Agent")

        HTTP.setuseragent!(old_user_agent)
    finally
        close(server)
        HTTP.ConnectionPool.closeall()
    end
end

import NetworkOptions, MbedTLS
@testset "NetworkOptions for host verification" begin
    # Set up server with self-signed cert
    server = listen(IPv4(0), 8443)
    try
        cert, key = joinpath.(dirname(pathof(HTTP)), "../test", "resources", ("cert.pem", "key.pem"))
        sslconfig = MbedTLS.SSLConfig(cert, key)
        tsk = @async HTTP.listen("0.0.0.0", 8443; server=server, sslconfig=sslconfig) do http
            HTTP.setstatus(http, 200)
            HTTP.startwrite(http)
            HTTP.write(http, "hello, world")
        end
        url = "https://localhost:8443"
        env = ["JULIA_NO_VERIFY_HOSTS" => nothing, "JULIA_SSL_NO_VERIFY_HOSTS" => nothing, "JULIA_ALWAYS_VERIFY_HOSTS" => nothing]
        withenv(env...) do
            @test NetworkOptions.verify_host(url)
            @test NetworkOptions.verify_host(url, "SSL")
            @test_throws HTTP.IOError HTTP.get(url; retries=1)
            @test_throws HTTP.IOError HTTP.get(url; require_ssl_verification=true, retries=1)
            @test HTTP.get(url; require_ssl_verification=false).status == 200
        end
        withenv(env..., "JULIA_NO_VERIFY_HOSTS" => "localhost") do
            @test !NetworkOptions.verify_host(url)
            @test !NetworkOptions.verify_host(url, "SSL")
            @test HTTP.get(url).status == 200
            @test_throws HTTP.IOError HTTP.get(url; require_ssl_verification=true, retries=1)
            @test HTTP.get(url; require_ssl_verification=false).status == 200
        end
        withenv(env..., "JULIA_SSL_NO_VERIFY_HOSTS" => "localhost") do
            @test NetworkOptions.verify_host(url)
            @test !NetworkOptions.verify_host(url, "SSL")
            @test HTTP.get(url).status == 200
            @test_throws HTTP.IOError HTTP.get(url; require_ssl_verification=true, retries=1)
            @test HTTP.get(url; require_ssl_verification=false).status == 200
        end
    finally
        close(server)
        HTTP.ConnectionPool.closeall()
    end
end

@testset "Public entry point of HTTP.request and friends (e.g. issue #463)" begin
    headers = Dict("User-Agent" => "HTTP.jl")
    query = Dict("hello" => "world")
    body = UInt8[1, 2, 3]
    stack = HTTP.stack()
    function test(r, m)
        @test r.status == 200
        d = JSON.parse(IOBuffer(HTTP.payload(r)))
        @test d["headers"]["User-Agent"] == "HTTP.jl"
        @test d["data"] == "\x01\x02\x03"
        @test endswith(d["url"], "?hello=world")
        @test d["method"] == m
    end
    for uri in ("https://httpbin.org/anything", HTTP.URI("https://httpbin.org/anything"))
        # HTTP.request
        test(HTTP.request("GET", uri; headers=headers, body=body, query=query), "GET")
        test(HTTP.request("GET", uri, headers; body=body, query=query), "GET")
        test(HTTP.request("GET", uri, headers, body; query=query), "GET")
        !isa(uri, HTTP.URI) && test(HTTP.request(stack, "GET", uri; headers=headers, body=body, query=query), "GET")
        test(HTTP.request(stack, "GET", uri, headers; body=body, query=query), "GET")
        test(HTTP.request(stack, "GET", uri, headers, body; query=query), "GET")
        # HTTP.get
        test(HTTP.get(uri; headers=headers, body=body, query=query), "GET")
        test(HTTP.get(uri, headers; body=body, query=query), "GET")
        test(HTTP.get(uri, headers, body; query=query), "GET")
        # HTTP.put
        test(HTTP.put(uri; headers=headers, body=body, query=query), "PUT")
        test(HTTP.put(uri, headers; body=body, query=query), "PUT")
        test(HTTP.put(uri, headers, body; query=query), "PUT")
        # HTTP.post
        test(HTTP.post(uri; headers=headers, body=body, query=query), "POST")
        test(HTTP.post(uri, headers; body=body, query=query), "POST")
        test(HTTP.post(uri, headers, body; query=query), "POST")
        # HTTP.patch
        test(HTTP.patch(uri; headers=headers, body=body, query=query), "PATCH")
        test(HTTP.patch(uri, headers; body=body, query=query), "PATCH")
        test(HTTP.patch(uri, headers, body; query=query), "PATCH")
        # HTTP.delete
        test(HTTP.delete(uri; headers=headers, body=body, query=query), "DELETE")
        test(HTTP.delete(uri, headers; body=body, query=query), "DELETE")
        test(HTTP.delete(uri, headers, body; query=query), "DELETE")
    end
end

@testset "HTTP CONNECT Proxy" begin
    @testset "Host header" begin
        # Stores the http request passed by the client
        req = String[]

        # Trivial implementation of a proxy server
        # We are only interested in the request passed in by the client
        # Returns 400 after reading the http request into req
        proxy = listen(IPv4(0), 8082)
        try
            @async begin
                sock = accept(proxy)
                while isopen(sock)
                    line = readline(sock)
                    isempty(line) && break

                    push!(req, line)
                end
                write(sock, "HTTP/1.1 400 Bad Request\r\n\r\n")
            end

            # Make the HTTP request
            HTTP.get("https://example.com"; proxy="http://localhost:8082", retry=false, status_exception=false)

            # Test if the host header exist in the request
            @test "Host: example.com:443" in req
        finally
            close(proxy)
            HTTP.ConnectionPool.closeall()
        end
    end
end

end # module


text:
import Pkg

neuralpde_dir = abspath(joinpath(@__DIR__, ".."))
@info "neuralpde_directory: $(neuralpde_dir)"

neuralpdelogging_dir = abspath(joinpath(neuralpde_dir, "lib", "NeuralPDELogging"))
@info "developing NeuralPDELogging subpackage at dir: $(neuralpdelogging_dir)"
neuralpdelogging = Pkg.PackageSpec(path = neuralpdelogging_dir)
@info "neuralpdelogging: $(string(neuralpdelogging))"

@info "testing NeuralPDELogging subpackage"
inner_is_CI = haskey(ENV,"CI")
@info "inner_is_CI: $(inner_is_CI)"
Pkg.activate(neuralpdelogging_dir)
@info Pkg.project()
Pkg.test(; coverage = inner_is_CI)


text:
struct ReLU{T, O} <: AbstractBlock
    input::O
    ReLU{T, O}() where {T, O} = new{T, O}()
end

ReLU(::Type{T}, ::Type{O}) where {T, O} = ReLU{T, O}()
ReLU(::Type{T};nbatch=1) where T = ReLU(T, nbatch > 1 ? Matrix{T} : Vector{T})
ReLU(;nbatch=1) = ReLU(Float64; nbatch=nbatch)

_relu(x::T) where {T <: Real} = x > 0 ? x : 0
_relu_grad(x::T, grad::T) where {T <: Real} = x > 0 ? grad : 0

function _relu(x::T) where {T <: Complex}
    if real(x) > 0 && imag(x) > 0
        return x
    elseif real(x) > 0 && imag(x) < 0
        return real(x)
    elseif real(x) < 0 && imag(x) > 0
        return imag(x)
    else
        return 0
    end
end

function _relu_grad(x::T, grad::T) where {T <: Complex}
    if real(x) > 0 && imag(x) > 0
        return grad
    elseif real(x) > 0 && imag(x) < 0
        return real(grad)
    elseif real(x) < 0 && imag(x) > 0
        return -im * real(grad)
    else
        return 0
    end
end

function forward(op::ReLU{T, O}, input::O) where {T, O}
    return _relu.(input)
end

function backward(op::ReLU{T, O}, grad::O) where {T, O}
    return _relu_grad.(op.input, grad)
end


text:
using PyCall
using PyCallUtils
using SparseArrays
using BSON
using Test
using PyCall: python

@imports numpy as np
@test np == pyimport("numpy")

@from datetime imports datetime
@test datetime == pyimport("datetime").datetime

if occursin("scipy", read(`$python -m pip list`, String))
    x = sprand(100, 100, 0.5)
    @test convert(SparseMatrixCSC, PyObject(x)) == x
end

@test pytypename(np.sin) == "ufunc"
@test pymodulename(np.sin) == "numpy"

BSON.@save "test.bson" x = np.sin
BSON.@load "test.bson" x
@test x == np.sin
rm("test.bson")

@test (np.sin ∘ np.cos)(1f0) == np.sin(np.cos(1f0))

if occursin("pandas", read(`$python -m pip list`, String))
    @imports pandas as pd

    df = pd.DataFrame(Dict(:name => ["a", "b"], :age => [27, 30]))
    @test df.loc[1, "age"] == 30

    sr = pd.Series([3, 5], index = [:a, :b])
    @test sr.a == 3
    @test sr["a"] == 3
    @test sr.loc["a"] == 3
    @test sr.b == 5
    @test sr.iloc[0] == 3
    @test sr.iloc[1] == 5
    @test length(sr) == 2
    @test sr.sum() == 8
    @test all(sr.iloc[0:1] == sr)
    @test all(sr.iloc[:] == sr)
end


text:
# auxiliary types and functions for Newton's method
# linear operator, newton data etc.

@doc """
  This type holds the data required by [`newtonInner`](@ref) as well as
  configuration settings.

  **Public Fields**

   * myrank: MPI rank
   * commsize: MPI communicator size
   * itr: number of iterations
   * res_norm_i: current iteration residual norm
   * res_norm_i_1: previous iteration residual norm
   * step_norm_i: current iteration newton step norm
   * step_norm_i_1: previous iteration newton step norm
   * res_norm_rel: norm of residual used as the reference point when computing
                   relative residuals.  If this is -1 on entry to newtonInner,
                   then the norm of the initial residual is used.
   * step_fac: factor used in step size limiter
   * res_reltol: nonlinear relative residual tolerance
   * res_abstol: nonlinear residual absolute tolerance
   * step_tol: step norm tolerance
   * itermax: maximum number of newton iterations
   * use_inexact_nk: true if inexact-NK should be used, false otherwise
   * krylov_gamma: parameter used by inexact newton-krylov
   * recalc_policy: a [`RecalculationPolicy`](@ref).
   * ls: a [`LinearSolver`](@ref)
   * fconv: convergence.dat file handle (or DevNull if not used)
   * verbose: how much logging/output to do

  **Options Keys**

  If `res_reltol0` is negative, the residual of the initial condition will be
  used for res_norm_rel

  `newton_verbosity` is used to the `verbose` field
"""->
mutable struct NewtonData{Tsol, Tres, Tsolver <: LinearSolver}

  # MPI info
  myrank::Int
  commsize::Int
  itr::Int  # newton iteration number

  # working variables
  res_norm_i::Float64  # current step residual norm
  res_norm_i_1::Float64  # previous step residual norm
  step_norm_i::Float64
  step_norm_i_1::Float64
  res_norm_rel::Float64  # norm of the residual used for the relative residual
                         # tolerance
  set_rel_norm::Bool # whether or not to use the initial residual as the
                     # res_norm_rel
  step_fac::Float64

  # tolerances (newton)
  res_reltol::Float64
  res_abstol::Float64
  step_tol::Float64
  itermax::Int

  # inexact Newton-Krylov parameters
  use_inexact_nk::Bool
  krylov_gamma::Float64  # update parameter for krylov tolerance

  recalc_policy::RecalculationPolicy
  ls::Tsolver
  res_0::Array{PetscScalar, 1}
  delta_q_vec::Array{PetscScalar, 1}
  fconv::IO  # convergence.dat, rank 0 only
  verbose::Int

end

#TODO: see if the static parameters are still needed
function NewtonData(mesh, sbp,
        eqn::AbstractSolutionData{Tsol, Tres}, opts,
        ls::LinearSolver) where {Tsol, Tres}

  myrank = mesh.myrank
  commsize = mesh.commsize
  itr = 0
  verbose = opts["newton_verbosity"]

  res_norm_i = 0.0
  res_norm_i_1 = 0.0
  step_norm_i = 0.0
  step_norm_i_1 = 0.0
  res_norm_rel = opts["res_reltol0"]
  set_rel_norm = res_norm_rel < 1
  step_fac = 1.0

  res_reltol = opts["res_reltol"]
  res_abstol = opts["res_abstol"]
  step_tol = opts["step_tol"]
  itermax = opts["itermax"]

  use_inexact_nk = opts["use_inexact_nk"]
  krylov_gamma = opts["krylov_gamma"]

  recalc_policy = getRecalculationPolicy(opts, "newton")

  # temporary vectors
  res_0 = zeros(PetscScalar, mesh.numDof)  # function evaluated at u0
  delta_q_vec = zeros(PetscScalar, mesh.numDof)  # newton update

  # convergence.dat
  if verbose >= 5
    if myrank == 0
      fconv = BufferedIO("convergence.dat", "a+")
    else
      fconv = DevNull
    end
  else
    fconv = DevNull
  end

  return NewtonData{Tsol, Tres, typeof(ls)}(myrank, commsize, itr,
                    res_norm_i, res_norm_i_1, step_norm_i, step_norm_i_1,
                    res_norm_rel, set_rel_norm, step_fac,
                    res_reltol, res_abstol, step_tol, itermax,
                    use_inexact_nk, krylov_gamma, recalc_policy, ls, res_0,
                    delta_q_vec, fconv, verbose)
end

@doc """
### NonlinearSolvers.setupNewton
  Performs setup work for [`newtonInner`](@ref), including creating a
  [`NewtonData`](@ref) object.

  This function also resets the implicit Euler globalization.

  alloc_rhs: keyword arg to allocate a new object or not for rhs_vec
                true (default) allocates a new vector
                false will use eqn.res_vec

  rhs_func: only used for Petsc in matrix-free mode to do Jac-vec products
            should be the rhs_func passed into [`newtonInner`](@ref)
  ctx_residual: ctx_residual passed into [`newtonInner`](@ref)

  Allocates Jac & RHS

  See [`cleanupNewton`](@ref) to the cleanup function

"""->
function setupNewton(mesh, pmesh, sbp,
         eqn::AbstractSolutionData{Tsol, Tres}, opts,
         ls::LinearSolver; alloc_rhs=true) where {Tsol, Tres}

  newton_data = NewtonData(mesh, sbp, eqn, opts, ls)

  clearEulerConstants(ls)
  # For simple cases, especially for Newton's method as a steady solver,
  #   having rhs_vec and eqn.res_vec pointing to the same memory
  #   saves us from having to copy back and forth
  if alloc_rhs
    rhs_vec = zeros(Tsol, size(eqn.res_vec))
  else
    rhs_vec = eqn.res_vec
  end

  return newton_data, rhs_vec

end   # end of setupNewton

"""
  Reinitialized the NewtonData object for a new solve.

  Note that this does not reset the linear solver, which might be a problem
  if inexact newton-krylov was used for the previous solve

  Uses `opts["newton_scale_euler"]` to determine how the `opts["euler_tau"]`
  value should be interpreted
"""
function reinitNewtonData(newton_data::NewtonData, mesh, sbp, eqn, opts)

  clearEulerConstants(newton_data.ls)
  newton_data.itr = 0
  newton_data.res_norm_i = 0
  newton_data.res_norm_i_1 = 0
  newton_data.step_norm_i = 0
  newton_data.step_norm_i_1 = 0
  newton_data.step_fac = 1.0
  resetRecalculationPolicy(newton_data.recalc_policy)

  # if using globalization, reset it
  if opts["newton_globalize_euler"]
    lo = getInnerLO(newton_data.ls.lo, NewtonLO)

    #TODO: it would be slightly better to update tau_vec inplace
    reinitImplicitEuler(mesh, opts, lo.idata)

    if !(typeof(newton_data.ls.pc) <: PCNone)
      # no need to do anything for PCNone
      pc = getInnerPC(newton_data.ls.pc, NewtonPC)
      reinitImplicitEuler(mesh, opts, lo.idata)
    end

    if opts["newton_scale_euler"]
      # in this mode, we say that the euler wants the pseudo time step to
      # be opts["euler_tau"] when the residual norm = 1 (therefore the residual
      # norm of the initial condition causes the time step to get scaled
      # up or down depending on the residual).
      recordEulerResidual(newton_data.ls, 1)
      # shift the residuals so they are in the i-1 position
      # This is caused by a bad interface to ImplicitEulerData:
      # recordEulerResidual should do this shifting, not useEulerConstants
      useEulerConstants(lo)
      if !(typeof(newton_data.ls.pc) <: PCNone)
        useEulerConstants(Pc)
      end
    end
  end

  return nothing
end

"""
  Cleans up after running Newton's method.

  **Inputs**

   * newton_data: the NewtonData object

"""
function free(newton_data::NewtonData)

  free(newton_data.ls)
  close(newton_data.fconv)

  return nothing
end

"""
  Records the most recent nonlinear residual norm in the NewtonData object.
  Also updates the implicit Euler globalization

  **Inputs**

   * newton_data: the NewtonData
   * res_norm: the residual norm
"""
function recordResNorm(newton_data::NewtonData, res_norm::Number)

  newton_data.res_norm_i_1 = newton_data.res_norm_i
  newton_data.res_norm_i = res_norm

  # update implicit Euler globalization
  recordEulerResidual(newton_data.ls, res_norm)

  return nothing
end

"""
  Records norm of the most recent newton step (ie. the norm of delta q)
  in the NewtonData object

  **Inputs**

   * newton_data: the NewtonData object
   * step_norm: the norm of the step
"""
function recordStepNorm(newton_data::NewtonData, step_norm::Number)

  newton_data.step_norm_i_1 = newton_data.step_norm_i
  newton_data.step_norm_i = step_norm

  return nothing
end

#------------------------------------------------------------------------------
# getter for PC and LO

"""
  Returns the Newton precondtioner and linear operator specified by the options
  dictionary

  **Inputs**

   * mesh
   * sbp
   * eqn
   * opts
   * rhs_func: rhs_func required by [`newtonInner`](@ref)
"""
function getNewtonPCandLO(mesh, sbp, eqn, opts,
                          jactype::Integer=opts["jac_type"])

  # get PC
  if jactype <= 2
    pc = PCNone(mesh, sbp, eqn, opts)
  else
    if opts["use_volume_preconditioner"]
      pc = NewtonBDiagPC(mesh, sbp, eqn, opts)
    else
      pc = NewtonMatPC(mesh, sbp, eqn, opts)
    end
  end

  if jactype == 1
    lo = NewtonDenseLO(pc, mesh, sbp, eqn, opts)
  elseif jactype == 2
    lo = NewtonSparseDirectLO(pc, mesh, sbp, eqn, opts)
  elseif jactype == 3
    lo = NewtonPetscMatLO(pc, mesh, sbp, eqn, opts)
  elseif jactype == 4
    lo = NewtonPetscMatFreeLO(pc, mesh, sbp, eqn, opts)
  end

  return pc, lo
end

import PDESolver.createLinearSolver

function createLinearSolver(mesh::AbstractMesh, sbp::AbstractOperator,
                            eqn::AbstractSolutionData, opts,
                            jac_type::Integer=opts["jac_type"])

  # general linear solvers should not have globilization
  val_orig = opts["setup_globalize_euler"]
  opts["setup_globalize_euler"] = false
  pc, lo = getNewtonPCandLO(mesh, sbp, eqn, opts, jac_type)
  ls = StandardLinearSolver(pc, lo, eqn.comm, opts)
  opts["setup_globalize_euler"] = val_orig

  return ls
end

#------------------------------------------------------------------------------
# preconditioner

abstract type NewtonMatFreePC <: AbstractPetscMatFreePC end

"""
  This function initializes the data needed to do Psudo-Transient Continuation
  globalization (aka. Implicit Euler) of Newton's method, using a spatially
  varying pseudo-timestep.

  Updates the jacobian with a diagonal term, as though the jac was the
  jacobian of this function:
  (u - u_i_1)/delta_t + f(u)
  where f is the original residual and u_i_1 is the previous step solution

  This globalization is activated using the option `newton_globalize_euler`.
  The initial value of the scaling factor tau is specified by the option
  `euler_tau`.

"""
mutable struct ImplicitEulerData
    use_implicit_euler::Bool  # whether or not use use implicit Euler
    res_norm_i::Float64  # current step residual norm
    res_norm_i_1::Float64  # previous step residual norm
    # Pseudo-transient continuation Euler
    tau_l::Float64  # current pseudo-timestep
    tau_vec::Array{Float64, 1}  # array of element-local time steps
end

"""
  Constructor for the case where implicit Euler will be used

  **Inputs**

   * mesh
   * opts
   * tau_l: the reference timestep
"""
function ImplicitEulerData(mesh::AbstractMesh, opts, tau_l::Number)

  use_implicit_euler = true
  res_norm_i = 0.0
  res_norm_i_1 = 0.0

  tau_vec = zeros(mesh.numDof)
  calcTauVec(mesh, opts, tau_l, tau_vec)

  return ImplicitEulerData(use_implicit_euler, res_norm_i, res_norm_i_1,
                           tau_l, tau_vec)
end

"""
  Constructor for the case where implicit Euler will not be used.
"""
function ImplicitEulerData()

  use_implicit_euler = false
  res_norm_i = 0.0
  res_norm_i_1 = 0.0
  tau_l = 0.0
  tau_vec = Array{Float64}(0)

  return ImplicitEulerData(use_implicit_euler, res_norm_i, res_norm_i_1,
                           tau_l, tau_vec)
end

"""
  Constructor from options dictionary.  Uses opts["setup_globalize_euler"]
  and opts["euler_tau"] to configure the returned object

  **Inputs**

   * mesh
   * opts
"""
function ImplicitEulerData(mesh::AbstractMesh, opts)

  if opts["setup_globalize_euler"]
    obj = ImplicitEulerData(mesh, opts, opts["euler_tau"])
  else
    obj = ImplicitEulerData()
  end

  return obj
end

"""
  Matrix-based Petsc preconditioner for Newton's method
"""
mutable struct NewtonMatPC <: AbstractPetscMatPC
  pc_inner::PetscMatPC
  idata::ImplicitEulerData
end

"""
  Outer constructor for [`NewtonMatPC`](@ref)
"""
function NewtonMatPC(mesh::AbstractMesh, sbp::AbstractOperator,
                    eqn::AbstractSolutionData, opts::Dict)

  pc_inner = PetscMatPC(mesh, sbp, eqn, opts)
  idata = ImplicitEulerData(mesh, opts)

  return NewtonMatPC(pc_inner, idata)
end

function calcPC(pc::NewtonMatPC, mesh::AbstractMesh, sbp::AbstractOperator,
                eqn::AbstractSolutionData, opts::Dict, ctx_residual, t)

  calcPC(pc.pc_inner, mesh, sbp, eqn, opts, ctx_residual, t)
  physicsJac(mesh, sbp, eqn, opts, getBasePC(pc).A, ctx_residual, t)

  if opts["newton_globalize_euler"]
    # TODO: updating the Euler parameter here is potentially wrong if we
    #       are not updating the Jacobian at every newton step
    updateEuler(pc)
    applyEuler(mesh, sbp, eqn, opts, pc)
  end

  return nothing
end

#------------------------------------------------------------------------------
# linear operator

# because Julia lack multiple inheritance, we have to define 4 of these
# make sure they share the same fields whenever needed

"""
  Dense linear operator for Newton's method.

  Subtype of [`AbstractDenseLO`](@ref)
"""
mutable struct NewtonDenseLO <: AbstractDenseLO
  lo_inner::DenseLO
  idata::ImplicitEulerData
  myrank::Int
  commsize::Int
end

"""
  Outer constructor for [`NewtonDenseLO`](@ref)
"""
function NewtonDenseLO(pc::PCNone, mesh::AbstractMesh,
                    sbp::AbstractOperator, eqn::AbstractSolutionData, opts::Dict)

  lo_inner = DenseLO(pc, mesh, sbp, eqn, opts)
  idata = ImplicitEulerData(mesh, opts)

  return NewtonDenseLO(lo_inner, idata, mesh.myrank, mesh.commsize)
end

"""
  Sparse direct linear operator for Newton's method.  Subtype of
  [`AbstractSparseDirectLO`](@ref)
"""
mutable struct NewtonSparseDirectLO <: AbstractSparseDirectLO
  lo_inner::SparseDirectLO
  idata::ImplicitEulerData
  myrank::Int
  commsize::Int
end

"""
  Outer constructor for [`NewtonSparseDirectLO`](@ref)
"""
function NewtonSparseDirectLO(pc::PCNone, mesh::AbstractMesh,
                    sbp::AbstractOperator, eqn::AbstractSolutionData, opts::Dict)

  lo_inner = SparseDirectLO(pc, mesh, sbp, eqn, opts)
  idata = ImplicitEulerData(mesh, opts)

  return NewtonSparseDirectLO(lo_inner, idata, mesh.myrank, mesh.commsize)
end

"""
  Petsc matrix based linear operator for Newton's method.

  Subtype of [`AbstractPetscMatLO`](@ref)
"""
mutable struct NewtonPetscMatLO <: AbstractPetscMatLO
  lo_inner::PetscMatLO
  idata::ImplicitEulerData
  myrank::Int
  commsize::Int
end

"""
  Outer constructor for [`NewtonPetscMatLO`](@ref)
"""
function NewtonPetscMatLO(pc::AbstractPetscPC, mesh::AbstractMesh,
                    sbp::AbstractOperator, eqn::AbstractSolutionData, opts::Dict)

  lo_inner = PetscMatLO(pc, mesh, sbp, eqn, opts)
  idata = ImplicitEulerData(mesh, opts)

  return NewtonPetscMatLO(lo_inner, idata, mesh.myrank, mesh.commsize)
end

"""
  Petsc matrix-free linear operator for Newton's method.

  Subtype of [`AbstractPetscMatFreeLO`](@ref)
"""
mutable struct NewtonPetscMatFreeLO <: AbstractPetscMatFreeLO
  lo_inner::PetscMatFreeLO
  idata::ImplicitEulerData
  myrank::Int
  commsize::Int
end

"""
  Newton mat-free linear operator constructor

  **Inputs**

   * pc
   * mesh
   * sbp
   * eqn
   * opts
   * rhs_func: rhs_func from [`newtonInner`](@ref)
"""
function NewtonPetscMatFreeLO(pc::AbstractPetscPC, mesh::AbstractMesh,
                    sbp::AbstractOperator, eqn::AbstractSolutionData, opts::Dict)

  lo_inner = PetscMatFreeLO(pc, mesh, sbp, eqn, opts)
  idata = ImplicitEulerData(mesh, opts)

  return NewtonPetscMatFreeLO(lo_inner, idata, mesh.myrank, mesh.commsize)
end

"""
  All Newton linear operators
"""
const NewtonLO = Union{NewtonDenseLO, NewtonSparseDirectLO, NewtonPetscMatLO, NewtonPetscMatFreeLO}

"""
  Newton matrix-explicit linear operators
"""
const NewtonMatLO = Union{NewtonDenseLO, NewtonSparseDirectLO, NewtonPetscMatLO}

"""
  Any PC or LO that has a matrix in the field `A`
"""
const NewtonHasMat = Union{NewtonMatPC, NewtonDenseLO, NewtonSparseDirectLO, NewtonPetscMatLO}

"""
  Any Newton PC or LO.
"""
const NewtonLinearObject = Union{NewtonDenseLO, NewtonSparseDirectLO, NewtonPetscMatLO, NewtonPetscMatFreeLO, NewtonMatPC}

"""
  Any Newton PC
"""
const NewtonPC = Union{NewtonMatPC, NewtonMatFreePC}

function calcLinearOperator(lo::NewtonMatLO, mesh::AbstractMesh,
                            sbp::AbstractOperator, eqn::AbstractSolutionData,
                            opts::Dict, ctx_residual, t)

  calcLinearOperator(lo.lo_inner, mesh, sbp, eqn, opts, ctx_residual, t)

  lo2 = getBaseLO(lo)
  physicsJac(mesh, sbp, eqn, opts, lo2.A, ctx_residual, t)

  if opts["newton_globalize_euler"]
    # TODO: updating the Euler parameter here is potentially wrong if we
    #       are not updating the Jacobian at every newton step
    updateEuler(lo)
    applyEuler(mesh, sbp, eqn, opts, lo)
  end

  return nothing
end

function calcLinearOperator(lo::NewtonPetscMatFreeLO, mesh::AbstractMesh,
                            sbp::AbstractOperator, eqn::AbstractSolutionData,
                            opts::Dict, ctx_residual, t)

  if opts["newton_globalize_euler"]
    updateEuler(lo)
  end

  setLOCtx(lo, mesh, sbp, eqn, opts, ctx_residual, t)

  return nothing
end

function applyLinearOperator(lo::NewtonPetscMatFreeLO, mesh::AbstractMesh,
                       sbp::AbstractOperator, eqn::AbstractSolutionData{Tsol},
                       opts::Dict, ctx_residual, t, x::AbstractVector,
                       b::AbstractVector) where Tsol

  @assert !(Tsol <: AbstractFloat)  # complex step only!

  epsilon =  opts["epsilon"]::Float64
  pert = Tsol(0, epsilon)

  # apply perturbation
  for i=1:mesh.numDof
    eqn.q_vec[i] += pert*x[i]
  end

  physicsRhs(mesh, sbp, eqn, opts, eqn.res_vec, ctx_residual, t)

  # calculate derivatives, store into b
  calcJacCol(b, eqn.res_vec, epsilon)

  if opts["newton_globalize_euler"]
    applyEuler(mesh, sbp, eqn, opts, x, lo, b)
  end

  # undo perturbation
  for i=1:mesh.numDof
    eqn.q_vec[i] -= pert*x[i]
  end

  return nothing
end

function applyLinearOperatorTranspose(lo::NewtonPetscMatFreeLO,
                             mesh::AbstractMesh,
                             sbp::AbstractOperator, eqn::AbstractSolutionData{Tsol},
                             opts::Dict, ctx_residual, t, x::AbstractVector,
                             b::AbstractVector) where Tsol

  error("applyLinearOperatorTranspose() not supported by NewtonPetscMatFreeLO")

end


text:
"""
    Maybe{X}

Type shortcut for "`X` or nothing" or "nullable `X`" in javaspeak. Name
got inspired by our functional friends.
"""
const Maybe{X} = Union{Nothing,X}

"""
Part of a measurement unit definition that corresponds to the SBML definition
of `Unit`. For example, the unit "per square megahour", Mh^(-2), is written as:

    UnitPart("second",  # base SI unit, this says we are measuring time
             -2,        # exponent, says "per square"
             6,         # log-10 scale of the unit, says "mega"
             1/3600)    # second-to-hour multiplier

Compound units (such as "volt-amperes" and "dozens of yards per ounce") are
built from multiple `UnitPart`s; see the definition of field `units` in
[`Model`](@ref).
"""
struct UnitPart
    kind::String
    exponent::Int
    scale::Int
    multiplier::Float64
    UnitPart(k, e, s, m) = new(k, e, s, m)
end

"""
Abstract type for all kinds of gene product associations
"""
abstract type GeneProductAssociation end

"""
Gene product reference in the association expression
"""
struct GPARef <: GeneProductAssociation
    gene_product::String
end

"""
Boolean binary "and" in the association expression
"""
struct GPAAnd <: GeneProductAssociation
    terms::Vector{GeneProductAssociation}
end

"""
Boolean binary "or" in the association expression
"""
struct GPAOr <: GeneProductAssociation
    terms::Vector{GeneProductAssociation}
end

"""
Reaction with stoichiometry that assigns reactants and products their relative
consumption/production rates (accessible in field `stoichiometry`), lower/upper
bounds (in tuples `lb` and `ub`, with unit names), and objective coefficient
(`oc`). Also may contains `notes` and `annotation`.
"""
struct Reaction
    stoichiometry::Dict{String,Float64}
    lb::Tuple{Float64,String}
    ub::Tuple{Float64,String}
    oc::Float64
    gene_product_association::Maybe{GeneProductAssociation}
    notes::Maybe{String}
    annotation::Maybe{String}
    Reaction(s, l, u, o, as, n = nothing, an = nothing) = new(s, l, u, o, as, n, an)
end

"""
Species metadata -- contains a human-readable `name`, a `compartment`
identifier, `formula`, `charge`, and additional `notes` and `annotation`.
"""
struct Species
    name::String
    compartment::String
    formula::Maybe{String}
    charge::Maybe{Int}
    notes::Maybe{String}
    annotation::Maybe{String}
    Species(na, co, f, ch, no = nothing, a = nothing) = new(na, co, f, ch, no, a)
end

"""
Gene product metadata.
"""
struct GeneProduct
    name::Maybe{String}
    label::Maybe{String}
    notes::Maybe{String}
    annotation::Maybe{String}
    GeneProduct(na, l, no = nothing, a = nothing) = new(na, l, no, a)
end

"""
Structure that collects the model-related data. Contains `parameters`, `units`,
`compartments`, `species` and `reactions` and `gene_products`, and additional
`notes` and `annotation` (also present internally in some of the data fields).
The contained dictionaries are indexed by identifiers of the corresponding
objects.
"""
struct Model
    parameters::Dict{String,Float64}
    units::Dict{String,Vector{UnitPart}}
    compartments::Vector{String}
    species::Dict{String,Species}
    reactions::Dict{String,Reaction}
    gene_products::Dict{String,GeneProduct}
    notes::Maybe{String}
    annotation::Maybe{String}
    Model(p, u, c, s, r, g, n = nothing, a = nothing) = new(p, u, c, s, r, g, n, a)
end


text:
#
# This is a script that re-generates the test cases used in experiment1.jl
# and experiment2.jl, that are save in MAT files. Having these tests already
# pre-computed saves an enourmous time while testing.
#

using MAT
using BivMatFun
using LinearAlgebra
using SpecialFunctions

include("experiments_common.jl")

println("Data directory: ", data_directory)

function generate_test_1()
    # This data needs to be kept in sync with experiment1.jl
    a = 1/2;
    f = (x,y,i,j) -> (x+y).^(-a-i-j) * gamma(a+i+j) ./ gamma(a) * (-1)^(i+j);
    ntests = 5;

    for i = 1 : ntests
        n = 32 * i;

        println("Generating experiment1, n = ", n, " random")

        A, B, C, Y = construct_random_benchmark(f, n)
        errest1 = BivMatFun.fun2m_condest(f, A, B, C) * eps();

        matwrite("$(data_directory)/experiment1_$(n)_random.mat", Dict(
        "A" => A,
        "B" => B,
        "C" => C,
        "Y" => Y,
        "errest" => errest1
        ));

        println("Generating experiment1, n = ", n, " grcar-random")

        A = grcar(n) + I

        A = convert(Matrix{ComplexF64}, A)
        B = convert(Matrix{ComplexF64}, B)
        C = convert(Matrix{ComplexF64}, C)

        Y2 = evaluate_reference_solution(f, A, B, C)
        errest2 = BivMatFun.fun2m_condest(f, A, B, C) * eps();

        matwrite("$(data_directory)/experiment1_$(n)_grcar.mat", Dict(
        "A" => A,
        "B" => B,
        "C" => C,
        "Y" => Y2,
        "errest" => errest2
        ));
    end
end

function generate_test_2()

    matrices = [ "jordbloc1", "grcar", "smoke", "kahan2", "lesp", "sampling", "grcar-randn" ];

    m = 64;
    n = m;
    a = .5;

    fncts = [
    Dict("name" => "sqrt", "f" => (x,y,i,j) -> sqrt(x+y)),
    Dict("name" => "invsqrt", "f" => (x,y,i,j)  -> 1 ./ sqrt(x+y)),
    Dict("name" => "phi1", "f" => (x,y,i,j) -> expm1(x+y) ./ (x+y)),
    Dict("name" => "expa", "f" => (x,y,i,j) -> exp((x+y).^a))
    ]

    shift = 1.00;

    lesp = n -> convert(Matrix{ComplexF64}, diagm(-1 => 1 ./ (2:n), 0 => -(5 : 2 : 2*n+3), 1 => 2 : n));
    smoke = n -> diagm(0 => exp.(2im * pi * (1 : n) ./ n), 1 => complex(ones(n-1)), -n+1 => [ 1.0 + 0im ]);

    for i = 1 : length(fncts)
        fname = fncts[i]["name"];
        f = fncts[i]["f"];

        for j = 1 : length(matrices);
            matrix_name = matrices[j];

            println("Generating experiment2, $(matrix_name), $(fname)")

            if matrix_name == "jordbloc1"
                k = min(m, n, 8);

                J = diagm(1 => ones(k-1)) + 0.1 * I;

                A = randn(m-k, m-k); A = A / opnorm(A); A = A + shift * I;
                A = cat(J, A, dims=[1,2]);
                Q = qr(randn(m, m)).Q;
                A = Q * A * Q';

                B = randn(n-k, n-k); B = B / opnorm(B); B = B + shift * I;
                B = cat(J, B, dims=[1,2]);
                Q = qr(randn(m, m)).Q;
                B = Q * B * Q';

            elseif matrix_name == "kahan2"
                k = min(min(m,n), 64);
                XA = randn(m-k, m-k)
                XA = XA / opnorm(XA) + shift * I;
                XB = randn(n-k, n-k)
                XB = XB / opnorm(XB) + shift * I;
                A = cat(kahan(k), XA, dims = [1,2]);
                B = cat(kahan(k), XB, dims = [1,2]);

            elseif matrix_name == "smoke"
                A = schur(smoke(m)).T;
                B = A;

            elseif matrix_name == "grcar-randn"
                A = grcar(m);
                B = randn(n, n);

            elseif matrix_name == "grcar"
                A = grcar(m);
                B = grcar(n);

            elseif matrix_name == "lesp"
                k = min(32, min(m,n));
                XA = randn(m-k, m-k)
                XA = XA / opnorm(XA) - shift * I;
                XB = randn(n-k, n-k)
                XB = XB / opnorm(XB) - shift * I;
                A = -cat(schur(lesp(k)).T, XA, dims = [1,2]);
                B = -cat(schur(lesp(k)).T, XB, dims = [1,2]);

            elseif matrix_name == "sampling"
                k = min(32, min(m,n));

                A = sampling(k) + I;
                B = sampling(k) + I;

                XA = randn(m-k, m-k)
                XA = XA / opnorm(XA) + shift * I;
                XB = randn(n-k, n-k)
                XB = XB / opnorm(XB) + shift * I;

                A = cat(A, XA, dims=[1,2])
                B = cat(B, XB, dims=[1,2])
            end

            A = complex(A);
            B = complex(B)

            C = randn(m, n); C = C / norm(C);
            C = complex(C)

            Y = evaluate_reference_solution(f, A, B, C)
            errest = BivMatFun.fun2m_condest(f, A, B, C) * eps();

            matwrite("$(data_directory)/experiment2_$(matrix_name)_$(fname).mat", Dict(
            "A" => A,
            "B" => B,
            "C" => C,
            "Y" => Y,
            "errest" => errest
            ))
        end
    end
end

generate_test_1()
generate_test_2()

max_stars_repo_name: numpi/BivMatFun

max_issues_repo_name: numpi/BivMatFun

max_forks_repo_name: numpi/BivMatFun


text:
function predict_nfolds(a_all_all,a_layer1_all,n_per_part_all,inds_all_all,x_new,i_new,vars,cc1,xx_fold_all,nof_folds,maetrs,neurons_all)
    predl1_all=Vector{Vector{Float64}}()
    for i=1:nof_folds
        println(Dates.format(now(), "HH:MM:SS"), " Predicting ",i, " of ",nof_folds)
        a_all_tmp=a_all_all[i]; a_layer1_tmp=a_layer1_all[i]; n_per_part_tmp=n_per_part_all[i]; inds_all_tmp=inds_all_all[i];xx_fold_tmp=xx_fold_all[i];
        predl1_new,layer_new=predict_new_rbf(a_all_tmp,a_layer1_tmp,x_new,i_new,neurons_all[i],n_per_part_tmp,inds_all_tmp,xx_fold_tmp,vars,cc1)
       push!(predl1_all,predl1_new) # .*(1 ./maetrs[i]))
    end
    return predl1_all
end


text:
using Test
using StringManipulation
using Markdown

@testset "Alignment" verbose = true begin
    include("./alignment.jl")
end

@testset "Cropping" verbose = true begin
    include("./crop.jl")
end

@testset "Decorations" verbose = true begin
    include("./decorations.jl")
end

@testset "Highlighting" verbose = true begin
    include("./highlighting.jl")
end

@testset "Search" verbose = true begin
    include("./search.jl")
end

@testset "Splitting" verbose = true begin
    include("./split.jl")
end

@testset "View" verbose = true begin
    include("./view.jl")
end

@testset "Width" verbose = true begin
    include("./width.jl")
end


text:
using Random
using Plots
include("../src/sample/direct.jl")
include("../src/vas.jl")

function sir_vas(β, c, γ)
    take = [
        1 0;
        1 1;
        0 0;
    ]
    give = [
        0 0;
        2 0;
        0 1;
    ]
    rates = [(state) -> Exponential(β*c*state[2]/sum(state)*state[1]), (state) -> Exponential(state[2] * γ)]
    (take, give, rates)
end

p = [0.05,10.0,0.25]; # β,c,γ
u0 = [990,10,0]; # S,I,R

take, give, rates = sir_vas(p...)

rng = MersenneTwister()

vas = VectorAdditionSystem(take, give, rates)
state = zero_state(vas)

vam = VectorAdditionModel(vas, state, 0.0)
fsm = VectorAdditionFSM(vam, DirectCall(Int))

initial_state = vas_initial(vas, u0)

out = Matrix{Float64}(undef, u0[2] + 2*u0[1] + 1, 4)

when, next_transition = simstep!(fsm, initial_state, rng)

event_cnt = 0
tnow = 0.0
out[event_cnt + 1, 2:4] = u0
out[event_cnt + 1, 1] = tnow

while next_transition !== nothing
    tnow += when
    out[event_cnt + 1, 1] = tnow
    out[event_cnt + 1, 2:4] = vam.state
    token = vas_delta(vas, next_transition)
    when, next_transition = simstep!(fsm, token, rng)
    event_cnt += 1
end

plot(
    out[1:event_cnt, 1],
    out[1:event_cnt, 2:4],
    label=["S" "I" "R"],
    xlabel="Time",
    ylabel="Number"
)


text:
module SparseInverseProblems
	module Util
	include("util/ip.jl")
	include("util/ip_lasso.jl")
	end
include("abstractTypes.jl")
include("lsLoss.jl")
include("BoxConstrainedDifferentiableModel.jl")
include("ADCG.jl")
end


text:
# tutorial on conventional 2D SLAM example

addprocs(3)

# This tutorial shows how to use some of the commonly used factor types
# This tutorial follows from the ContinuousScalar example from IncrementalInference

@everywhere begin
  using RoME, Distributions
  using IncrementalInference
  using TransformUtils
end # everywhere

# Start with an empty graph
fg = initfg(sessionname="SLAM2D_TUTORIAL")

# also add a PriorPose2 to pin the first pose at a fixed location
addNode!(fg, :x0, Pose2, labels=["POSE"])
addFactor!(fg, [:x0], Prior(MvNormal([0.0;0.0;0], diagm([1.0;1.0;0.01].^2))) )

# Drive around in a hexagon
for i in 0:5
  psym = Symbol("x$i")
  nsym = Symbol("x$(i+1)")
  addNode!(fg, nsym, Pose2, labels=["POSE"])
  addFactor!(fg, [psym;nsym], Pose2Pose2_NEW(MvNormal([10.0;0;pi/3], diagm([0.1;0.1;0.1].^2))))
end

# Graphs.plot(fg.g)

# Add landmarks with Bearing range measurements
addNode!(fg, :l1, Point2, labels=["LANDMARK"])
p2br = Pose2DPoint2DBearingRange(Normal(0,0.1),Normal(20.0,1.0))
addFactor!(fg, [:x0; :l1], p2br)

# Add landmarks with Bearing range measurements
p2br2 = Pose2DPoint2DBearingRange(Normal(0,0.1),Normal(20.0,1.0))
addFactor!(fg, [:x6; :l1], p2br2)

isInitialized(fg, :l1)

ensureAllInitialized!(fg)

using RoMEPlotting, Gadfly

pl = plotKDE(fg, [:x0; :x1; :x2; :x3; :x4; :x5; :x6]);

Gadfly.draw(PDF("tmpX0123456.pdf", 15cm, 15cm), pl)

@async run(`evince tmpX0123456.pdf`)

# pl = drawPoses(fg)
pl = drawPosesLandms(fg)
Gadfly.draw(PDF("tmpPosesFg.pdf", 15cm, 15cm), pl)
@async run(`evince tmpPosesFg.pdf`)

tree = wipeBuildNewTree!(fg)

@async Graphs.plot(tree.bt)

@time inferOverTree!(fg, tree)

# Graphs.plot(tree.bt)

@async Graphs.plot(fg.g)

# These functions need more work

using KernelDensityEstimate

function plot(::Type{Pose2}, p::Vector{BallTreeDensity})

  xy = plotKDE(p,dims=[1;2])
  th = plotKDE(p,dims=[3])

  hstack(xy, th)
end

vert0 = getVert(fg, :x0)
vert1 = getVert(fg, :x1)
pl = plot(typeof(getData(vert).softtype), [getKDE(vert0); getKDE(vert1)])

Gadfly.draw(PDF("tmpX01.pdf", 15cm, 15cm), pl)
@async run(`evince tmpX01.pdf`)

#


text:
using ComplexBigMatrices
using LinearAlgebra
using Test

@testset "ComplexBigMatrices.jl" begin
    # Create Float64 and BigFloat Matrices
    setprecision(512)
    a = [[1 1.0im];[1.0im 1]];
    biga = Complex{BigFloat}.(a)

    # @test sum(Float64.(abs.(eigen(a).vectors - eigen(biga).vectors))) = 0.0
    @test sum(Float64.(abs.(exp(a) - exp(biga)))) < 1e-12
    @test sum(Float64.(abs.(log(a) - log(biga)))) < 1e-12
end


text:
################################################################################
#
#   module LatPhysBandstructures
#   -> LatticePhysics_Base
#   -> LatticePhysics_ReciprocalSpace
#   -> LinearAlgebra
#
#   --> Definition of Bond Hamiltonians
#   --> Definition of Hamiltonians
#   --> Definition of Bandstructures
#
################################################################################

module LatPhysBandstructures

# Dependencies
using LatPhysBase
using LatPhysReciprocal
using LinearAlgebra
using Distributed
using Random
using HDF5

# BOND HAMILTONIANS

# Abstract type definition
include("bond_hamiltonian/abstract_bond_hamiltonian.jl")

# Concrete Hamiltonian: Sum of two Hamiltonians
include("bond_hamiltonian/concrete_bond_hamiltonian_sum.jl")

# Concrete Hamiltonian: Heisenberg
include("bond_hamiltonian/concrete_bond_spin_hamiltonian_heisenberg.jl")
# Concrete Hamiltonian: Kitaev
include("bond_hamiltonian/concrete_bond_spin_hamiltonian_kitaev.jl")

# Concrete Hamiltonian: Simple Hopping Hamiltonian (independent on label of bonds)
include("bond_hamiltonian/concrete_bond_hopping_hamiltonian_simple.jl")
# Concrete Hamiltonian: Simple Hopping Hamiltonian (NN)
include("bond_hamiltonian/concrete_bond_hopping_hamiltonian_simple_NN.jl")
# Concrete Hamiltonian: Simple Hopping Hamiltonian (NN)
include("bond_hamiltonian/concrete_bond_hopping_hamiltonian_dict.jl")

#= generate fancy syntax of passing the bond to the object directly
for H in Base.subtypes(AbstractBondHamiltonian)
    function (hb::H)(b :: AbstractBond)
        # return the bond term function
        return bondterm(hb, b)
    end
end =#

# HAMILTONIANS

# Abstract type definition
include("hamiltonian/abstract_hamiltonian.jl")

# Concrete type definition
include("hamiltonian/concrete_hamiltonian.jl")

# BANDSTRUCTURE

# Abstract type definition
include("bandstructure/abstract_bandstructure.jl")

# Concrete type definition
include("bandstructure/concrete_bandstructure.jl")

# CONSTANT ENERGY MANIFOLD

# Abstract type definition
include("energy_manifold/abstract_energy_manifold.jl")

# Concrete type definition
include("energy_manifold/concrete_energy_manifold.jl")

end # module


text:
# broadcasting

using Base.Broadcast: BroadcastStyle, Broadcasted

struct ROCArrayStyle{N} <: AbstractGPUArrayStyle{N} end
ROCArrayStyle(::Val{N}) where N = ROCArrayStyle{N}()
ROCArrayStyle{M}(::Val{N}) where {N,M} = ROCArrayStyle{N}()

BroadcastStyle(::Type{ROCArray{T,N}}) where {T,N} = ROCArrayStyle{N}()

Base.similar(bc::Broadcasted{ROCArrayStyle{N}}, ::Type{T}) where {N,T} =
    similar(ROCArray{T}, axes(bc))

Base.similar(bc::Broadcasted{ROCArrayStyle{N}}, ::Type{T}, dims) where {N,T} =
    ROCArray{T}(undef, dims)

## replace base functions with libdevice alternatives

rocfunc(f) = f
rocfunc(::Type{T}) where T = (x...) -> T(x...) # broadcasting type ctors isn't GPU compatible

Broadcast.broadcasted(::ROCArrayStyle{N}, f, args...) where {N} =
    Broadcasted{ROCArrayStyle{N}}(rocfunc(f), args, nothing)

const device_intrinsics = :[
    cos, cospi, sin, sinpi, tan, acos, asin, atan,
    cosh, sinh, tanh, acosh, asinh, atanh, angle,
    log, log10, log1p, log2, logb, ilogb,
    exp, exp2, exp10, expm1, ldexp,
    erf, erfinv, erfc, erfcinv, erfcx,
    brev, clz, ffs, byte_perm, popc,
    isfinite, isinf, isnan, nearbyint,
    nextafter, signbit, copysign, abs,
    sqrt, rsqrt, cbrt, rcbrt, pow,
    ceil, floor, saturate,
    lgamma, tgamma,
    j0, j1, jn, y0, y1, yn,
    normcdf, normcdfinv, hypot,
    fma, sad, dim, mul24, mul64hi, hadd, rhadd, scalbn].args

for f in device_intrinsics
    isdefined(Base, f) || continue
    @eval rocfunc(::typeof(Base.$f)) = $f
end

# broadcast ^

rocliteral_pow(::typeof(^), x::T, ::Val{0}) where {T<:Real} = one(x)
rocliteral_pow(::typeof(^), x::T, ::Val{1}) where {T<:Real} = x
rocliteral_pow(::typeof(^), x::T, ::Val{2}) where {T<:Real} = x * x
rocliteral_pow(::typeof(^), x::T, ::Val{3}) where {T<:Real} = x * x * x
rocliteral_pow(::typeof(^), x::T, ::Val{p}) where {T<:Real,p} = pow(x, Int32(p))

rocfunc(::typeof(Base.literal_pow)) = rocliteral_pow
rocfunc(::typeof(Base.:(^))) = pow

using MacroTools

const _rocfuncs = [copy(device_intrinsics); :^]
rocfuncs() = (global _rocfuncs; _rocfuncs)

_rocint(x::Int) = Int32(x)
_rocint(x::Expr) = x.head == :call && x.args[1] == :Int32 && x.args[2] isa Int ? Int32(x.args[2]) : x
_rocint(x) = x

function _rocpowliteral(x::Expr)
    if x.head == :call && x.args[1] == :(AMDGPU.rocfunc(^)) && x.args[3] isa Int32
        num = x.args[3]
        if 0 <= num <= 3
            sym = gensym(:x)
            new_x = Expr(:block, :($sym = $(x.args[2])))

            if iszero(num)
                push!(new_x.args, :(one($sym)))
            else
                unroll = Expr(:call, :*)
                for x = one(num):num
                    push!(unroll.args, sym)
                end
                push!(new_x.args, unroll)
            end

            x = new_x
        end
    end
    x
end
_rocpowliteral(x) = x

function replace_device(ex)
    global _rocfuncs
    MacroTools.postwalk(ex) do x
        x = x in _rocfuncs ? :(AMDGPU.rocfunc($x)) : x
        x = _rocint(x)
        x = _rocpowliteral(x)
        x
    end
end

macro rocfunc(ex)
    global _rocfuncs
    def = MacroTools.splitdef(ex)
    f = def[:name]
    def[:name] = Symbol(:cu, f)
    def[:body] = replace_device(def[:body])
    push!(_rocfuncs, f)
    quote
        $(esc(MacroTools.combinedef(def)))
        AMDGPU.rocfunc(::typeof($(esc(f)))) = $(esc(def[:name]))
    end
end


text:
# pointer intrinsics

export @typed_ccall

# TODO: can we use constant propagation instead of passing the alignment as a Val?

using Core: LLVMPtr

@generated function pointerref(ptr::LLVMPtr{T,A}, i::Int, ::Val{align}) where {T,A,align}
    sizeof(T) == 0 && return T.instance
    Context() do ctx
        eltyp = convert(LLVMType, T; ctx)

        T_int = convert(LLVMType, Int; ctx)
        T_ptr = convert(LLVMType, ptr; ctx)

        T_typed_ptr = LLVM.PointerType(eltyp, A)

        # create a function
        param_types = [T_ptr, T_int]
        llvm_f, _ = create_function(eltyp, param_types)

        # generate IR
        Builder(ctx) do builder
            entry = BasicBlock(llvm_f, "entry"; ctx)
            position!(builder, entry)

            typed_ptr = bitcast!(builder, parameters(llvm_f)[1], T_typed_ptr)
            typed_ptr = inbounds_gep!(builder, typed_ptr, [parameters(llvm_f)[2]])
            ld = load!(builder, typed_ptr)

            if A != 0
                metadata(ld)[LLVM.MD_tbaa] = tbaa_addrspace(A; ctx)
            end
            alignment!(ld, align)

            ret!(builder, ld)
        end

        call_function(llvm_f, T, Tuple{LLVMPtr{T,A}, Int}, :ptr, :(Int(i-one(i))))
    end
end

@generated function pointerset(ptr::LLVMPtr{T,A}, x::T, i::Int, ::Val{align}) where {T,A,align}
    sizeof(T) == 0 && return
    Context() do ctx
        eltyp = convert(LLVMType, T; ctx)

        T_int = convert(LLVMType, Int; ctx)
        T_ptr = convert(LLVMType, ptr; ctx)

        T_typed_ptr = LLVM.PointerType(eltyp, A)

        # create a function
        param_types = [T_ptr, eltyp, T_int]
        llvm_f, _ = create_function(LLVM.VoidType(ctx), param_types)

        # generate IR
        Builder(ctx) do builder
            entry = BasicBlock(llvm_f, "entry"; ctx)
            position!(builder, entry)

            typed_ptr = bitcast!(builder, parameters(llvm_f)[1], T_typed_ptr)
            typed_ptr = inbounds_gep!(builder, typed_ptr, [parameters(llvm_f)[3]])
            val = parameters(llvm_f)[2]
            st = store!(builder, val, typed_ptr)

            if A != 0
                metadata(st)[LLVM.MD_tbaa] = tbaa_addrspace(A; ctx)
            end
            alignment!(st, align)

            ret!(builder)
        end

        call_function(llvm_f, Cvoid, Tuple{LLVMPtr{T,A}, T, Int},
                      :ptr, :(convert(T,x)), :(Int(i-one(i))))
    end
end

Base.unsafe_load(ptr::Core.LLVMPtr, i::Integer=1, align::Val=Val(1)) =
    pointerref(ptr, Int(i), align)

Base.unsafe_store!(ptr::Core.LLVMPtr{T}, x, i::Integer=1, align::Val=Val(1)) where {T} =
    pointerset(ptr, convert(T, x), Int(i), align)

# pointer operations

# NOTE: this is type-pirating; move functionality upstream

LLVMPtr{T,A}(x::Union{Int,UInt,Ptr}) where {T,A} = reinterpret(LLVMPtr{T,A}, x)
LLVMPtr{T,A}() where {T,A} = LLVMPtr{T,A}(0)

# conversions from and to integers
Base.UInt(x::LLVMPtr) = reinterpret(UInt, x)
Base.Int(x::LLVMPtr) = reinterpret(Int, x)
Base.convert(::Type{LLVMPtr{T,A}}, x::Union{Int,UInt}) where {T,A} =
    reinterpret(LLVMPtr{T,A}, x)

Base.isequal(x::LLVMPtr, y::LLVMPtr) = (x === y)
Base.isless(x::LLVMPtr{T,A}, y::LLVMPtr{T,A}) where {T,A} = x < y

Base.:(==)(x::LLVMPtr{<:Any,A}, y::LLVMPtr{<:Any,A}) where {A} = UInt(x) == UInt(y)
Base.:(<)(x::LLVMPtr{<:Any,A},  y::LLVMPtr{<:Any,A}) where {A} = UInt(x) < UInt(y)
Base.:(==)(x::LLVMPtr, y::LLVMPtr) = false

Base.:(-)(x::LLVMPtr{<:Any,A},  y::LLVMPtr{<:Any,A}) where {A} = UInt(x) - UInt(y)

Base.:(+)(x::LLVMPtr, y::Integer) = oftype(x, Base.add_ptr(UInt(x), (y % UInt) % UInt))
Base.:(-)(x::LLVMPtr, y::Integer) = oftype(x, Base.sub_ptr(UInt(x), (y % UInt) % UInt))
Base.:(+)(x::Integer, y::LLVMPtr) = y + x

Base.unsigned(x::LLVMPtr) = UInt(x)
Base.signed(x::LLVMPtr) = Int(x)

# pointer type-preserving ccall

@generated function _typed_llvmcall(::Val{intr}, rettyp, argtt, args...) where {intr}
    # make types available for direct use in this generator
    rettyp = rettyp.parameters[1]
    argtt = argtt.parameters[1]
    argtyps = DataType[argtt.parameters...]
    argexprs = Expr[:(args[$i]) for i in 1:length(args)]

    # build IR that calls the intrinsic, casting types if necessary
    Context() do ctx
        T_ret = convert(LLVMType, rettyp; ctx)
        T_args = LLVMType[convert(LLVMType, typ; ctx) for typ in argtyps]

        llvm_f, _ = create_function(T_ret, T_args)
        mod = LLVM.parent(llvm_f)

        Builder(ctx) do builder
            entry = BasicBlock(llvm_f, "entry"; ctx)
            position!(builder, entry)

            # Julia's compiler strips pointers of their element type.
            # reconstruct those so that we can accurately look up intrinsics.
            T_actual_args = LLVMType[]
            actual_args = LLVM.Value[]
            for (arg, argtyp) in zip(parameters(llvm_f),argtyps)
                if argtyp <: LLVMPtr
                    # passed as i8*
                    T,AS = argtyp.parameters
                    actual_typ = LLVM.PointerType(convert(LLVMType, T; ctx), AS)
                    actual_arg = bitcast!(builder, arg, actual_typ)
                elseif argtyp <: Ptr
                    # passed as i64
                    T = eltype(argtyp)
                    actual_typ = LLVM.PointerType(convert(LLVMType, T; ctx))
                    actual_arg = inttoptr!(builder, arg, actual_typ)
                else
                    actual_typ = convert(LLVMType, argtyp; ctx)
                    actual_arg = arg
                end
                push!(T_actual_args, actual_typ)
                push!(actual_args, actual_arg)
            end

            # same for the return type
            if rettyp <: LLVMPtr
                T,AS = rettyp.parameters
                T_ret_actual = LLVM.PointerType(convert(LLVMType, T; ctx), AS)
            elseif rettyp <: Ptr
                T = eltype(rettyp)
                T_ret_actual = LLVM.PointerType(convert(LLVMType, T; ctx))
            else
                T_ret_actual = T_ret
            end

            intr_ft = LLVM.FunctionType(T_ret_actual, T_actual_args)
            intr_f = LLVM.Function(mod, String(intr), intr_ft)

            rv = call!(builder, intr_f, actual_args)

            if T_ret_actual == LLVM.VoidType(ctx)
                ret!(builder)
            else
                # also convert the return value
                if rettyp <: LLVMPtr
                    rv = bitcast!(builder, rv, T_ret)
                elseif rettyp <: Ptr
                    rv = ptrtoint!(builder, rv, T_ret)
                end

                ret!(builder, rv)
            end
        end

        call_function(llvm_f, rettyp, argtt, argexprs...)
    end
end

# perform a `ccall(intrinsic, llvmcall)` with accurate pointer types for calling intrinsic.
# this may be needed when selecting LLVM intrinsics, to avoid assertion failures, or when
# the back-end actually emits different code depending on types (e.g. SPIR-V and atomics).
#
# NOTE: this will become unnecessary when LLVM switches to typeless pointers,
#       or if and when Julia goes back to emitting exact types when passing pointers.
macro typed_ccall(intrinsic, cc, rettyp, argtyps, args...)
    # destructure and validate the arguments
    cc == :llvmcall || error("Can only use @typed_ccall with the llvmcall calling convention")
    Meta.isexpr(argtyps, :tuple) || error("@typed_ccall expects a tuple of argument types")

    # assign arguments to variables and unsafe_convert/cconvert them as per ccall behavior
    vars = Tuple(gensym() for arg in args)
    var_exprs = map(zip(vars, args)) do (var,arg)
        :($var = $arg)
    end
    arg_exprs = map(zip(vars,argtyps.args)) do (var,typ)
        :(Base.unsafe_convert($typ, Base.cconvert($typ, $var)))
    end

    esc(quote
        $(var_exprs...)
        GC.@preserve $(vars...) begin
            $_typed_llvmcall($(Val(Symbol(intrinsic))), $rettyp, Tuple{$(argtyps.args...)}, $(arg_exprs...))
        end
    end)
end


text:
@testset "cb" begin
    @test cb() == []

    f = @formula(y ~ cb(1))
    @test f.rhs.forig(f.rhs.args_parsed...) == [1]

    f = @formula(y ~ cb(1,2))
    @test f.rhs.forig(f.rhs.args_parsed...) == [1,2]

    f = @formula(y ~ cb(1,x))
    @test_throws ArgumentError f.rhs.forig(f.rhs.args_parsed...)
end

@testset "unpack" begin
    @test unpack(term(1)) == 1
    @test unpack(term(:x)) == :x
    @test unpack(term(:Unconditional)) == Unconditional()
    f = @formula(y ~ cb)
    @test unpack(f.rhs) == []
    f = @formula(y ~ cb(1,2))
    @test unpack(f.rhs) == [1,2]
end

@testset "exampledata" begin
    @test exampledata() == (:hrs, :nsw, :mpdta)
    hrs = exampledata("hrs")
    @test size(hrs) == (3280, 11)
    @test hrs[1] isa Vector
    @test size(exampledata(:nsw)) == (32834, 11)
    @test size(exampledata(:mpdta)) == (2500, 5)
end

@testset "checktable" begin
    @test_throws ArgumentError checktable(rand(3,3))
    @test_throws ArgumentError checktable([(a=1, b=2), (a=1, b=2)])
end


text:
# !((A & B) | (!A & !B))
#
# A B | A & B
# T T | T
# T F | F
# F T | F
# F F | F
#
# A B | !A & !B
# T T | F
# T F | F
# F T | F
# F F | T
#
# neural net layer 2[1] is going to be (A & B)
# neural net layer 2[2] is going to be (!A & !B)
#
# 2[1]: -7.5*bias + 5.0*A + 5.0*B
# 2[2]: 7.5*bias - 10.0*A - 10.0*B
#
# neural net layer 3[1] is going to be !(3[1] | 3[2])
# 3[1]: 5.0*bias + -10.0*A + -10.0*B
#
# A B | !(A | B)
# T T | F
# T F | F
# F T | F
# F F | T

#generate a multilayer perceptron
xornet = GenML.MLP.MultilayerPerceptron{Float64, (2,2,1)}()

#hand-written transition matrices.
xornet.layers[1].bias[:] = [-7.5, 7.5]
xornet.layers[1].transition[:] = [5.0 5.0; -10.0 -10.0]
xornet.layers[2].bias[:] = [5.0]
xornet.layers[2].transition[:] = [-10.0 -10.0]

@test xornet([true, true])[1] < 0.5
@test xornet([true, false])[1] > 0.5
@test xornet([false, true])[1] > 0.5
@test xornet([false, false])[1] < 0.5

#testing the allocations for evaluation.
dataset = rand(Bool, 2, 500)
results = Array{Float64}(1, 500)

GenML.evaluate!(results, xornet, dataset)
@time GenML.evaluate!(results, xornet, dataset)

max_stars_repo_name: interplanetary-robot/GenML

max_issues_repo_name: interplanetary-robot/GenML

max_forks_repo_name: interplanetary-robot/GenML


text:
import PyCall
import POMDPs: MDP, Policy

#
# In case you are planning to use the code for anything else than ToyText from
# OpenAI gym please pay attention to `initial_state`, `isterminal` and
# `action_index` functions as they are OpenAI specific. For example, OpenAI
# gym assume indices start from 0, but we change it to 1 etc.
#

#
# Define a new type responsible for working with ToyText from OpenAi gym
# It also supports other environments with a fixed number of states and actions
#
struct ToyTextMDP <: MDP{Int64, Int64}
    env::PyCall.PyObject
    state_space::Vector{Int64}
    action_space::Vector{Int64}
end

#
# Initializer function that is using PyObject with ToyText environment
# to initialize the variables
#
function ToyTextMDP(env::PyCall.PyObject)

    states = Vector{Int64}(1:env[:observation_space][:n])
    actions = Vector{Int64}(1:env[:action_space][:n])

    return ToyTextMDP(env, states, actions)
end

#
# We define number of function required to play the game
#
states(mdp::ToyTextMDP) = mdp.state_space
actions(mdp::ToyTextMDP) = mdp.action_space
n_states(mdp::ToyTextMDP) = length(mdp.state_space)
n_actions(mdp::ToyTextMDP) = length(mdp.action_space)
initial_state(mdp::ToyTextMDP, rng::AbstractRNG = MersenneTwister(0)) = mdp.env[:reset]()[] + 1
isterminal(mdp::ToyTextMDP, s::Int64) = s == 0
action_index(mdp::ToyTextMDP, a::Int64) = a - 1

#
# The only function we can use from OpenAIGym is getting a reward when taking a
# particular action in a specific step
#
function generate_sr(mdp::ToyTextMDP, s::Int64, a::Int64)

    state, reward, done = mdp.env[:step](a)

    if done
        return (zero(state), reward)
    end

    return (state + 1, reward)
end

# Play a single episode
function run_experiment(env::ToyTextMDP, policy::Policy, max_frame_iterations::Int64; keep_history = false)

    reward::Float64 = .0
    previous_state::Int64 = initial_state(env)

    #
    # We might decide to keep the history for some of our learning algorithms
    # therefore we need to initialize lists of a fixed size
    #
    list_size = if keep_history max_frame_iterations else zero(max_frame_iterations) end
    frames_played = 0

    states = zeros(Int64, list_size)
    new_states = zeros(Int64, list_size)
    actions = zeros(Int64, list_size)
    rewards = zeros(Float64, list_size)

    for i=1:max_frame_iterations

        current_action_i = action(policy, previous_state)
        current_action = action_index(env, current_action_i)
        current_state, current_reward = generate_sr(env, previous_state, current_action)
        reward += current_reward
        frames_played += 1

        #
        # In case we keep history lets track it
        #
        if keep_history
            states[i] = previous_state
            new_states[i] = current_state
            actions[i] = current_action_i
            rewards[i] = current_reward
        end

        if isterminal(env, current_state)
            break
        end

        previous_state = current_state

    end

    #
    # In case we keep history we should also cut in based on a number of frames required
    #
    if keep_history
        return states[1:frames_played], actions[1:frames_played], reward, new_states[1:frames_played], rewards[1:frames_played]
    else
        return states, actions, reward, new_states, current_reward
    end
end

#
# Executes policy on an environment for a number of repeats with a
# limitation on how many actions/ frames the episode can handle
#
function execute_policy(env::ToyTextMDP, policy::Policy, experiment_repeats::Int64, max_frame_iterations::Int64; experiment_index = nothing, verbose = false)

    if experiment_index != nothing && verbose == true
        println("Policy execution #$experiment_index")
    end

    return mean(map(x -> run_experiment(env, policy, max_frame_iterations)[3], zeros(experiment_repeats)))
end

max_stars_repo_name: dmitrijsc/practical-rl

max_issues_repo_name: dmitrijsc/practical-rl

max_forks_repo_name: dmitrijsc/practical-rl