text:
@inline calc_V(Y::Vector{<:Number}, p::AbstractModel)       = @inbounds Y[p.ind.Φ_s[1]] - Y[p.ind.Φ_s[end]]
@inline calc_I(Y::Vector{<:Number}, p::AbstractModel)       = @inbounds Y[p.ind.I[1]]
@inline calc_P(Y::Vector{<:Number}, p::AbstractModel)       = calc_I(Y,p)*p.θ[:I1C]*calc_V(Y,p)
@inline calc_c_e(Y::Vector{<:Number}, p::AbstractModel)     = @inbounds @views Y[p.ind.c_e]
@inline calc_c_s_avg(Y::Vector{<:Number}, p::AbstractModel) = @inbounds @views Y[p.ind.c_s_avg]
@inline calc_SOH(Y::Vector{<:Number}, p::AbstractModel)     = @inbounds Y[p.ind.SOH[1]]
@inline calc_j(Y::Vector{<:Number}, p::AbstractModel)       = @inbounds @views Y[p.ind.j]
@inline calc_Φ_e(Y::Vector{<:Number}, p::AbstractModel)     = @inbounds @views Y[p.ind.Φ_e]
@inline calc_Φ_s(Y::Vector{<:Number}, p::AbstractModel)     = @inbounds @views Y[p.ind.Φ_s]
@inline calc_film(Y::Vector{<:Number}, p::AbstractModel)    = @inbounds @views Y[p.ind.film]
@inline calc_j_s(Y::Vector{<:Number}, p::AbstractModel)     = @inbounds @views Y[p.ind.j_s]
@inline calc_Q(Y::Vector{<:Number}, p::AbstractModel)       = @inbounds @views Y[p.ind.Q]

@inline calc_T(Y::Vector{<:Number}, p::AbstractModelTemp{true}) = @inbounds @views Y[p.ind.T]
@inline calc_T(::Vector{<:Number}, p::AbstractModelTemp{false}) = repeat([p.θ[:T₀]], p.N.a+p.N.p+p.N.s+p.N.n+p.N.z)
@inline calc_K_eff(Y::Vector{<:Number}, p::AbstractModelTemp{true})  = @inbounds @views p.numerics.K_eff(Y[p.ind.c_e.p], Y[p.ind.c_e.s], Y[p.ind.c_e.n], Y[p.ind.T.p], Y[p.ind.T.s], Y[p.ind.T.n], p)
@inline calc_K_eff(Y::Vector{<:Number}, p::AbstractModelTemp{false}) = @inbounds @views p.numerics.K_eff(Y[p.ind.c_e.p], Y[p.ind.c_e.s], Y[p.ind.c_e.n], repeat([p.θ[:T₀]], p.N.p), repeat([p.θ[:T₀]], p.N.s), repeat([p.θ[:T₀]], p.N.n), p)

calc_η_plating(Y::Vector{<:Number},p::AbstractModel) = @views @inbounds Y[p.ind.Φ_s.n[1]] - Y[p.ind.Φ_e.n[1]]
calc_η_plating(t,Y,YP,p) = calc_η_plating(Y,p)

function calc_OCV(Y::AbstractVector{<:Number}, p::AbstractModel)
    """
    Calculate the open circuit voltages for the positive & negative electrodes
    """
    p_indices = c_s_indices(p, :p; surf=true, offset=true)
    n_indices = c_s_indices(p, :n; surf=true, offset=true)

    c_s_star_p = @views @inbounds Y[p_indices]
    c_s_star_n = @views @inbounds Y[n_indices]

    T = calc_T(Y,p)
    T_p = T[(1:p.N.p) .+ (p.N.a)]
    T_n = T[(1:p.N.n) .+ (p.N.a+p.N.p+p.N.s)]

    # Put the surface concentration into a fraction
    θ_p = c_s_star_p./p.θ[:c_max_p]
    θ_n = c_s_star_n./p.θ[:c_max_n]

    # Compute the OCV for the positive & negative electrodes.
    U_p = p.numerics.OCV_p(θ_p, T_p, p)[1]
    U_n = p.numerics.OCV_n(θ_n, T_n, p)[1]

    return return U_p, U_n
end

function calc_R_internal(Y::AbstractVector{<:Number}, p::AbstractModel)
    I = calc_I(Y,p)*calc_I1C(p)
    V = calc_V(Y,p)

    U_p, U_n = calc_OCV(Y,p)
    OCV = @inbounds U_p[1] - U_n[end]

    R_internal = abs((V - OCV)/I)

    return R_internal
end

# Metaprogramming to broadcast calc_x
for x in (:I,:V,:P,:c_e,:c_s_avg,:SOH,:j,:Φ_e,:Φ_s,:film,:j_s,:Q,:T,:K_eff,:η_plating,:OCV,:R_internal)
    name = "calc_$x"
    str = "Base.broadcasted(f::typeof($name),  Y::T, p::AbstractModel) where T<:VectorOfArray{Float64, 2, Vector{Vector{Float64}}} = [f(y,p) for y in Y]"

    str = remove_module_name(str)
    eval(Meta.parse(str))
end

@inline method_I(Y, p)   = calc_I(Y,p)
@inline method_V(Y, p)   = calc_V(Y,p)
@inline method_P(Y, p)   = calc_P(Y,p)

@inline scalar_residual!(res::Vector{T},t,Y,YP,p,run::run_constant{method,<:Any}) where {method<:AbstractMethod,T<:Number} = @inbounds (res[end] = method(Y,p) - value(run))

@inline scalar_residual!(res::Vector{T},t,Y,YP,p,run::run_function{method,func}) where {method<:AbstractMethod,T<:Num,func<:Function}     = @inbounds (res[end] = method(Y,p) - run.func(t,Y,YP,p))
@inline scalar_residual!(res::Vector{T},t,Y,YP,p,run::run_function{method,func}) where {method<:AbstractMethod,T<:Float64,func<:Function} = @inbounds (val = run.func(t,Y,YP,p); run.value .= val; res[end] = method(Y,p) - val)

@inline function scalar_jacobian(t,Y,YP,γ,p::AbstractModel,run::T) where T<:Union{run_constant,run_function}
    J = get_jacobian_sparsity(p,run;jac_type=typeof(t))
    J_nzval = @views @inbounds J.nzval[collect(1:length(J.nzval))]
    scalar_jacobian!(J_nzval,t,Y,YP,γ,p,run)
    return J
end
@inline @inbounds function scalar_jacobian!(J::SM,::T,Y::T2,::T2,γ::T,p::AbstractModel,::AbstractRun{method_I,<:Any}) where {T<:Number,SM<:SubArray{T,1,Vector{T},Tuple{Vector{Int64}},false},T2<:Vector{T}}
    J[1] = 1.0
    return nothing
end
@inline @inbounds function scalar_jacobian!(J::SM,::T,Y::T2,::T2,γ::T,p::AbstractModel,::AbstractRun{method_V,<:Any}) where {T<:Number,SM<:SubArray{T,1,Vector{T},Tuple{Vector{Int64}},false},T2<:Vector{T}}
    J[1] = 1.0
    J[2] = -1.0
    return nothing
end
@inline @inbounds function scalar_jacobian!(J::SM,::T,Y::T2,::T2,γ::T,p::AbstractModel,::AbstractRun{method_P,<:Any}) where {T<:Number,SM<:SubArray{T,1,Vector{T},Tuple{Vector{Int64}},false},T2<:Vector{T}}
    I1C = p.θ[:I1C]
    I = calc_I(Y,p)*I1C
    V = calc_V(Y,p)
    J[1] = I
    J[2] = -I
    J[3] = V*I1C
    return nothing
end

@inbounds function get_jacobian_sparsity(p::AbstractModel, ::AbstractRun{method_I,<:Any};jac_type::DataType=Float64)
    J = spzeros(jac_type,p.N.tot)
    J[p.ind.I[1]] = 1

    return J
end
@inbounds function get_jacobian_sparsity(p::AbstractModel, ::AbstractRun{method_V,<:Any};jac_type::DataType=Float64)
    J = spzeros(jac_type,p.N.tot)
    J[p.ind.Φ_s[[1,end]]] .= 1

    return J
end
@inbounds function get_jacobian_sparsity(p::AbstractModel, ::AbstractRun{method_P,<:Any};jac_type::DataType=Float64)
    J = spzeros(jac_type,p.N.tot)
    J[p.ind.I[1]] = 1
    J[p.ind.Φ_s[[1,end]]] .= 1

    return J
end

@inline function numargs(f::T,maxargs::Int64) where T<:Function
    maxval = 0
    maxargs += 1
    @inbounds for m in methods(f)
        val = num_types_in_tuple(m.sig)
        if val == maxargs
            maxval = val
            break
        elseif val ≤ maxargs && val > maxval
            maxval = val
        end
    end
    return maxval - 1
end
@inline num_types_in_tuple(sig::DataType) = length(sig.parameters)
@inline num_types_in_tuple(sig::UnionAll) = length(Base.unwrap_unionall(sig).parameters)

@inline function redefine_func(f::T;
    args::Int64=numargs(f,4),
    ) where T<:Function
    """
    Redefine the input function if its arguments aren't (t,Y,YP,p)
    """

    if     args == 4
        f_new = f
        # inputs = (:t,:Y,:YP,:p)
    elseif args == 3
        f_new = (t,Y,YP,p) -> f(t,Y,p)
        # inputs = (:t,:Y,:p)
    elseif args == 2
        f_new = (t,Y,YP,p) -> f(t,p)
        # inputs = (:t,:p)
    elseif args == 1
        f_new = (t,Y,YP,p) -> f(t)
        # inputs = (:t)
    else
        error("Input function must have at least one argument.")
    end

    return f_new#, inputs
end

function _get_method_funcs(p::model, run::run_function)
    θ_sym, Y, YP, t, SOC, I, γ, p_sym, θ_keys = get_symbolic_vars(p; original_keys=p.cache.θ_keys)
    res = similar(Y) .= 0.0

    scalar_residual!(res,t,Y,YP,p_sym,run)
    Jac_constant = scalar_jacobian(t,Y,YP,γ,p_sym,run)

    J_Y = J_YP = J_vec = nothing
    is_differentiable = true
    try
        J_Y  = @inbounds sparsejacobian([res[end]], Y)[:]
        J_YP = @inbounds sparsejacobian([res[end]], YP)[:]
        J_vec = J_Y .+ γ.*J_YP
    catch
        is_differentiable = false
    end

    if is_differentiable && !(Jac_constant == J_vec)
        return differentiate_residual_func(p,run,J_vec,J_Y,J_YP,res,θ_sym,Y,YP,t,SOC,I,γ,p_sym,θ_keys)
    else
        return _get_method_funcs_no_differentiation(p,run)
    end
end

function differentiate_residual_func(p::model,run::T,J_vec,J_Y,J_YP,res,θ_sym,Y,YP,t,SOC,I,γ,p_sym,θ_keys) where T<:run_function
    scalar_contains_Y_diff = @inbounds !isempty(J_Y[1:p.N.diff].nzval)
    scalar_contains_YP     = !isempty(J_YP.nzval)

    J_sp_scalar = spzeros(Float64,p.N.tot)
    @inbounds J_sp_scalar[J_vec.nzind] .= 1

    """
    Updating the theta vector to ensure any new parameters will be accounted for
    """
    θ_keys_scalar = @inbounds get_only_θ_used_in_sol(θ_sym, θ_keys, res[end])[2]

    θ_keys = deepcopy(p.cache.θ_keys)
    @inbounds for key in θ_keys_scalar
        if !(key ∈ θ_keys)
            push!(θ_keys, key)
        end
    end

    θ_tot      = [p.θ[key] for key in θ_keys]
    θ_sym_slim = [p_sym.θ[key] for key in θ_keys]
    update_θ!(θ_tot,θ_keys,p)

    J_scalar_func = eval(build_function(J_vec.nzval,t,Y,YP,γ,θ_sym_slim; expression=Val{false})[2])

    J_sp_base = p.funcs.J_y!.sp
    J_base_func = p.funcs.J_y!

    """
    Creating the scalar function
    """
    if scalar_contains_YP
        vars_in_residual = @inbounds get_variables(res[end])
        # Find all indices where YP is used in the residual equation
        ind_differential = findall(in(vars_in_residual), @views @inbounds YP[1:p.N.diff])

        @assert !isempty(ind_differential)

        res_algebraic = @inbounds res[end]

        try
            if !options[:SAVE_SYMBOLIC_FUNCTIONS] error() end
            f_new! = truncated_res_diff(p,ind_differential)
            f_new!(res,t,Y,YP,θ_sym_slim)
        catch
            residuals_PET!(res,t,Y,YP,p_sym)
        end
        @inbounds for ind in ind_differential
            res_algebraic = substitute(res_algebraic, Dict(YP[ind] => res[ind] + YP[ind]))
        end

        scalar_residal_alg! = eval(build_function(res_algebraic,t,Y,YP,θ_sym_slim; expression=Val{false}))

        J_alg_vec = @inbounds sparsejacobian([res_algebraic], Y[p.N.diff+1:end])[:]

        J_scalar_alg_func = eval(build_function(J_alg_vec.nzval,t,Y,YP,γ,θ_sym_slim; expression=Val{false})[2])

        J_sp_alg_scalar = spzeros(Float64,p.N.alg)
        @inbounds J_sp_alg_scalar[J_alg_vec.nzind] .= 1
    else
        scalar_residal_alg! = scalar_residual!

        if scalar_contains_Y_diff
            J_scalar_alg_func = eval(build_function((@inbounds J_vec[p.N.diff+1:end].nzval),t,Y,YP,γ,θ_sym_slim; expression=Val{false})[2])
        else
            J_scalar_alg_func = J_scalar_func
        end
        J_sp_alg_scalar = @inbounds J_sp_scalar[p.N.diff+1:end]
    end

    return combine_Jac_and_res(p,J_sp_base,J_base_func,J_sp_scalar,J_scalar_func,θ_tot,θ_keys,scalar_residual!,scalar_residal_alg!,J_scalar_alg_func,J_sp_alg_scalar)
end

function truncated_res_diff(p::AbstractModel,ind_differential)
    path = strings_directory_func(p; create_dir=false) * "/f_diff.jl"
    str_old = join(readlines(path), "\n");

    # If the function doesn't have rearranged `if` statements
    if !contains(str_old, ";")
        str_old = rearrange_if_statements(str_old)
    end

    str_old = PETLION.convert_to_ifelse(str_old)

    str = collect(str_old)

    # Remove all the unnecessary lines
    out_char(i::Int) = collect("ˍ₋out[$i]")
    ind = PETLION.find_next(str, 1, out_char(1))

    @inbounds for i in 1:p.N.diff
        ind = PETLION.find_next(str, ind[1], out_char(i))
        ind_last = PETLION.find_next(str, ind[end], '\n')
        if i ∈ ind_differential
            # keep this line and skip over it
            ind = ind_last[end]
        else
            deleteat!(str, ind[1]:ind_last[end])
        end
    end

    f_new = eval(Meta.parse(join(str)))

    return f_new
end

function _get_method_funcs_no_differentiation(p::model, run::AbstractRun)
    J_sp_scalar = get_jacobian_sparsity(p,run)
    θ_tot = p.cache.θ_tot
    θ_keys = p.cache.θ_keys
    update_θ!(θ_tot,θ_keys,p)

    J_scalar_func = scalar_jacobian!
    J_sp_base = p.funcs.J_y!.sp
    J_base_func = p.funcs.J_y!

    scalar_residal_alg! = scalar_residual!
    J_scalar_alg_func = J_scalar_func
    J_sp_alg_scalar = J_sp_scalar[p.N.diff+1:end]

    return combine_Jac_and_res(p,J_sp_base,J_base_func,J_sp_scalar,J_scalar_func,θ_tot,θ_keys,scalar_residual!,scalar_residal_alg!,J_scalar_alg_func,J_sp_alg_scalar)
end

function combine_Jac_and_res(p,J_sp_base,J_base_func,J_sp_scalar,J_scalar_func,θ_tot,θ_keys,scalar_residual!,scalar_residal_alg!,J_scalar_alg_func,J_sp_alg_scalar)
    J_full = _get_jacobian_combined(
        J_sp_base,
        J_base_func,
        J_sp_scalar,
        J_scalar_func,
        θ_tot,
        θ_keys;
        lu_decomposition=true,
        )

    R_full = residual_combined(
        p.funcs.f_diff!,
        p.funcs.f_alg!,
        scalar_residual!,
        1:p.N.diff,
        p.N.diff+1:p.N.tot,
        θ_tot,
        θ_keys,
    )

    J_alg = _get_jacobian_combined(
        (@inbounds J_sp_base[p.N.diff+1:p.N.tot-1,p.N.diff+1:p.N.tot]),
        p.funcs.J_y_alg!,
        J_sp_alg_scalar,
        J_scalar_alg_func,
        θ_tot,
        θ_keys;
        lu_decomposition=true,
        )

    R_diff = residual_combined(
        p.funcs.f_diff!,
        emptyfunc,
        emptyfunc,
        1:p.N.diff,
        1:0,
        θ_tot,
        θ_keys,
    )

    R_alg = residual_combined(
        emptyfunc,
        p.funcs.f_alg!,
        scalar_residal_alg!,
        1:0,
        1:p.N.alg,
        θ_tot,
        θ_keys,
    )

    return Jac_and_res(J_full,R_full,J_alg,R_diff,R_alg,Sundials.IDAIntegrator[])
end

function factorize(x...;kw...)
    method =  Symbol(lowercase(String(options[:FACTORIZATION_METHOD])))
    if method == :lu
        factor = LinearAlgebra.lu(x...;kw...)
    elseif method == :klu
        factor = klu(x...;kw...)
    else
        error("FACTORIZATION_METHOD $(options[:FACTORIZATION_METHOD]) is not supported.")
    end
    return factor
end
function _get_jacobian_combined(J_sp_base,J_base_func,J_sp_scalar,J_scalar_func,θ_tot,θ_keys;lu_decomposition=false)
    J_sp = [J_sp_base; J_sp_scalar']

    if lu_decomposition
        J_sp.nzval .= rand(length(J_sp.nzval))
        factor = factorize(J_sp)
    else
        factor = factorize(sparse([1],[1],[1.0]))
    end

    ind_base   = findall(J_sp.rowval .< size(J_sp,1))
    ind_scalar = findall(J_sp.rowval .== size(J_sp,1))
    if ind_base == 1:length(ind_base)
        ind_base = 1:length(ind_base)
    end

    J_base   = @views @inbounds J_sp.nzval[ind_base]
    J_scalar = @views @inbounds J_sp.nzval[ind_scalar]

    J = jacobian_combined(J_sp,J_base_func,J_base,J_scalar_func,J_scalar,θ_tot,θ_keys,factor)

    return J
end

function get_method_funcs!(p::model,run::run_constant{method,<:Any}) where method<:AbstractMethod
    p.funcs.Dict_constant[method] = _get_method_funcs_no_differentiation(p,run)
end
function get_method_funcs!(p::model,run::run_function{method,func}) where {method<:AbstractMethod,func<:Function}
    if !haskey(p.funcs.Dict_function,method)
        p.funcs.Dict_function[method] = Dict{DataType,jacobian_combined}()
    end

    func_only_uses_time = true
    try
        # Seeing if the residual function uses Y or YP. If there's an error or if it does not output a number,
        # get a Jacobian which is differentiable. Otherwise, the Jacobian is known
        if !(run.func(0.0,nothing,nothing,p) isa Number)
            func_only_uses_time = false
        end
    catch
        func_only_uses_time = false
    end

    if func_only_uses_time
        p.funcs.Dict_function[method][func] = _get_method_funcs_no_differentiation(p,run)
    else
        p.funcs.Dict_function[method][func] = _get_method_funcs(p,run)
    end

    return nothing
end

"""
Multiple dispatch for the residuals function
"""
@inline function (r::residual_combined{T1,T2,T3})(res::T,t,Y,YP,p,run) where {T1<:Function,T2<:Function,T3<:Function,T<:AbstractVector{Float64}}
    r.f_diff!((@views @inbounds res[r.ind_diff]), t, Y, YP, r.θ_tot)
    r.f_alg!( (@views @inbounds res[r.ind_alg]),  t, Y, YP, r.θ_tot)
    @inbounds res[end] = r.f_scalar!(t,Y,YP,r.θ_tot)
    return nothing
end
@inline function (r::residual_combined{T1,T2,T3})(res::T,t,Y,YP,p,run) where {T1<:typeof(emptyfunc),T2<:Function,T3<:Function,T<:AbstractVector{Float64}}
    r.f_alg!(res,t,Y,YP,r.θ_tot)
    @inbounds res[end] = r.f_scalar!(t,Y,YP,r.θ_tot)
    return nothing
end
@inline function (r::residual_combined{T1,T2,T3})(res::T,t,Y,YP,p,run) where {T1<:Function,T2<:Function,T3<:typeof(scalar_residual!),T<:AbstractVector{Float64}}
    r.f_diff!((@views @inbounds res[r.ind_diff]), t, Y, YP, r.θ_tot)
    r.f_alg!( (@views @inbounds res[r.ind_alg]),  t, Y, YP, r.θ_tot)
    scalar_residual!(res,t,Y,YP,p,run)
    return nothing
end
@inline function (r::residual_combined{T1,T2,T3})(res::T,t,Y,YP,p,run) where {T1<:typeof(emptyfunc),T2<:Function,T3<:typeof(scalar_residual!),T<:AbstractVector{Float64}}
    r.f_alg!(res,t,Y,YP,r.θ_tot)
    scalar_residual!(res,t,Y,YP,p,run)
    return nothing
end
@inline function (r::residual_combined{T1,T2,T3})(res::T,t,Y,YP,p,run) where {T1<:Function,T2<:typeof(emptyfunc),T3<:typeof(emptyfunc),T<:AbstractVector{Float64}}
    r.f_diff!(res, t, Y, YP, r.θ_tot)
    return nothing
end

"""
Multiple dispatch for the Jacobian function
"""
@inline function (J::jacobian_combined)(J_new::SparseMatrixCSC{Float64,Int64},x...)
    J(x...)
    @inbounds J_new.nzval .= J.sp.nzval
    return nothing
end
@inline function (J::jacobian_combined{T1,T2,T3})(t,Y,YP,γ,p,run) where {T1<:Function,T2,T3<:Function}
    J.base_func(J.J_base,t,Y,YP,γ,J.θ_tot)
    J.scalar_func(J.J_scalar,t,Y,YP,γ,J.θ_tot)
    return nothing
end
@inline function (J::jacobian_combined{T1,T2,T3})(t,Y,YP,γ,p,run) where {T1<:Function,T2,T3<:typeof(scalar_jacobian!)}
    J.base_func(J.J_base,t,Y,YP,γ,J.θ_tot)
    scalar_jacobian!(J.J_scalar,t,Y,YP,γ,p,run)
    return nothing
end
@inline function (J::jacobian_combined{T1,T2,T3})(t,Y::Vector{Float64},YP::AbstractVector{Float64},γ,p::model{<:jacobian_AD},run) where {T1<:Function,T2,T3<:Function}
    J.scalar_func(J.J_scalar,t,Y,YP,γ,J.θ_tot)
    res_FD = J.base_func.f!
    if size(J.sp) == (p.N.alg,p.N.alg)
        @inbounds @views res_FD.Y_cache[1:res_FD.N.diff] .= Y[1:res_FD.N.diff]
        Y_new = @views @inbounds Y[p.N.diff+1:end]
        @inbounds res_FD.YP_cache .= 0.0
    else
        @inbounds res_FD.YP_cache .= YP
        Y_new = Y
    end
    J.base_func(t,Y_new,YP,γ,p,run)
    J.J_base .= J.base_func.sp.nzval
    return nothing
end
@inline function (J::jacobian_combined{T1,T2,T3})(t,Y::Vector{Float64},YP::AbstractVector{Float64},γ,p::model{<:jacobian_AD},run) where {T1<:Function,T2,T3<:typeof(scalar_jacobian!)}
    scalar_jacobian!(J.J_scalar,t,Y,YP,γ,p,run)
    res_FD = J.base_func.f!
    if size(J.sp) == (p.N.alg,p.N.alg)
        @inbounds @views res_FD.Y_cache[1:res_FD.N.diff] .= Y[1:res_FD.N.diff]
        Y_new = @views @inbounds Y[p.N.diff+1:end]
        @inbounds res_FD.YP_cache .= 0.0
    else
        @inbounds res_FD.YP_cache .= YP
        Y_new = Y
    end
    J.base_func(t,Y_new,YP,γ,p,run)
    J.J_base .= J.base_func.sp.nzval
    return nothing
end

"""
model_funcs definitions
"""
(f::model_funcs)(::run_constant{method,input}) where {method<:AbstractMethod,input<:Any}     = f.Dict_constant[method]
(f::model_funcs)(::run_function{method,func})  where {method<:AbstractMethod,func<:Function} = f.Dict_function[method][func]

Base.haskey(f::model_funcs,::run_constant{method,input})  where {method<:AbstractMethod,input<:Any}     = haskey(f.Dict_constant,method)
Base.haskey(f::model_funcs,::run_function{method,func})   where {method<:AbstractMethod,func<:Function} = haskey(f.Dict_function,method) && haskey(f.Dict_function[method],func)


text:
export pred_ared_computation
function pred_ared_computation(gₖ :: Float64,
                               fₖ :: Float64,
                               dersec :: Float64,
                               d :: Float64,
                               ftestTR :: Float64,
                               gtestTR :: Float64;
                               seuil :: Float64=-1e-10)

  pred = gₖ * d + 0.5 * dersec * d^2

  if pred > seuil
    #print_with_color(:yellow,"!!")
    ared = (gₖ + gtestTR) * d / 2
  else
    #print_with_color(:green,"!!")
    ared = ftestTR - fₖ
  end

  ratio = ared / pred

  return (pred, ared, ratio)
end

ext: jl


text:
using MarketData

facts("time-related utilities") do

  context("parsedatetime") do
      @fact 1 => roughly(1.0000000001)
  end

  context("make datetime from array of Int") do
      @fact 1 => roughly(1.0000000001)
  end

  context("convert date to datetime with the last second of the day") do
    @fact 1.0000000000000001 => roughly(1)
  end
end

facts("signal utilities") do

  context("discretesignal") do
    @fact_throws log(-1)
    @fact_throws sum(foo, bar)
  end
end


text:
export
    create_database,
    create_table,
    create_all,

    from_sql,
    apply!

# ######
# # db #
# ######

const QSDB = string(QuantumStatesData)

function create_database(; dbname::String=QSDB)
    create_database(dbname)
    enable_uuid(dbname)
end

function create_table(table::Type{<:QuantumStatesData}; dbname::String=QSDB)
    dbconfig = copy(current_dbconfig())
    dbconfig[:dbname] = dbname

    sql = gen_table_schema(table)

    return create_table(string(table), sql, dbconfig=dbconfig)
end

function create_all(; dbname::String=QSDB)
    create_database(dbname=dbname)

    for qs in subtypes(QuantumStatesData)
        create_table(qs, dbname=dbname)
    end
end

# ##############
# # preprocess #
# ##############

# ###############
# # postprocess #
# ###############

function nrow(table::Type{<:QuantumStatesData}; dbname::String=QSDB)
    dbconfig = copy(current_dbconfig())
    dbconfig[:dbname] = dbname

    df = from_sql("select count(*) from $(string(table));", dbconfig=dbconfig)

    return Int(df[1, 1])
end

function from_sql(
    table::Type{<:QuantumStatesData},
    n::Integer;
    offset=-1,
    order=:id,
    dbname::String=QSDB
)
    dbconfig = copy(current_dbconfig())
    dbconfig[:dbname] = dbname

    (offset < 0) && (offset = rand(0:nrow(table)-n))

    sql = """
        SELECT *
            FROM $(string(table))
            ORDER BY $order
        LIMIT $n OFFSET $offset;
    """

    return from_sql(sql, dbconfig=dbconfig)
end

function apply!(df, pairs...)
    for (col, f) in pairs
        df[!, col] = f.(df[!, col])
    end

    return df
end


text:
using Nefarious
using Test

@testset "Nefarious.jl" begin
    # Write your own tests here.
end

max_stars_repo_name: StefanKarpinski/Nefarious

max_issues_repo_name: StefanKarpinski/Nefarious

max_forks_repo_name: StefanKarpinski/Nefarious


text:
#=
Purpose: time utilities that depend on custom Types go here

Difference from CoreUtils/time.jl functions here require SeisIO Types to work
=#
function mk_t!(C::GphysChannel, nx::Integer, ts_new::Int64)
  T = Array{Int64, 2}(undef, 2, 2)
  setindex!(T, one(Int64), 1)
  setindex!(T, nx, 2)
  setindex!(T, ts_new, 3)
  setindex!(T, zero(Int64), 4)
  setfield!(C, :t, T)
  return nothing
end

function check_for_gap!(S::GphysData, i::Integer, ts_new::Int64, nx::Integer, v::Integer)
  Δ = round(Int64, sμ / getindex(getfield(S, :fs), i))
  T = getindex(getfield(S, :t), i)
  T1 = t_extend(T, ts_new, nx, Δ)
  if T1 != nothing
    if v > 1
      te_old = endtime(T, Δ)
      δt = ts_new - te_old
      (v > 1) && println(stdout, lpad(S.id[i], 15), ": time difference = ", lpad(δt, 16), " μs (old end = ", lpad(te_old, 16), ", new start = ", lpad(ts_new, 16), ", gap = ", lpad(δt-Δ, 16), " μs)")
    end
    S.t[i] = T1
  end
  return nothing
end


text:
using Utils
using Underscores

nums(s) = @_ split(s) |> map(parse(Int, _), __)

struct Rule
    refs::Array{Int}
end

struct Alternative
    left::Array{Int}
    right::Array{Int}
end

const RulePossibilities = Union{AbstractString,Rule,Alternative}
const Rules = Dict{Int,RulePossibilities}

struct Messages
    rules::Rules
    messages::Array{SubString}
end

function parse_file(f)
    rules, messages = split_blocks(f)

    rls = Rules()
    for line in eachline(IOBuffer(rules))
        head, content = split(line, ": ")
        head = parse(Int, head)
        if startswith(content, '"')
            rls[head] = strip(content, '"')
        else
            if occursin('|', content)
                a, b = split(content, " | ")
                rls[head] = Alternative(nums(a), nums(b))
            else
                rls[head] = Rule(nums(content))
            end
        end
    end
    Messages(rls, collect(eachline(IOBuffer(messages))))
end

function load()
    open(joinpath(@__DIR__, "19.txt"), "r") do f
        parse_file(f)
    end
end

messages = load()

function build_regex(rules::Rules, current)
    current = rules[current]
    if isa(current, Alternative)
        lefts = @_ current.left |> map(build_regex(rules, _), __) |> join
        rights = @_ current.right |> map(build_regex(rules, _), __) |> join
        "(?:$lefts|$rights)"
    elseif isa(current, Rule)
        @_ current.refs |> map(build_regex(rules, _), __) |> join
    else
        current
    end
end

build_regex(rules) = '^' * build_regex(rules, 0) * '$' |> Regex

function ex1(messages)
    re = build_regex(messages.rules)
    @_ count(occursin(re, _), messages.messages)
end

function rule_walker(_, str, current::AbstractString, pos)::Array{Int,1}
    [startswith(str[pos:end], current) ? pos + length(current) : 0]
end

function rule_walker(rules::Rules, str, current::Int, pos)::Array{Int,1}
    rule_walker(rules, str, rules[current], pos)
end

function rule_walker(rules::Rules, str, current::Alternative, pos)::Array{Int,1}
    re = rule_walker(rules, str, Rule(current.left), pos)
    append!(re, rule_walker(rules, str, Rule(current.right), pos))
    filter(p -> p != 0, re)
end

function rule_walker(rules::Rules, str, current::Rule, pos)::Array{Int,1}
    positions = [pos]
    for (i,ref) in enumerate(current.refs)
        positions = @_ map(rule_walker(rules, str, ref, _), positions) |> Iterators.flatten |> collect
        if i < lastindex(current.refs)
            filter!(p -> p in eachindex(str), positions)
        end
    end
    positions
end

rule_walker(rules, str) = (length(str) + 1) in rule_walker(rules, str, 0, 1)

function ex2(messages, update=true)
    changed = copy(messages.rules)
    if update
        changed[8] = Alternative([42], [42, 8])
        changed[11] = Alternative([42, 31], [42, 11, 31])
    end

    @_ count(rule_walker(changed, _), messages.messages)
end

println("Number of matching messages: ", ex1(messages))
println("Matching after change: ", ex2(messages))

using Test

@testset "Monster Messages" begin
    input1 = """
    0: 4 1 5
    1: 2 3 | 3 2
    2: 4 4 | 5 5
    3: 4 5 | 5 4
    4: "a"
    5: "b"

    ababbb
    bababa
    abbbab
    aaabbb
    aaaabbb
    """

    example1 = parse_file(IOBuffer(input1))

    input2 = """
    42: 9 14 | 10 1
    9: 14 27 | 1 26
    10: 23 14 | 28 1
    1: "a"
    11: 42 31
    5: 1 14 | 15 1
    19: 14 1 | 14 14
    12: 24 14 | 19 1
    16: 15 1 | 14 14
    31: 14 17 | 1 13
    6: 14 14 | 1 14
    2: 1 24 | 14 4
    0: 8 11
    13: 14 3 | 1 12
    15: 1 | 14
    17: 14 2 | 1 7
    23: 25 1 | 22 14
    28: 16 1
    4: 1 1
    20: 14 14 | 1 15
    3: 5 14 | 16 1
    27: 1 6 | 14 18
    14: "b"
    21: 14 1 | 1 14
    25: 1 1 | 1 14
    22: 14 14
    8: 42
    26: 14 22 | 1 20
    18: 15 15
    7: 14 5 | 1 21
    24: 14 1

    abbbbbabbbaaaababbaabbbbabababbbabbbbbbabaaaa
    bbabbbbaabaabba
    babbbbaabbbbbabbbbbbaabaaabaaa
    aaabbbbbbaaaabaababaabababbabaaabbababababaaa
    bbbbbbbaaaabbbbaaabbabaaa
    bbbababbbbaaaaaaaabbababaaababaabab
    ababaaaaaabaaab
    ababaaaaabbbaba
    baabbaaaabbaaaababbaababb
    abbbbabbbbaaaababbbbbbaaaababb
    aaaaabbaabaaaaababaa
    aaaabbaaaabbaaa
    aaaabbaabbaaaaaaabbbabbbaaabbaabaaa
    babaaabbbaaabaababbaabababaaab
    aabbbbbaabbbaaaaaabbbbbababaaaaabbaaabba
    """

    example2 = parse_file(IOBuffer(input2))

    @testset "interna" begin
        @testset "string rule" begin
            @test rule_walker(Rules(), "abcd", "ab", 1) == [3]
            @test rule_walker(Rules(), "abcd", "abc", 1) == [4]
            @test rule_walker(Rules(), "abcd", "d", 4) == [5]
            @test rule_walker(Rules(), "abcd", "bc", 2) == [4]
            @test rule_walker(Rules(), "abcd", "bc", 1) == [0]
        end

        @testset "int rule" begin
            @assert rule_walker(Rules(0 => "ab"), "abcd", 0, 1) == [3]
            @assert rule_walker(Rules(1 => "abc"), "abcd", 1, 1) == [4]
            @assert rule_walker(Rules(2 => "d"), "abcd", 2, 4) == [5]
            @assert rule_walker(Rules(3 => "bc"), "abcd", 3, 2) == [4]
            @assert rule_walker(Rules(4 => "bc"), "abcd", 4, 1) == [0]
        end

        @testset "alternative rule" begin
            @assert 5 in rule_walker(Rules(0 => Alternative([1, 2], [2, 1]), 1 => "ab", 2 => "cd"), "abcd", 0, 1)
            @assert rule_walker(Rules(0 => Alternative([1, 2], [2, 1]), 1 => "ab", 2 => "cde"), "abcd", 0, 1) == []
            @assert 5 in rule_walker(Rules(0 => Alternative([1, 2], [2, 3]), 1 => "ab", 2 => "cd", 3 => "cde"), "abcd", 0, 1)
            @assert 6 in rule_walker(Rules(0 => Alternative([1, 2], [1, 3]), 1 => "ab", 2 => "cd", 3 => "cde"), "abcde", 0, 1)
            @assert 7 in rule_walker(Rules(0 => Rule([9, 10]), 9 => Alternative([1, 2], [1, 3]), 1 => "ab", 2 => "cd", 3 => "cde", 10 => "z"), "abcdez", 0, 1)
        end

        @testset "rule" begin
            @assert !(4 in rule_walker(Rules(0 => Rule([1, 2]), 1 => "a", 2 => "b"), "abc", 0, 1))
            @assert rule_walker(Rules(0 => Rule([1, 2]), 1 => "b", 2 => "a"), "ab", 0, 1) == []
            @assert 3 in rule_walker(Rules(0 => Rule([1, 2]), 1 => "a", 2 => "b"), "ab", 0, 1)
            @assert 4 in rule_walker(Rules(0 => Rule([1, 2]), 1 => "a", 2 => Alternative([3], [3, 4]), 3 => "b", 4 => "c"), "abc", 0, 1)
        end

        @testset "shell" begin
            @assert !rule_walker(Rules(0 => Rule([1, 2]), 1 => "b", 2 => "a"), "ab")
            @assert rule_walker(Rules(0 => Rule([1, 2]), 1 => "a", 2 => "b"), "ab")
            @assert !rule_walker(Rules(0 => Rule([1, 2]), 1 => "a", 2 => "b"), "abc")
            @assert rule_walker(Rules(0 => Rule([1, 2]), 1 => "a", 2 => Alternative([3], [3, 4]), 3 => "b", 4 => "c"), "abc")
        end
    end

    @testset "example 1" begin
        @test ex1(example1) == 2
        @test ex1(example2) == 3
    end

    @testset "example 2" begin
        example3 = Messages(copy(example2.rules), example2.messages)
        example3.rules[8] = Alternative([42], [42, 8])
        example3.rules[11] = Alternative([42, 31], [42, 11, 31])

        @test ex2(example1, false) == 2
        @test ex2(example2, false) == 3
        @test ex2(example3, false) == 12
    end

    @testset "results" begin
        @test ex1(messages) == 205
        @test ex2(messages) == 329
    end
end

max_stars_repo_name: CmdQ/AoC


text:
@testset "1108.defanging-an-ip-address.jl" begin
    @test defanging_ip_address("1.1.1.1") == "1[.]1[.]1[.]1"
    @test defanging_ip_address("255.100.50.0") == "255[.]100[.]50[.]0"
end


text:
#=
Copyright (c) 2015, Intel Corporation
All rights reserved.

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:
- Redistributions of source code must retain the above copyright notice,
  this list of conditions and the following disclaimer.
- Redistributions in binary form must reproduce the above copyright notice,
  this list of conditions and the following disclaimer in the documentation
  and/or other materials provided with the distribution.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
THE POSSIBILITY OF SUCH DAMAGE.
=#

# This code was ported from a Python implementation available at
# http://nengo.ca/docs/html/nef_algorithm.html
# Original Python implementation is copyright (c) the author.

using DocOpt
if "--plot" in ARGS
    using PyPlot
end

typealias Float Float64

const dt = Float(0.001)    # simulation time step
const end_t = Float(10.0)  # simulation duration
const t_rc = Float(0.02)   # membrane RC time constant
const t_ref = Float(0.002) # refractory period
const t_pstc = Float(0.1)  # post-synaptic time constant
const N_samples = 100      # number of sample points to use when finding decoders
const rate_A = [25, 75]    # range of maximum firing rates for population A
const rate_B = [50, 100]   # range of maximum firing rates for population B

# the input to the system over time
input(t)=sin(t)

# the function to compute between A and B
funcAB(x)=x*x

uniform(low, high) = rand() * (high - low) + low

function generate_gain_and_bias(count, intercept_low, intercept_high, rate_low, rate_high)
    gain = Array(Float, count)
    bias = Array(Float, count)
    for i in 1:count
        # desired intercept (x value for which the neuron starts firing
        intercept = uniform(intercept_low, intercept_high)
        # desired maximum rate (firing rate when x is maximum)
        rate = uniform(rate_low, rate_high)

        # this algorithm is specific to LIF neurons, but should
        #  generate gain and bias values to produce the desired
        #  intercept and rate
        z = 1.0 / (1-exp((t_ref-(1.0/rate))/t_rc))
        g = (1 - z)/(intercept - 1.0)
        b = 1 - g*intercept
        gain[i] = g
        bias[i] = b
    end
    return gain, bias
end

# a simple leaky integrate-and-fire model, scaled so that v=0 is resting
#  voltage and v=1 is the firing threshold
function run_neurons(input, v, ref)
    lenv = length(v)
    spikes = Array(Bool, lenv)
    # v is updated in the following loop
    for i in 1:lenv
        dV = dt * (input[i]-v[i]) / t_rc    # the LIF voltage change equation
        v[i] += dV
        if v[i] < 0 v[i] = 0 end            # don't allow voltage to go below 0

        if ref[i] > 0                       # if we are in our refractory period
            v[i] = 0                        #   keep voltage at zero and
            ref[i] -= dt                    #   decrease the refractory period
        end

        if v[i] > 1                         # if we have hit threshold
            spikes[i] = true                #   spike
            v[i] = 0                        #   reset the voltage
            ref[i] = t_ref                  #   and set the refractory period
        else
            spikes[i] = false
        end
    end
    return spikes

end

# measure the spike rate of a whole population for a given represented value x
function compute_response(x, encoder, gain, bias, time_limit=0.5)
    N = length(encoder)          # number of neurons
    v = zeros(Float, N)          # voltage
    ref = zeros(Float, N)        # refractory period

    # compute input corresponding to x
    input = Array(Float, N)
    for i in 1:N
        input[i] = x*encoder[i]*gain[i]+bias[i]
        v[i]=uniform(0,1)  # randomize the initial voltage level
    end

    count = zeros(Int, N)    # spike count for each neuron

    # feed the input into the population for a given amount of time
    t = 0
    while t<time_limit
        spikes=run_neurons(input, v, ref)
        for i in 1:length(spikes)
            if spikes[i] count[i]+=1 end
        end
        t += dt
    end
    return [convert(Float, c)/time_limit for c in count]  # return the spike rate (in Hz)
end

# compute the tuning curves for a population
function compute_tuning_curves(encoder, gain, bias)
    N = length(encoder)
    # generate a set of x values to sample at
    x_values = Float[i*2.0/N_samples - 1.0 for i in 0:N_samples-1]

    # build up a matrix of neural responses to each input (i.e. tuning curves)
    A = Array(Float, N, N_samples)
    for i in 1:length(x_values)
        x = x_values[i]
        response = compute_response(x, encoder, gain, bias)
        A[:, i] = response
    end
    return x_values, A
end

# compute decoders
function compute_decoder(encoder, gain, bias, func=x->x)
    # get the tuning curves
    x_values, A = compute_tuning_curves(encoder, gain, bias)

    # get the desired decoded value for each sample point
    # result is N_samples * 1
    value = Float[ func(x) for x in x_values, i in 1:1 ]
    println("size(A) = ", size(A))
    println("size(value) = ", size(value))
    # find the optimum linear decoder
    At = transpose(A)
    # Gamma is N_samples * N_samples
    Gamma = A * At
    # Upsilon is N_samples * 1
    Upsilon = A * value
    Ginv = pinv(Gamma)
    println("size(Ginv) = ", size(Ginv))
    println("size(Upsilon) = ", size(Upsilon))
    decoder = Ginv * Upsilon ./ dt
    return decoder
end

function simulate(NT, inputs, weights, pstc_scale, decoder_B, gain_B, bias_B, encoder_A, gain_A, bias_A, N_A, N_B)

    v_A = zeros(Float, N_A)       # voltage for population A
    ref_A = zeros(Float, N_A)     # refractory period for population A
    input_A = zeros(Float, N_A)   # input for population A

    v_B = zeros(Float, N_B)       # voltage for population B
    ref_B = zeros(Float, N_B)     # refractory period for population B
    input_B = zeros(Float, N_B)   # input for population B
    outputs = zeros(Float, NT)    # outputs
    output = 0.0                  # the decoded output value from population B
    for i = 1:NT
        # call the input function to determine the input value
        x = inputs[i]
        # convert the input value into an input for each neuron
        input_A = x .* encoder_A .* gain_A .+ bias_A

        # run population A and determine which neurons spike
        spikes_A=run_neurons(input_A, v_A, ref_A)
        # for each neuron that spikes, increase the input current
        #  of all the neurons it is connected to by the synaptic
        #  connection weight
        for j = 1:N_A
            if spikes_A[j]
                for k in 1:N_B
                    input_B[k] += weights[k,j] * pstc_scale
                end
            end
        end

        # compute the total input into each neuron in population B
        #  (taking into account gain and bias)
        total_B = gain_B .* input_B .+ bias_B
        input_B = input_B .* (1.0 - pstc_scale)

        # run population B and determine which neurons spike
        spikes_B=run_neurons(total_B, v_B, ref_B)
        # for each neuron in B that spikes, update our decoded value
        #  (also applying the same post-synaptic filter)
        output *= (1.0-pstc_scale)
        for j in 1:length(spikes_B)
            if spikes_B[j]
                output += decoder_B[j] * pstc_scale
            end
        end

        outputs[i] = output
    end
    return outputs
end

function main()
    doc = """nengo-nef.jl

Nengo neural engineering framework.

Usage:
  nengo-nef.jl -h | --help
  nengo-nef.jl [--plot] [--n_a=<n_a>] [--n_b=<n_b>]

Options:
  -h --help    Show this screen.
  --plot       Enable PyPlot drawing.
  --n_a=<n_a>  Number of neurons in first population [default: 500].
  --n_b=<n_b>  Number of neurons in second population [default: 400].
"""
    arguments = docopt(doc)

    N_A = parse(Int, arguments["--n_a"])
    N_B = parse(Int, arguments["--n_b"])
    plotting = arguments["--plot"]

    # gc_enable(false)

    srand(0)

    #################################################
    # Step 1: Initialization
    #################################################

    # create random encoders for the two populations
    encoder_A = Int[(rand(Bool) ? -1 : 1) for i in 1:N_A]
    encoder_B = Int[(rand(Bool) ? -1 : 1) for i in 1:N_B]

    # random gain and bias for the two populations
    gain_A, bias_A = generate_gain_and_bias(N_A, -1, 1, rate_A[1], rate_A[2])
    gain_B, bias_B = generate_gain_and_bias(N_B, -1, 1, rate_B[1], rate_B[2])

    # find the decoders for A and B
    decoder_A=compute_decoder(encoder_A, gain_A, bias_A, funcAB)
    decoder_B=compute_decoder(encoder_B, gain_B, bias_B)

    # compute the weight matrix
    weights = encoder_B * transpose(decoder_A)

    #################################################
    # Step 2: Running the simulation
    #################################################

    # scaling factor for the post-synaptic filter
    pstc_scale=1.0-exp(-dt/t_pstc)

    # for storing simulation data to plot afterward
    NT = convert(Int, end_t / dt)

    times=Float[ dt * t for t = 0:NT-1 ]
    inputs=Float[ input(x) for x in times ]
    ideal=Float[ funcAB(x) for x in inputs ]

    println("Starting simulation")
    tic()
    outputs=simulate(NT, inputs, weights, Float(pstc_scale), decoder_B, gain_B, bias_B, encoder_A, gain_A, bias_A, N_A, N_B)
    println("SELFTIMED ", toq())

    # gc_enable(true)

    if plotting
        #################################################
        # Step 3: Plot the results
        #################################################
        x,A = compute_tuning_curves(encoder_A, gain_A, bias_A)
        x,B = compute_tuning_curves(encoder_B, gain_B, bias_B)
        PyPlot.figure()
        PyPlot.plot(x, transpose(A))
        PyPlot.title("Tuning curves for population A")
        PyPlot.savefig("A.png")

        PyPlot.figure()
        PyPlot.plot(x, transpose(B))
        PyPlot.title("Tuning curves for population B")
        PyPlot.savefig("B.png")

        PyPlot.figure()
        PyPlot.plot(times, inputs, label = "input")
        PyPlot.plot(times, ideal, label= "ideal")
        PyPlot.plot(times, outputs, label = "output")
        PyPlot.title("Simulation results")
        PyPlot.legend()
        PyPlot.savefig("result.png")
    end

    println("checksum: ", sum(outputs))
    return outputs

end

main()


text:
@testset "Spatial" begin

    fileshp = joinpath(dirname(@__FILE__), "data", "SudQC_GCM.shp")
    P = extractpoly(fileshp, n=1)
    @test size(P) == (2,6)
    @test isnan(P[1,1])
    @test isnan(P[2,1])

end


text:
### A Pluto.jl notebook ###
# v0.12.18

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 6b732070-3cb2-11eb-15a0-65fd51108ed7
begin
	using HDF5
	using Plots
	using PlutoUI
	using LaTeXStrings
	using Printf
	import Pda
end

# ╔═╡ 2ea55a90-7e03-11eb-2849-c52e0bfcc500

# ╔═╡ d54000a4-3cc1-11eb-3eb9-e57ecd5f138b
@bind frame PlutoUI.Slider(0:499, show_value=true)

# ╔═╡ 550980e0-3ccb-11eb-0884-03a9d8e523f7
struct ConfigFrame{T}
	id::Int
	x::Vector{T}
	y::Vector{T}
	theta::Vector{T}
	t::T
	tauR::T
	function ConfigFrame(file::String, id::Int)
		t = Pda.get_h5data(file, "config/$(id)/t")
		x = Pda.get_h5data(file, "config/$(id)/x")
		y = Pda.get_h5data(file, "config/$(id)/y")
		theta = Pda.get_h5data(file, "config/$(id)/theta")
		T = eltype(x)
		tauR = Pda.get_h5attr(file, "tauR")
		return new{T}(id, x, y, theta, t, tauR)
	end
end

# ╔═╡ ac910d60-3cc1-11eb-03a7-cb017e973698
begin
	file = "lH1.000alpha1.200.h5"
	t = Pda.get_h5data(file, "config/$(frame)/t")
	y = Pda.get_h5data(file, "config/$(frame)/y")
	loc, n = Pda.density1d(y, bins=100, normalize=true)
	tauR = Pda.get_h5attr(file, "tauR")
end

# ╔═╡ 442a778a-3d06-11eb-1f3f-b5d1e42ffaca
function jump_max(frame1, frame2)
	x0 = Pda.get_h5data(file, "config/$(frame1)/x")
	x1 = Pda.get_h5data(file, "config/$(frame2)/x")
	y0 = Pda.get_h5data(file, "config/$(frame1)/y")
	y1 = Pda.get_h5data(file, "config/$(frame2)/y")
	dx = sum((x1-x0).^2 +(y1-y0).^2)/length(x0)
	#dxm, i = findmax(dx)
	return dx
end

# ╔═╡ d5f96e5a-3cfc-11eb-26de-414f6ec905a8
begin
	x1 = Float64[]
	y1 = Float64[]
	t2 =Float64[]
	i = 91466
	for frame in 0:199
			append!(x1, Pda.get_h5data(file, "config/$(frame)/x")[i])
			append!(y1, Pda.get_h5data(file, "config/$(frame)/y")[i])
			append!(t2, Pda.get_h5data(file, "config/$(frame)/t"))
	end
end

# ╔═╡ 4d58dede-3d04-11eb-03ef-69cfd329dbff
begin
	scatter(x1, y1, label=nothing, framestyle=:box)
	plot!(x1, y1, lw=2, label=nothing)
	xlabel!(L"x")
	ylabel!(L"y")
end

# ╔═╡ 535eb2b6-3d06-11eb-16b1-e3ebcc4e86f2
x1[end]-x1[1]

# ╔═╡ c36512f2-3d04-11eb-3afb-f3f2e62d07da
begin
	scatter(t2, x1, label=nothing, framestyle=:box)
	plot!(t2, x1, lw=2, label=nothing)
	xlabel!(L"t")
	ylabel!(L"x")
end

# ╔═╡ 127bd560-3cc2-11eb-10d0-3f01f597307c
begin
	scatter(loc[2:end-1],n[2:end-1], framestyle=:box, label=nothing)#,yscale=:log10)
	xlabel!(L"y/H")
	ylabel!("number density")
	time_label = @sprintf "t/\\tau_R=%.2f" t/tauR
	title!(latexstring(time_label))
end

# ╔═╡ 92e63448-3ccc-11eb-1acf-3b4309d37ae2
begin
	plot(Pda.get_h5data(file,"x/t")/tauR,
		Pda.get_h5data(file, "x/m"),
		framestyle=:box,
		label=nothing)
	xlabel!(L"t/\tau_R")
	ylabel!(L"<x>")
end

# ╔═╡ edc909f8-3ccc-11eb-2a4c-31536c489723
begin
	scatter(Pda.get_h5data(file,"x/t")[2:end]/tauR,
		Pda.get_h5data(file, "x/v")[2:end],
		framestyle=:box,
		label=nothing,
		xscale=:log10,
		yscale=:log10)
	t1 =Pda.get_h5data(file,"x/t")[2:end]/tauR
	v1 = Pda.get_h5data(file, "x/v")[2:end]
	plot!(t1, t1.^2/t1[1].^2*v1[1], linestyle=:dash, lw=2, label="slope=2")
	plot!(t1, t1.^1/t1[1].^1*v1[1],
		linestyle=:dash,
		lw=2,
		label="slope=1",
		legend=:topleft)
	xlabel!(L"t/\tau_R")
	ylabel!(L"<x^2>")
end

# ╔═╡ 4980b4d6-3d0a-11eb-2340-ed91b8805fc1
v1[751]/v1[748]

# ╔═╡ 2294b9ec-3d0a-11eb-1a1f-7bb34b614a9c
t1[750]

# ╔═╡ 43bdbdfe-3d09-11eb-3e64-497f67bf7e71
function pareto_runtime(U0, tauR)
	taum = (alpha - 1.0) * tauR / alpha
	return taum / pow(rand(), 1.0 / alpha)
end

# ╔═╡ f95566a2-3cfb-11eb-3202-bf60fca4728a
function pareto_update!(x::Vector{Float64},
		y::Vector{Float64},
		theta::Vector{Float64},
		U0::Real,
		tauR::Real,
		alpha::Real)

end

# ╔═╡ Cell order:
# ╠═6b732070-3cb2-11eb-15a0-65fd51108ed7
# ╠═2ea55a90-7e03-11eb-2849-c52e0bfcc500
# ╠═d54000a4-3cc1-11eb-3eb9-e57ecd5f138b
# ╠═442a778a-3d06-11eb-1f3f-b5d1e42ffaca
# ╠═d5f96e5a-3cfc-11eb-26de-414f6ec905a8
# ╠═4d58dede-3d04-11eb-03ef-69cfd329dbff
# ╠═535eb2b6-3d06-11eb-16b1-e3ebcc4e86f2
# ╠═c36512f2-3d04-11eb-3afb-f3f2e62d07da
# ╟─550980e0-3ccb-11eb-0884-03a9d8e523f7
# ╠═ac910d60-3cc1-11eb-03a7-cb017e973698
# ╠═127bd560-3cc2-11eb-10d0-3f01f597307c
# ╠═92e63448-3ccc-11eb-1acf-3b4309d37ae2
# ╠═edc909f8-3ccc-11eb-2a4c-31536c489723
# ╠═4980b4d6-3d0a-11eb-2340-ed91b8805fc1
# ╠═2294b9ec-3d0a-11eb-1a1f-7bb34b614a9c
# ╠═43bdbdfe-3d09-11eb-3e64-497f67bf7e71
# ╠═f95566a2-3cfb-11eb-3202-bf60fca4728a


text:
using DataFrames

"""
    LabelEncoder()
LabelEncoder structure. `LE(label; count=false, decode=false)` Convert labels(like string) to class numbers(encode), and convert class numbers to labels(decode).

# Example
```jldoctest
julia> label = ["Apple", "Apple", "Pear", "Pear", "Lemon", "Apple", "Pear", "Lemon"]
8-element Vector{String}
 "Apple"
 "Apple"
 "Pear"
 "Pear"
 "Lemon"
 "Apple"
 "Pear"
 "Lemon"

julia> LE = LabelEncoder()
LabelEncoder(Dict{Any, Any}())

julia> classes, count = LE(label, count=true) #Encode
([3.0 3.0 … 1.0 2.0], [3.0 2.0 3.0])

julia> LE(classes, decode=true) #Decode
1×18 Matrix{String}:
 "Apple" "Apple" "Pear" "Pear" "Lemon" "Apple" "Pear" "Lemon"
```
"""
mutable struct LabelEncoder
    d::Dict
    LabelEncoder() = new(Dict())
end

function (LE::LabelEncoder)(data; count::Bool=false, decode::Bool=false)
    if decode
        return map(x->LE.d[x], data)
    end
    target = ones(Int64, length(data))
    s = unique(data)
    if count
        c = zeros(Int64, length(s))
    end
    for i in 1 : length(data)
        j = findfirst(isequal(data[i]), s)
        target[i] = j
        if count
            c[j]+=1
        end
    end
    d = Dict()
    for i in 1 : length(s)
        d[i] = s[i]
    end
    LE.d = d
    if count
        return target, c
    end
    return target
end

"""
    OneHotEncoder()
convert data to Ont-Hot format. If you specified `decode`, `data` will be decoded.

# Example
```jldoctest classification
julia> x = [4.9 3.0; 4.6 3.1; 4.4 2.9; 4.8 3.4; 5.1 3.8; 5.4 3.4; 4.8 3.4; 5.2 4.1; 5.5 4.2; 5.5 3.5; 4.8 3.0; 5.1 3.8; 5.0 3.3; 6.4 3.2; 5.7 2.8; 6.1 2.9; 6.7 3.1; 5.6 2.5; 6.3 2.5; 5.6 3.0; 5.6 2.7; 7.6 3.0; 6.4 2.7; 6.4 3.2; 6.5 3.0; 7.7 3.8; 7.2 3.2; 7.2 3.0; 6.3 2.8; 6.1 2.6; 6.3 3.4; 6.0 3.0; 6.9 3.1; 6.7 3.1; 5.8 2.7; 6.8 3.2; 6.3 2.5];#These data are also used to explanations of other functions.

julia> t = [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2];

julia> OHE = OneHotEncoder()
OneHotEncoder()

julia> ct = OHE(t)
37×3 Matrix{Float64}:
 1.0 0.0 0.0
 1.0 0.0 0.0
 1.0 0.0 0.0
 1.0 0.0 0.0
 1.0 0.0 0.0
 1.0 0.0 0.0
 1.0 0.0 0.0
 1.0 0.0 0.0
 1.0 0.0 0.0
 1.0 0.0 0.0
 1.0 0.0 0.0
 1.0 0.0 0.0
 0.0 1.0 0.0
 0.0 1.0 0.0
 0.0 1.0 0.0
 0.0 1.0 0.0
 0.0 1.0 0.0
 ⋮
 0.0 1.0 0.0
 0.0 0.0 1.0
 0.0 0.0 1.0
 0.0 0.0 1.0
 0.0 0.0 1.0
 0.0 0.0 1.0
 0.0 0.0 1.0
 0.0 0.0 1.0
 0.0 0.0 1.0
 0.0 0.0 1.0
 0.0 0.0 1.0
 0.0 0.0 1.0
 0.0 0.0 1.0
 0.0 0.0 1.0
 0.0 0.0 1.0
 0.0 0.0 1.0
 0.0 0.0 1.0

julia> OHE(ct, decode = true)
37-element Vector{Int64}:
 1
 1
 1
 1
 1
 1
 1
 1
 1
 1
 1
 1
 1
 2
 2
 2
 2
 ⋮
 2
 3
 3
 3
 3
 3
 3
 3
 3
 3
 3
 3
 3
 3
 3
 3
 3
```
"""
struct OneHotEncoder end

function (OHE::OneHotEncoder)(data::AbstractVector{T}; decode::Bool=false) where {T}
    if !decode
        data = Int.(data)
        0 in data && @. data+=1
        unqs = unique(data)
        out = zeros(Int64, length(data), length(unqs))
        for (i, unq) in enumerate(unqs)
            out[findall(isequal(unq), data), i] .= 1
        end
        return out
    else
        return [argmax(data[i], dims=2) for i in 1 : size(data, 1)]
    end
end

function dataframe_ohe(data::AbstractVector{T}; prifex::String="") where {T}
    if !isempty(prifex)
        prifex *= "_"
    else
        prifex = "OHE_"
    end
    unqs = unique(data)
    out = DataFrame(Dict([string(prifex,unq)=>zeros(length(data)) for unq in unqs]...))
    for (ind, k) in enumerate(unqs)
        out[findall(isequal(k), data),ind] .= 1
    end
    return out
end

function (OHE::OneHotEncoder)(df::DataFrame, cols::Vector{T} ) where T
    out = deepcopy(df)
    for col in cols
        if T == Int
            prifex = names(df)[col]
        elseif T == Symbol
            prifex = string(col)
        else
            prifex = col
        end
        data = out[!,col]
        out = select(out, Not([col]))
        out = hcat(out, dataframe_ohe(data, prifex=prifex))
    end
    return out
end


text:
#!/usr/bin/env julia

#=
----------------------------------------------------------------------------
timeline generation function forms the basis of any load-generation workflow
----------------------------------------------------------------------------
=#

"""
    timeline(r)

Generate `RheoTimeData` struct with only the time data based on the range provided.

# Example:

`timeline(0:0.1:10)` returns a RheoTimeData with only time data, with the following series: [0.  0.1   0.2  ...  9.9  10].
"""
function timeline(r::R; savelog = true) where R <: AbstractRange
    log = savelog ? [RheoLogItem( (type=:source, funct=:timeline, params=(r,), keywords=()), (comment="timeline created", type=time_only))] : nothing
    RheoTimeData([], [], Vector{RheoFloat}(r), log)
end

"""
    timeline(;t_start::Real=0., t_end::Real, step::Real=(t_end - t_start)/250.)

Generate `RheoTimeData` struct with only the time data based on the keywords provided.

# Arguments

- `t_start`: Starting time, typically 0
- `t_end`: End time
- `step`: Time between sample. Default set in order to provide ≈ 250 time points.
"""
function timeline(;t_start::Real = 0., t_end::Real, step::Real = (t_end - t_start)/250., savelog = true)
    timeline(t_start:step:t_end, savelog = savelog)
end

"""
    frequencyspec(r::R; logscale=true)

Generate `RheoFreqData` struct with only the frequency data distributed on a linear or log scale.

# Arguments

- `r`: range, e.g. 1:0.1:100
- `logscale`: if true, set the range in logscale, e.g. -2:0.1:20  --> 10 values per decade between 10^-2 and 10^2
"""
function frequencyspec(r::R; logscale=true, savelog = true) where R <: AbstractRange
    if logscale
        ω=[10^o for o in r]
    else
        ω=Vector{RheoFloat}(r)
    end
    log = savelog ? [RheoLogItem( (type=:source, funct=:frequencyspec, params=(r,), keywords=(logscale=logscale,)), (comment="frequency range created", type=freq_only))] : nothing
    RheoFreqData([], [], ω, log)
end

"""
    frequencyspec(;ω_start::Real, ω_end::Real, logstep::Real= 0.05)

Generate `RheoFreqData` struct with only the frequency data distributed on a log scale.

# Arguments

- `ω_start`: min frequency
- `ω_end`: max frequency
- `logstep`: step between frequencies in log scale, e.g. 0.1 --> 10 values per decade.
"""
function frequencyspec(;ω_start::Real, ω_end::Real, logstep::Real = 0.05, savelog = true)
    frequencyspec(log10(ω_start):logstep:log10(ω_end), logscale=true, savelog = savelog)
end

#=
------------------------------------------------------------------------
strainfunction and stressfunction can be used to impose strain or stress
respectively. The subsequent convenience functions are passed to them to
generate the appropriate loading.
------------------------------------------------------------------------
=#
"""
    strainfunction(data::RheoTimeData, f::T) where T<:Function
    strainfunction(f::T, data::RheoTimeData) where T<:Function

Accepts a `RheoTimeData` and outputs a new `RheoTimeData` with a strain imposed.
The strain signal is determined by the function provided, which should take
time as its only argument. The original data's time signal is used.

Normally used with a RheoTimeData generated using the `timeline` function.
"""
function strainfunction(data::RheoTimeData, f::T) where T<:Function
    log = data.log === nothing ? nothing : [data.log; RheoLogItem( (type=:process, funct=:strainfunction, params=(f=f,), keywords=()),
                                    (comment="strain function applied to timeline",) ) ]

    return RheoTimeData(data.σ, rheoconvert(map(f, data.t)), data.t, log)
end

function strainfunction(f::T, data::RheoTimeData) where T<:Function
    strainfunction(data,f)
end

"""
    strainfunction!(data::RheoTimeData, f::T) where T<:Function
    strainfunction!(f::T, data::RheoTimeData) where T<:Function

In-place version of `strainfunction`.
"""
function strainfunction!(data::RheoTimeData, f::T) where T<:Function
    data.log === nothing ? nothing : push!(data.log, RheoLogItem( (type=:process, funct=:strainfunction, params=(f=f,), keywords=()),
                                    (comment="strain function applied to timeline",) ) )
    _setstrain!(data, map(f, data.t))
    return data
end

function strainfunction!(f::T, data::RheoTimeData) where T<:Function
    strainfunction!(data,f)
end

"""
    stressfunction(data::RheoTimeData, f::T) where T<:Function
    stressfunction(f::T, data::RheoTimeData) where T<:Function

Accepts a `RheoTimeData` and outputs a new `RheoTimeData` with a stress imposed.
The stress signal is determined by the function provided, which should take
time as its only argument. The original data's time signal is used.

Normally used with a `RheoTimeData` generated using the `timeline` function.
"""
function stressfunction(data::RheoTimeData, f::T) where T<:Function
    log = data.log === nothing ? nothing : [data.log; RheoLogItem( (type=:process, funct=:stressfunction, params=(f=f,), keywords=()),
                                    (comment="stress function applied to timeline",) ) ]
    return RheoTimeData(convert(Vector{RheoFloat}, map(f, data.t)), data.ϵ, data.t, log)
end

function stressfunction(f::T, data::RheoTimeData) where T<:Function
    stressfunction(data,f)
end

"""
    stressfunction!(data::RheoTimeData, f::T) where T<:Function
    stressfunction!(f::T, data::RheoTimeData) where T<:Function

In-place version of `stressfunction`.
"""
function stressfunction!(data::RheoTimeData, f::T) where T<:Function
    data.log === nothing ? nothing : push!(data.log, RheoLogItem( (type=:process, funct=:stressfunction, params=(f=f,), keywords=()),
                                    (comment="stress function applied to timeline",) ) )
    _setstress!(data, map(f, data.t))
    return data
end

function stressfunction!(f::T, data::RheoTimeData) where T<:Function
    stressfunction!(data,f)
end

#=
------------------------------------------------------------------------------
These functions provide convenient tools to design steps, ramps and
other input patterns for the stress and strain data

Functions may also return an anonymous function if the parameter t is omitted.
------------------------------------------------------------------------------
=#
"""
    hstep(t; offset=0., amp=1.)

Step generation function for use with `stressfunction` or `strainfunction`.
`offset` keyword arguent determines start of step. `amp` argument determines
amplitude (height) of step.
"""
function hstep(t; offset=0., amp=1.)
    return (t<offset) ? 0 : amp
end

function hstep(; kwargs...)
    return t -> hstep(t; kwargs...)
end

"""
    ramp(t; offset=0., gradient=1.)

Ramp signal generation function for use with `stressfunction` or `strainfunction`.
`offset` keyword argument determines start of ramp. `gradient` argument determines
the linear gradient of the ramp.
"""
function ramp(t; offset=0., gradient=1.)
    return (t<offset) ? 0 : (t-offset) * gradient
end

function ramp(; kwargs...)
    return t -> ramp(t; kwargs...)
end

"""
    stairs(t; offset=0., amp=1., width=1.)

Stairs signal generation function for use with `stressfunction` or `strainfunction`.
Equivalent to additional steps being added every `width` seconds. `offset` keyword
argument determines start of stairs signal. `amp` argument determines the height
of each additional step.
"""
function stairs(t; offset=0., amp=1., width=1.)
    return (t<offset) ? 0 : amp * floor((t - offset)/width + 1)
end

function stairs(; kwargs...)
    return t -> stairs(t; kwargs...)
end

"""
    square(t; offset=0., amp=1., period=1., width=0.5*period)

Square signal generation function for use with `stressfunction` or `strainfunction`.
`offset` keyword argument determines start of square signal. `amp` argument determines
the height of each square pulse. `period` determines the period of one off/on section
of the square wave signal. `width` determines the width of each square pulse.
"""
function square(t; offset=0., amp=1., period=1., width=0.5*period)
    return t<offset ? 0. : ( ((t-offset)%period) < width ? amp : 0.)
end

function square(; kwargs...)
    return t -> square(t; kwargs...)
end

"""
    sawtooth(t; offset=0., amp=1., period=1.)

Sawtooth signal generation function for use with `stressfunction` or `strainfunction`.
`offset` keyword argument determines start of sawtooth signal. `amp` argument determines
the height of each sawtooth pulse. `period` determines the period of the sawtooth wave signal.
"""
function sawtooth(t; offset=0., amp=1., period=1.)
    return t<offset ? 0. : amp*((t-offset)%period)/period
end

function sawtooth(; kwargs...)
    return t -> sawtooth(t; kwargs...)
end

"""
    triangle(t; offset=0., amp=1., period=1.)

Triangle signal generation function for use with `stressfunction` or `strainfunction`.
`offset` keyword argument determines start of triangle signal. `amp` argument determines
the height of each triangle pulse. `period` determines the period of the triangle wave signal.
`width` determines the width of the triangles.
"""
function triangle(t; offset=0., amp=1., period=1., width=0.5*period)
    if t<offset
        return 0.
    end
    t0=(t-offset)%period
    if t0 <width
        return amp*t0/width
    else
        return amp*(period-t0)/(period-width)
    end
end

function triangle(; kwargs...)
    return t -> triangle(t; kwargs...)
end

#=
------------------------------------------------------------------------------------
Functionality not yet integrated with new data generation format, see GitHub issues:
-----------------------------------------------------------------------------------
"""
    stepgen(t_total::Real, t_on::Real; t_trans::Real = 0.0, stepsize::Real = 1.0)

Generate RheologyData struct with a step loading of height 1.0. If `t_trans` is 0.0 then
the step is instantaneous, otherwise the step is approximated by a logistic function
approximately centered at `t_on`.

# Arguments

- `t_total`: Total time length of data
- `t_on`: Step on time
- `t_trans`: Step transition time
- `stepsize`: Time sampling period
"""
function stepgen(t_total::Real,
                  t_on::Real;
                  t_trans::Real = 0.0,
                  stepsize::Real = 1.0)

    t = collect(0.0:stepsize:t_total)

    data = similar(t)
    # for smooth transition
    if t_trans>0.0
        k = 10.0/t_trans
        data = 1.0./(1 .+ exp.(-k*(t.-t_on)))
    # discrete jump
    elseif t_trans==0.0
        for (i, t_i) in enumerate(t)
            if t_i>=t_on
                data[i] = 1.0
            elseif t_i<t_on
                data[i] = 0.0
            end
        end
    end

    RheologyData(data, t, ["stepgen: t_total: $t_total, t_on: $t_on, t_trans: $t_trans, stepsize: $stepsize"])

end

"""
    noisegen(t_total::Real; seed::Union{Int, Nothing} = nothing, stepsize::Real = 1.0)

Generate uniform random noise of maximum amplitude +/- 1.0. If reproducibility is required,
always use the same number in the `seed` keyword argument with the same non-negative integer.

# Arguments

- `t_total`: Total time length of data
- `seed`: Seed used for random number generation
- `baseval`: Initial amplitude before oscillation started
- `stepsize`: Time sampling period
"""
function noisegen(t_total::Real; seed::Union{Int, Nothing} = nothing, stepsize::Real = 1.0)

    # conditional loading of random as this is the only RHEOS function which requires it
    @eval import Random: seed!, rand

    t = collect(0.0:stepsize:t_total)

    if seed!=nothing
        # use specified seed
        @assert seed>=0 "Seed integer must be non-negative"
        log = ["noisegen: seed: $seed, stepsize: $stepsize"]
        seed!(seed)
    else
        log = ["noisegen: seed: nothing, stepsize: $stepsize"]
    end

    data = (2*rand(eltype(t), length(t)) .- 1)

    RheologyData(data, t, log)

end

"""
    repeatdata(self::RheologyData, n::Integer)

Repeat a given RheologyData data set `n` times.
"""
function repeatdata(self::RheologyData, n::Integer; t_trans = 0.0)

    @assert self.σ==self.ϵ "Repeat data only works when σ==ϵ which is the state of RheologyData after being generated using the built-in RHEOS data generation functions"

    dataraw = self.σ

    @assert constantcheck(self.t) "Data sample-rate must be constant"

    step_size = self.t[2] - self.t[1]

    t = collect(0.0:step_size:(self.t[end]*n))

    # smooth transition between repeats
    if t_trans>0.0
        elbuffer = round(Int, (1/2)*(t_trans/step_size))
        selflength = length(self.t)

        # ensure smooth transition from end of data to new beginning so no discontinuities
        data_smooth_end = stepgen(self.t[end], self.t[end]; t_trans = t_trans, amplitude = dataraw[1] - dataraw[end], baseval = dataraw[selflength - elbuffer], stepsize = step_size)
        data_smooth_start = stepgen(self.t[end], 0.0; t_trans = t_trans, amplitude = dataraw[1] - dataraw[end], baseval = dataraw[selflength - elbuffer], stepsize = step_size)
        data_smoother = data_smooth_end + data_smooth_start

        data_single = self + data_smoother

        data = repeat(data_single.data[1:end] - dataraw[selflength - elbuffer], outer=[n])

        for i = 1:(n-1)
            deleteat!(data, i*length(self.t) + (2-i))
        end

        # fix first repeat
        for (i, v) in enumerate(dataraw)
            if i<(elbuffer*100) && i<round(Int, length(self.t)/2)
                data[i] = v
            end
        end

        log = vcat(self.log, ["repeated data $n times with transition time $t_trans"])

        return RheologyData(data, t, log)

    # discrete jump
    elseif t_trans==0.0

        data = repeat(dataraw, outer=[n])

        for i = 1:(n-1)
            deleteat!(data, i*length(self.t) + (2-i))
        end

        log = vcat(self.log, ["repeated data $n times with transition time $t_trans"])

        return RheologyData(data, t, log)

    end

end
=#


text:
using DFTK
using Test
using LinearAlgebra: norm

include("testcases.jl")

@testset "Computing χ0" begin
    Ecut=3
    fft_size = [10, 1, 10]
    tol=1e-10
    ε = 1e-8
    kgrid = [3, 1, 1]
    testcase = silicon
    n_bands = 12

    spec = ElementPsp(testcase.atnum, psp=load_psp(testcase.psp))
    for temperature in (0, 0.03)
        model = model_LDA(testcase.lattice, [spec => testcase.positions], temperature=temperature)
        basis = PlaneWaveBasis(model, Ecut; kgrid=kgrid, fft_size=fft_size, use_symmetry=false)

        ρ0 = guess_density(basis)
        energies, ham0 = energy_hamiltonian(basis, nothing, nothing; ρ=ρ0)
        ρ1 = DFTK.next_density(ham0, tol=tol, eigensolver=diag_full, n_bands=n_bands).ρ
        χ0 = compute_χ0(ham0)

        # Now we make the same model, but add an artificial external potential ε * dV
        dV = randn(eltype(basis), basis.fft_size)
        term_builder = basis -> DFTK.TermExternal(basis, ε .* dV)
        model = model_LDA(testcase.lattice, [spec => testcase.positions], temperature=temperature, extra_terms=[term_builder])
        basis = PlaneWaveBasis(model, Ecut; kgrid=kgrid, fft_size=fft_size, use_symmetry=false)
        energies, ham = energy_hamiltonian(basis, nothing, nothing; ρ=ρ0)
        ρ2 = DFTK.next_density(ham, tol=tol, eigensolver=diag_full, n_bands=n_bands).ρ
        diff = (ρ2.real - ρ1.real)/ε

        predicted_diff = real(reshape(χ0*vec(dV), basis.fft_size))
        @test norm(diff - predicted_diff) < sqrt(ε)

        for droptol in (0, Inf)
            # Test the diagonal_only option
            χ0_diag = compute_χ0(ham0; droptol=droptol)
            diff_diag_1 = real(reshape(χ0_diag*vec(dV), basis.fft_size))

            EVs = [eigen(Hermitian(Array(Hk))) for Hk in ham0.blocks]
            Es = [EV.values[1:n_bands] for EV in EVs]
            Vs = [EV.vectors[:, 1:n_bands] for EV in EVs]
            occ, εF = find_occupation(basis, Es)
            diff_diag_2 = apply_χ0(ham0, dV, Vs, εF, Es; droptol=droptol,
                                   sternheimer_contribution=(droptol == 0))
            @test norm(diff_diag_1 - diff_diag_2) < sqrt(ε)
        end
    end
end


text:
"""
`Rotate180 <: ImageTransformation`

Description
============

Rotates the image 180 degrees. This is a special case rotation
because it can be performed very efficiently by simply rearranging
the existing pixels. Furthermore, the output images is guaranteed
to have the same dimensions as the input image.

If created using the parameter `chance`, the operation will be
lifted into a `Either`. See the documentation of `Either` for more
information.

Usage
======

    Rotate180()

    Rotate180(chance)

Arguments
==========

- **`chance`** : The probability of the operation to be applied to
the image if the operation is used in `transform`. Must be within
the set [0, 1].

Methods
========

- **`transform`** : Applies the transformation to the given Image
or set of images

Author(s)
==========

- Christof Stocker (Github: https://github.com/Evizero)

Examples
========

    using Images
    using Augmentor

    # Create a test image with 4 pixel.
    img = grayim(UInt8[200 150; 50 1])

    # Apply the image operation to the image
    img_new = transform(Rotate180(), img)

    # Check manually if the result is as expected
    @assert img_new == grayim(UInt8[1 50; 150 200])

see also
=========

`ImageTransformation`, `Either`, `transform`
"""
immutable Rotate180 <: ImageTransformation
end

Rotate180(chance) = ProbableOperation(Rotate180(); chance = chance)

Base.showcompact(io::IO, op::Rotate180) = print(io, "Rotate 180 degrees.")
multiplier(::Rotate180) = 1

function transform{T<:AbstractImage}(op::Rotate180, img::T)
    result = copyproperties(img, rotate_expand(img, 1π))::T
    _log_operation!(result, op)
end


text:
mutable struct Unit
	color
	ir::Integer #influence radius
	pl::Array{Int,1} #permitted layers
	passover::Bool
	passoverself::Bool
	inclusive::Bool
	loc
	groundlevel::Bool
	live!::Function
	baselife::Number
	harvest::Function
	name::String
	costfun::Function
	canspawn::Bool
	harvested
	graphic::Array #of vertices (x,y) where x is distance in radians from origin
	extra #for arbitrary functionality, define your own units and add them to units in units.jl
end
mutable struct Group
	spawns
	body
	units
	#lifmap #I want this
	points
	#allowdebt::Bool #should require sudo password irl, finally got that degree because of pressure from society and now you want to relax to get the creative juices going? Well tough luck, into the hamsterwheel, slave! For being self-proclaimed planetary rulers humen don't seem to have all that much time, baboons can be more with their children. But hey, they have money, gotta give em that, they cut down the trees they need to breathe and caused the largest mass extinction since the dinosaurs to play a game of who has most...
	harvested::Bool
end
mutable struct Board
	shells::Integer #layers of locations to add to the initial ones
	initlocs #initial locations
	grid
	c #GtkCanvas
	sizemod::Number #zoom
	size::Number
	offsetx::Number
	offsety::Number
	panx::Number
	pany::Number
	bgcolor #background
	gridcolor
	#expandbasecost::Number #deprecated because creating new land where our units can live free is not something that should be hindered... You want to make more baby units then go right ahead, we will even help you!
end
mutable struct Game
	name::String #defaults to a timestamp, a good name makes Game feel special and unique
	map #add locations here without changing board.grid to go offgrid!
	groups #units like grouping around babymakers
	spawns #units that can make babies
	units #All units!
	unitparams #for unit to be placed, this is just a string in an array now...
	color #(r,g,b), do not let r+g+b exceed 1, unless you want superpowered units
	colors
	colind::Integer #index of colors -> color
	colmax::Integer #return to first color after max
	colock::Bool #place a single color
	delete::Bool #delete units (by clicking on them)
	sequence::Array{Any} #placed units and performed (harvests/)expands/deletions
	board::Board #changing this after game is initialized causes segfaults
	printscore::Bool #who needs dis when we 'ave fancy labels
	points #not used anymore
	season::Integer #number of harvests, seasons and harvests will become relevant again in phase 2
	win #GtkWindow
	window #initial aspect ratio, funnily named...
	lifemap #deprecated. No! Reprecated
	g #GtkGrid
	gui::Dict #Gtk placeholder
	autoharvest::Bool #costs have been removed so no longer relevant, return to the source
end


text:
import Distributions: pdf, logpdf, NoncentralChisq, LogNormal, Normal

dim(::AbstractSDE{D,M}) where {D,M} = (D,M)
dim(::Type{AbstractSDE{D,M}}) where {D,M} = (D,M)
dim(::Type{T}) where {T<:AbstractSDE} = dim(supertype(T))

model_dim(::AbstractSDE{D,M}) where {D,M} = D
noise_dim(::AbstractSDE{D,M}) where {D,M} = M

function drift end
# function drift_jacobian end
function diffusion end
function jump end
function mark_distribution end
function variables end

_reshape_or_nothing(x::Number, dims) = x
_reshape_or_nothing(x::AbstractArray, dims) = reshape(x, dims)

function corrected_drift(model::AbstractSDE{D,M}, t0, x0, α) where {D,M}
  μ = drift(model, t0, x0)
  σ = diffusion(model, t0, x0)
  corr = @MArray zeros(D)
  ∂σ = _reshape_or_nothing(_diffusion_jacobian(model, t0, x0), (D, M, D))
  # μ - α*vec(sum(permutedims(∂σ, (3,2,1)).*σ, dims=(1,2)))
  for j=1:M, i=1:D, k=1:D
    corr[i] += σ[k,j] * ∂σ[i,j,k]
  end
  μ - α*corr
end

include("codegen.jl")

include("black_scholes.jl")
include("cox_ingersoll_ross.jl")
include("merton.jl")
include("ornstein_uhlenbeck.jl")

@sde_model FitzHughNagumo begin
  dX = ɛ*(X-X^3-Y+s)*dt + σ1*dW1
  dY =     (γ*X-Y+β)*dt + σ2*dW2
end

@sde_model Heston begin
  dS =     r*S*dt + sqrt(V)*S*dW1
  dV = κ*(θ-V)*dt + σ*sqrt(V)*dW2
  dW1*dW2 = ρ*dt
end

@sde_model Lorenz begin
  dx = σ*(y-x)*dt + σ1*dw1
  dy = (x*(ρ-z)-y)*dt + σ2*dw2
  dz = (x*y-β*z)*dt + σ3*dw3
end

@sde_model Bates begin
  dS =     α*S*dt + sqrt(V)*S*dW1 + S*(ξ-1.0)*dN
  dV = κ*(θ-V)*dt + σ*sqrt(V)*dW2
  dW1*dW2 = ρ*dt
  ξ ~ LogNormal(μ, δ)
end α κ θ σ ρ λ μ δ

@sde_model CompoundPoisson begin
  dS = S*(ξ-1.0)*dN
  ξ ~ LogNormal(μ, δ)
end λ μ δ


text:
"""
    simulate(model::Model, rule::Rule, replicates::Int, db::Sqlite.DB; batchname::String, seed::Int)

Run a model `replicates` number of times with a given `rule` for every possible sequence of meetings.
"""
function simulate(
    model::Model, rule::Rule, replicates::Int, db::SQLite.DB;
    batchname::String, seed::Int = 1
)
    register_model!(db, model, batchname)
    Random.seed!(seed)
    sequences = permutations(collect(combinations(model.parameter_set.parties, 2)))
    @showprogress 1 "Running simulations..." for (seq_idx, seq) in enumerate(sequences)
        results_data = @chain begin
            run_model_on_sequence(model, rule, seq, replicates)
            snap(_, :seq, seq_idx)
            snap(_, :batchname, batchname)
        end
        SQLite.load!(results_data, db, "results")
    end
    sequences_data = @chain begin
        format_sequences_for_db(sequences)
        snap(_, :batchname, batchname)
    end
    SQLite.load!(sequences_data, db, "sequences")
    return true
end

"""
    register_model!(db::SQLite.DB, model::Model, batchname::String)

Register a given model in the simulation database.
"""
function register_model!(db::SQLite.DB, model::Model, batchname::String)
    parliament_table = stack(DataFrame(sort(collect(model.parameter_set.parliament))), Not([]))
    rename!(parliament_table, [:party_shorthand, :seats])
    parliament_table[!, :batchname] .= batchname
    parameters_table = DataFrame(
        batchname = batchname,
        group_size = model.parameter_set.group_size,
        parliament_majority = model.parameter_set.parliament_majority,
        required_consensus = model.parameter_set.required_consensus
    )
    SQLite.load!(parliament_table, db, "parliament")
    SQLite.load!(parameters_table, db, "parameters")
    return true
end

"""
    run_sequence(sequence::AbstractArray, replicates::Int)

Run the model `replicates` times for a given `sequence`.
"""
function run_model_on_sequence(model::Model, rule::Rule, sequence::AbstractArray, replicates::Int)
    replicate_list = DataFrame[]
    # TODO: figure out how much can be parallelized and how
    @sync @distributed for rep in 1:replicates
        model_tracker = deepcopy(model)
        rep_data = snap(DataFrame(deepcopy(model_tracker.agents)), :step, 0)  # track initial configuration
        for (step, comb) in enumerate(sequence)
            meeting = Meeting(model_tracker, comb)  # we are in a meeting now
            rule.update!(meeting)
            step_data = DataFrame(deepcopy(model_tracker.agents))
            rep_data = reduce(vcat, [rep_data, snap(step_data, :step, step)])
        end
        push!(replicate_list, snap(deepcopy(rep_data), :rep, rep))
    end
    seq_data = reduce(vcat, replicate_list)
    seq_data = format_results_for_db(seq_data)
    return seq_data
end

"""
    format_results_for_db(data::DataFrame)

Format the simulation data for storage in the database.
"""
function format_results_for_db(data::DataFrame)
    reshaped_array = []
    for i in 1:length(data.opinions[1])
        current_statement = []
        for j in 1:length(data.opinions)
            push!(current_statement, data.opinions[j][i])
        end
        push!(reshaped_array, deepcopy(current_statement))
    end
    left_side = select(data, Not(:opinions))
    right_side = @chain begin
        DataFrame(reshaped_array, :auto)
        rename(_, Symbol.(1:ncol(_)))
    end
    data_formatted = @chain begin
        hcat(left_side, right_side)
        stack(_, 5:ncol(_))  # TODO: not ideal, better with pattern matching by column name?
        rename(_, Dict(:id => :agent_id, :variable => :statement_id, :value => :position))
    end
    return data_formatted
end

"""
    format_sequences_for_db(sequences)

Formats a sequence generator for storage in a result database.
"""
function format_sequences_for_db(sequences)
    df = DataFrame()
    for (i, seq) in enumerate(sequences)
        for (j, mtg) in enumerate(seq)
            push!(df, (seq_id = i, step = j, party_1 = mtg[1], party_2 = mtg[2]))
        end
    end
    return df
end

"""
    snap(data::DataFrame, scope::Symbol, val::Union{Int, String})

Make a snapshot of the simulation data at `val` at `scope`.
"""
function snap(data::DataFrame, scope::Symbol, val::Union{Int, String})
    data[!, scope] .= val
    return data
end


text:
using DataDeps
using BSON

@testset "Custom layers" begin
    @testset "WeightDroppedLSTM" begin
        wd = ULMFiT.WeightDroppedLSTM(4, 5, 0.3)
        @test all(wd.init .== wd.state)
        @test size(wd.cell.Wi) == size(wd.cell.maskWi)
        @test size(wd.cell.Wh) == size(wd.cell.maskWh)
        @test wd.cell.active
        @test_throws DimensionMismatch wd(rand(5, 3))
        x = rand(Float32, 4, 3)
        h = wd(x)
        @test size(h) == size(wd.state[1]) == (5, 3)
        maski = deepcopy(wd.cell.maskWi)
        maskh = deepcopy(wd.cell.maskWh)
        ULMFiT.reset_masks!(wd)
        @test maski != wd.cell.maskWi
        @test maskh != wd.cell.maskWh
        Flux.testmode!(wd)
        @test !(wd.cell.active)
        @test length(params(wd)) == 5
    end

    @testset "AWD_LSTM" begin
        awd = ULMFiT.AWD_LSTM(4, 3, 0.3)
        @test awd.layer.cell isa ULMFiT.WeightDroppedLSTMCell
        @test length(awd.accum) == 0
        ULMFiT.set_trigger!(4, awd)
        @test awd.T == 4
        ULMFiT.asgd_step!(4, awd)
        @test length(awd.accum) == 3
        temp = deepcopy(awd.accum[1][1])
        @test temp == Tracker.data(awd.layer.cell.Wi[1])
        ULMFiT.asgd_step!(5, awd)
        temp += temp
        @test temp == Tracker.data(awd.accum[1][1])
        @test length(params(awd)) == 5
    end

    @testset "VarDrop" begin
        vd = ULMFiT.VarDrop(0.3)
        @test vd.active
        @test vd.reset
        x = rand(10, 10)
        @test size(vd(x)) == size(x) == size(vd.mask)
        x = rand(5, 5)
        @test_throws DimensionMismatch vd(x)
        ULMFiT.reset_masks!(vd)
        @test vd.reset
        Flux.testmode!(vd)
        @test ~vd.active
        x = rand(10, 10)
        @test vd(x) == x
        @test size(vd(x)) == size(x) == size(vd.mask)
    end

    @testset "DroppedEmbeddings" begin
        de = ULMFiT.DroppedEmbeddings(6, 4, 0.2)
        @test size(de.emb) == (6, 4)
        @test size(de.mask) == (6,)
        x = [2,4,6,0.1]
        @test_throws BoundsError de(x)
        x = [2,4,6]
        @test size(de(x)) == (4, 3)
        x = rand(5, 3)
        @test_throws DimensionMismatch de(x, true)
        x = rand(4, 3)
        @test size(de(x, true)) == (6, 3)
        mask = de.mask
        ULMFiT.reset_masks!(de)
        mask != de.mask
        Flux.testmode!(de)
        @test ~de.active
        @test length(params(de)) == 1
    end

    @testset "PooledDense" begin
        pd = ULMFiT.PooledDense(10, 5)
        @test size(pd.W) == (5, 30)
        @test length(pd.b) == 5
        x = rand(Float32, 10, 3)
        @test_throws DimensionMismatch pd(x)
        @test size(pd([x])) == (5, 3)
        @test length(params(pd)) == 2
    end
end

@testset "Language model" begin
    lm = ULMFiT.LanguageModel()
    @test typeof(lm.vocab) == Vector{String}
    @test length(lm.vocab) == size(lm.layers[1].emb, 1)
    @test length(lm.layers) == 10
    @test length(params(lm)) == 16
    @test length(ULMFiT.get_trainable_params(lm.layers)) == 10

    pretrained_weights = BSON.load(datadep"Pretrained ULMFiT Language Model/ulmfit_lm_en.bson")
    @test length(pretrained_weights[:weights]) == 16
    @test all(size.(params(lm)) .== size.(pretrained_weights[:weights]))
end

@testset "Text Classifier" begin
    lm = ULMFiT.LanguageModel()
    tc = ULMFiT.TextClassifier(lm)
    @test tc.rnn_layers == lm.layers[1:8]
    @test length(tc.linear_layers) == 6
end


text:
# Modeling
@batchdef function time_modeling_azure(_model, srcG, q, recG, D, _dm, op, mode, opt)
    model = fetch(_model)
    dm = isnothing(_dm) ? _dm : fetch(_dm)
    argout = time_modeling(model, srcG, q, recG, D, dm, op, mode, opt)
    return argout
end

#FWI
@batchdef function fwi_objective_azure(_model, q, D, opt)
    model = fetch(_model)
    argout = fwi_objective(model, q, D, opt)
    return argout
end

#LSRTM
@batchdef function lsrtm_objective_azure(_model, q, D, _dm, opt; nlind=false)
    model = fetch(_model)
    dm = isnothing(_dm) ? _dm : fetch(_dm)
    argout = lsrtm_objective(model,q, D, dm, opt; nlind=nlind)
    return argout
end