text:
using DiffEqBase
import DiffEqBase: AbstractODESolution, HermiteInterpolation,AbstractDiffEqInterpolation
using DiffEqSensitivity
using ForwardDiff

abstract type AbstractSimulation end
export AbstractSimulation

struct Simulation{Q<:AbstractODESolution,W<:AbstractDomain,M,L<:AbstractArray,G<:Function,G2<:AbstractArray,G3,G4,G5} <: AbstractSimulation
    sol::Q
    domain::W
    interfaces::M
    names::L
    N::G
    Ns::G2
    species::G3
    reactions::G4
    p::G5
end

function Simulation(sol::Q,domain::W,interfaces=[],p=nothing) where {Q<:AbstractODESolution,W<:AbstractDomain}
    names = getfield.(domain.phase.species,:name)
    Ns = sum(hcat(sol.interp.u...)[domain.indexes[1]:domain.indexes[2],:],dims=1)
    if hasproperty(sol.interp,:du)
        Nderivs = sum(hcat(sol.interp.du...)[domain.indexes[1]:domain.indexes[2],:],dims=1)
    else
        Nderivs = sum(hcat([sol(t,Val{1}) for t in sol.t]...)[domain.indexes[1]:domain.indexes[2],:],dims=1)
    end
    N = HermiteInterpolation(sol.interp.t,Ns,Nderivs)
    F(t::T) where {T<:Real} = N(t,nothing,Val{0},sol.prob.p,:left)
    return Simulation(sol,domain,interfaces,names,F,Ns,domain.phase.species,domain.phase.reactions,p)
end

export Simulation

struct SystemSimulation{Q,B<:AbstractODESolution,X,Y,Z}
    sol::B
    sims::Q
    interfaces::X
    names::Array{String,1}
    species::Array{Z,1}
    reactions::Array{Y,1}
    p::Array{Float64,1}
end

function SystemSimulation(sol,domains,interfaces,p)
    sims = Tuple([Simulation(sol,domain) for domain in domains])
    names = Array{String,1}()
    reactions = Array{ElementaryReaction,1}()
    species = Array{Species,1}()
    for sim in sims
        append!(names,sim.names)
        append!(species,sim.domain.phase.species)
        append!(reactions,sim.domain.phase.reactions)
    end
    for inter in interfaces
        if hasproperty(inter,:reactions)
            append!(reactions,inter.reactions)
        end
    end
    return SystemSimulation(sol,sims,interfaces,names,species,reactions,p)
end
export SystemSimulation

length(p::T) where {T<:AbstractSimulation} = 1
export length

iterate(p::T) where {T<:AbstractSimulation} = p
export iterate

Broadcast.broadcastable(p::T) where {T<:AbstractSimulation} = Ref(p)
export broadcastable

spcindex(bsol::Z,name::Q) where {Z<:Simulation,Q<:AbstractString} = findfirst(isequal(name),getfield.(bsol.domain.phase.species,:name))
export spcindex

function molefractions(bsol::Q,name::W,t::E) where {Q<:AbstractSimulation, W<:String, E<:Real}
    @assert name in bsol.names
    ind = findfirst(isequal(name),bsol.names)+bsol.domain.indexes[1]-1
    return bsol.sol(t)[ind]/bsol.N(t)
end

function molefractions(bsol::Q, t::E) where {Q<:AbstractSimulation,E<:Real}
    return bsol.sol(t)[bsol.domain.indexes[1]:bsol.domain.indexes[2]]./bsol.N(t)
end

function molefractions(bsol::Q) where {Q<:AbstractSimulation}
    @views return hcat(bsol.sol.u...)[bsol.domain.indexes[1]:bsol.domain.indexes[2],:]./bsol.Ns
end

export molefractions

function concentrations(bsol::Q,name::W,t::E) where {Q<:AbstractSimulation, W<:String, E<:Real}
    @assert name in bsol.names
    ind = findfirst(isequal(name),bsol.names)+bsol.domain.indexes[1]-1
    return bsol.sol(t)[ind]/getdomainsize(bsol,t)
end

function concentrations(bsol::Q, t::E) where {Q<:AbstractSimulation,E<:Real}
    return bsol.sol(t)[bsol.domain.indexes[1]:bsol.domain.indexes[2]]./getdomainsize(bsol,t)
end

function concentrations(bsol::Q) where {Q<:AbstractSimulation}
    @views return hcat(bsol.sol.u...)[bsol.domain.indexes[1]:bsol.domain.indexes[2],:]./getdomainsize.(bsol,bsol.sol.t)'
end

function concentrations(ssys::Q,name::W,t::E) where {Q<:SystemSimulation, W<:String, E<:Real}
    @assert name in ssys.names
    for sim in ssys.sims
        if name in sim.names
            return concentrations(sim,name,t)
        end
    end
end

function concentrations(ssys::Q, t::E) where {Q<:SystemSimulation,E<:Real}
    cstot = zeros(length(ssys.species))
    for sim in ssys.sims
        cstot[sim.domain.indexes[1]:sim.domain.indexes[2]] .= concentrations(sim,t)
    end
    return cstot
end

function concentrations(ssys::Q) where {Q<:SystemSimulation}
    cstots = zeros((length(ssys.species),length(ssys.sol.t)))
    for (i,t) in enumerate(ssys.sol.t)
        for sim in ssys.sims
            cstots[sim.domain.indexes[1]:sim.domain.indexes[2],i] .= concentrations(sim,t)
        end
    end
    return cstots
end

export concentrations

getT(bsol::Simulation{Q,W,L,G}, t::K) where {W<:Union{ConstantTPDomain,ConstantTVDomain},K<:Real,Q,G,L} = bsol.domain.T
getT(bsol::Simulation{Q,W,L,G}, t::K) where {W<:Union{ConstantVDomain,ParametrizedVDomain,ConstantPDomain,ParametrizedPDomain},K<:Real,Q,G,L} = bsol.sol(t)[bsol.domain.indexes[3]]
getT(bsol::Simulation{Q,W,L,G}, t::K) where {W<:Union{ParametrizedTConstantVDomain,ParametrizedTPDomain},K<:Real,Q,G,L} = bsol.domain.T(t)
export getT
getV(bsol::Simulation{Q,W,L,G}, t::K) where {W<:Union{ConstantVDomain,ConstantTVDomain,ParametrizedTConstantVDomain},K<:Real,Q,G,L} = bsol.domain.V
getV(bsol::Simulation{Q,W,L,G}, t::K) where {W<:ParametrizedVDomain,K<:Real,Q,G,L} = bsol.domain.V(t)
getV(bsol::Simulation{Q,W,L,G}, t::K) where {W<:Union{ConstantTPDomain,ParametrizedTPDomain},K<:Real,Q,G,L} = bsol.sol(t)[bsol.domain.indexes[3]]
getV(bsol::Simulation{Q,W,L,G}, t::K) where {W<:Union{ParametrizedPDomain,ConstantPDomain},K<:Real,Q,G,L} = bsol.sol(t)[bsol.domain.indexes[4]]
export getV
getP(bsol::Simulation{Q,W,L,G}, t::K) where {W<:Union{ConstantTPDomain,ConstantPDomain},K<:Real,Q,G,L} = bsol.domain.P
getP(bsol::Simulation{Q,W,L,G}, t::K) where {W<:Union{ConstantTVDomain,ParametrizedTConstantVDomain},K<:Real,Q,G,L} = 1.0e6
getP(bsol::Simulation{Q,W,L,G}, t::K) where {W<:Union{ParametrizedTPDomain,ParametrizedPDomain},K<:Real,Q,G,L} = bsol.domain.P(t)
getP(bsol::Simulation{Q,W,L,G}, t::K) where {W<:Union{ConstantVDomain, ParametrizedVDomain},K<:Real,Q,G,L} = bsol.sol(t)[bsol.domain.indexes[4]]
export getP
getC(bsol::Simulation{Q,W,L,G}, t::K) where {W<:ConstantTPDomain,K<:Real,Q,G,L} = bsol.domain.P/(R*bsol.domain.T)
getC(bsol::Simulation{Q,W,L,G}, t::K) where {W<:Union{ConstantVDomain,ConstantTVDomain,ParametrizedTConstantVDomain},K<:Real,Q,G,L} = bsol.N(t)/bsol.domain.V
getC(bsol::Simulation{Q,W,L,G}, t::K) where {W<:ParametrizedVDomain,K<:Real,Q,G,L} = bsol.N(t)/bsol.domain.V(t)
getC(bsol::Simulation{Q,W,L,G}, t::K) where {W<:ParametrizedTPDomain,K<:Real,Q,G,L} = bsol.domain.P(t)/(R*bsol.domain.T(t))
getC(bsol::Simulation{Q,W,L,G}, t::K) where {W<:ConstantPDomain,K<:Real,Q,G,L} = bsol.domain.P/(R*getT(bsol,t))
getC(bsol::Simulation{Q,W,L,G}, t::K) where {W<:ParametrizedPDomain,K<:Real,Q,G,L} = bsol.domain.P(t)/(R*getT(bsol,t))
export getC
getdomainsize(bsol,t) = getV(bsol,t)
getdomainsize(bsol::Simulation{Q,W,L,G},t) where {W<:ConstantTAPhiDomain,Q,L,G} = bsol.domain.A

"""
calculates the rates of production/loss at a given time point
this outputs a sparse matrix of  num reactions xnum species containing the production/loss
rate of that species associated with that reaction
"""
function rops(bsol::Q,t::X) where {Q<:Simulation,X<:Real}
    ropmat = spzeros(length(bsol.domain.phase.reactions),length(bsol.domain.phase.species))
    cs,kfs,krevs = calcthermo(bsol.domain,bsol.sol(t),t)[[2,9,10]]
    V = getdomainsize(bsol,t)
    @simd for i in 1:length(bsol.domain.phase.reactions)
        rxn = bsol.domain.phase.reactions[i]
        R = getrate(rxn,cs,kfs,krevs)*V
        for ind in rxn.productinds
            ropmat[i,ind] += R
        end
        for ind in rxn.reactantinds
            ropmat[i,ind] -= R
        end
    end
    return ropmat
end

function rops(ssys::SystemSimulation,t)
    domains = getfield.(ssys.sims,:domain)
    Nrxns = sum([length(sim.domain.phase.reactions) for sim in ssys.sims])+sum([length(inter.reactions) for inter in ssys.interfaces if hasproperty(inter,:reactions)])
    Nspcs = sum([length(sim.domain.phase.species) for sim in ssys.sims])
    cstot = zeros(Nspcs)
    vns = Array{Any,1}(undef,length(domains))
    vcs = Array{Any,1}(undef,length(domains))
    vT = Array{Any,1}(undef,length(domains))
    vP = Array{Any,1}(undef,length(domains))
    vV = Array{Any,1}(undef,length(domains))
    vC = Array{Any,1}(undef,length(domains))
    vN = Array{Any,1}(undef,length(domains))
    vmu = Array{Any,1}(undef,length(domains))
    vkfs = Array{Any,1}(undef,length(domains))
    vkrevs = Array{Any,1}(undef,length(domains))
    vHs = Array{Any,1}(undef,length(domains))
    vUs = Array{Any,1}(undef,length(domains))
    vGs = Array{Any,1}(undef,length(domains))
    vdiffs = Array{Any,1}(undef,length(domains))
    vCvave = Array{Any,1}(undef,length(domains))
    vphi = Array{Any,1}(undef,length(domains))
    ropmat = spzeros(Nrxns,Nspcs)
    start = 1
    for (k,sim) in enumerate(ssys.sims)
        vns[k],vcs[k],vT[k],vP[k],vV[k],vC[k],vN[k],vmu[k],vkfs[k],vkrevs[k],vHs[k],vUs[k],vGs[k],vdiffs[k],vCvave[k],vphi[k] = calcthermo(sim.domain,ssys.sol(t),t)
        cstot[sim.domain.indexes[1]:sim.domain.indexes[2]] = vcs[k]
        rops!(ropmat,sim.domain.rxnarray,cstot,vkfs[k],vkrevs[k],vV[k],start)
        start += length(vkfs[k])
    end
    for inter in ssys.interfaces
        if hasproperty(inter,:reactions)
            kfs,krevs=getkfskrevs(inter,vT[inter.domaininds[1]],vT[inter.domaininds[2]],vphi[inter.domaininds[1]],vphi[inter.domaininds[2]],vGs[inter.domaininds[1]],vGs[inter.domaininds[2]],cstot)
            rops!(ropmat,inter.rxnarray,cstot,kfs,krevs,inter.A,start)
            start += length(kfs)
        end
    end
    return ropmat
end

"""
calculates the rates of production/loss at a given time point for a single species
this outputs a sparse vector of length num reactions containing the production/loss
rate associated with that reaction for the given species
"""
function rops(bsol::Y,name::X,t::Z) where {Y<:Simulation, X<:AbstractString, Z<:Real}
    rop = spzeros(length(bsol.domain.phase.reactions))
    cs,kfs,krevs = calcthermo(bsol.domain,bsol.sol(t),t)[[2,9,10]]
    V = getdomainsize(bsol,t)
    ind = findfirst(isequal(name),getfield.(bsol.domain.phase.species,:name))
    @assert !isa(ind,Nothing) "species $name not in species array"
    for (i,rxn) in enumerate(bsol.domain.phase.reactions)
        c = 0
        R = getrate(rxn,cs,kfs,krevs)*V
        c -= count(isequal(ind),rxn.reactantinds)
        c += count(isequal(ind),rxn.productinds)
        if c != 0
            rop[i] = c*R
        end
    end
    return rop
end

function rops(ssys::SystemSimulation,name,t)
    domains = getfield.(ssys.sims,:domain)
    ind = findfirst(isequal(name),ssys.names)
    Nrxns = sum([length(sim.domain.phase.reactions) for sim in ssys.sims])+sum([length(inter.reactions) for inter in ssys.interfaces if hasproperty(inter,:reactions)])
    Nspcs = sum([length(sim.domain.phase.species) for sim in ssys.sims])
    cstot = zeros(Nspcs)
    vns = Array{Any,1}(undef,length(domains))
    vcs = Array{Any,1}(undef,length(domains))
    vT = Array{Any,1}(undef,length(domains))
    vP = Array{Any,1}(undef,length(domains))
    vV = Array{Any,1}(undef,length(domains))
    vC = Array{Any,1}(undef,length(domains))
    vN = Array{Any,1}(undef,length(domains))
    vmu = Array{Any,1}(undef,length(domains))
    vkfs = Array{Any,1}(undef,length(domains))
    vkrevs = Array{Any,1}(undef,length(domains))
    vHs = Array{Any,1}(undef,length(domains))
    vUs = Array{Any,1}(undef,length(domains))
    vGs = Array{Any,1}(undef,length(domains))
    vdiffs = Array{Any,1}(undef,length(domains))
    vCvave = Array{Any,1}(undef,length(domains))
    vphi = Array{Any,1}(undef,length(domains))
    ropvec = spzeros(Nrxns)
    start = 0
    for (k,sim) in enumerate(ssys.sims)
        vns[k],vcs[k],vT[k],vP[k],vV[k],vC[k],vN[k],vmu[k],vkfs[k],vkrevs[k],vHs[k],vUs[k],vGs[k],vdiffs[k],vCvave[k],vphi[k] = calcthermo(sim.domain,ssys.sol(t),t)
        cstot[sim.domain.indexes[1]:sim.domain.indexes[2]] = vcs[k]
        rops!(ropvec,sim.domain.rxnarray,cstot,vkfs[k],vkrevs[k],vV[k],start,ind)
        start += length(vkfs[k])
    end
    for inter in ssys.interfaces
        if hasproperty(inter,:reactions)
            kfs,krevs=getkfskrevs(inter,vT[inter.domaininds[1]],vT[inter.domaininds[2]],vphi[inter.domaininds[1]],vphi[inter.domaininds[2]],vGs[inter.domaininds[1]],vGs[inter.domaininds[2]],cstot)
            rops!(ropvec,inter.rxnarray,cstot,kfs,krevs,inter.A,start,ind)
            start += length(kfs)
        end
    end
    return ropvec
end

export rops

function rops!(ropmat,rarray,cs,kfs,krevs,V,start)
    for i = 1:length(kfs)
        if @inbounds rarray[2,i] == 0
            @inbounds @fastmath fR = kfs[i]*cs[rarray[1,i]]
        elseif @inbounds rarray[3,i] == 0
            @inbounds @fastmath fR = kfs[i]*cs[rarray[1,i]]*cs[rarray[2,i]]
        else
            @inbounds @fastmath fR = kfs[i]*cs[rarray[1,i]]*cs[rarray[2,i]]*cs[rarray[3,i]]
        end
        if @inbounds rarray[5,i] == 0
            @inbounds @fastmath rR = krevs[i]*cs[rarray[4,i]]
        elseif @inbounds rarray[6,i] == 0
            @inbounds @fastmath rR = krevs[i]*cs[rarray[4,i]]*cs[rarray[5,i]]
        else
            @inbounds @fastmath rR = krevs[i]*cs[rarray[4,i]]*cs[rarray[5,i]]*cs[rarray[6,i]]
        end
        @fastmath R = (fR - rR)*V

        @inbounds @fastmath ropmat[i+start,rarray[1,i]] -= R
        if @inbounds rarray[2,i] != 0
            @inbounds @fastmath ropmat[i+start,rarray[2,i]] -= R
            if @inbounds rarray[3,i] != 0
                @inbounds @fastmath ropmat[i+start,rarray[3,i]] -= R
            end
        end
        @inbounds @fastmath ropmat[i+start,rarray[4,i]] += R
        if @inbounds rarray[5,i] != 0
            @inbounds @fastmath ropmat[i+start,rarray[5,i]] += R
            if @inbounds rarray[6,i] != 0
                @inbounds @fastmath ropmat[i+start,rarray[6,i]] += R
            end
        end
    end
end

function rops!(ropvec,rarray,cs,kfs,krevs,V,start,ind)
    for i = 1:length(kfs)
        c = count(isequal(ind),rarray[4:6,i])-count(isequal(ind),rarray[1:3,i])
        if c != 0.0
            if @inbounds rarray[2,i] == 0
                @inbounds @fastmath fR = kfs[i]*cs[rarray[1,i]]
            elseif @inbounds rarray[3,i] == 0
                @inbounds @fastmath fR = kfs[i]*cs[rarray[1,i]]*cs[rarray[2,i]]
            else
                @inbounds @fastmath fR = kfs[i]*cs[rarray[1,i]]*cs[rarray[2,i]]*cs[rarray[3,i]]
            end
            if @inbounds rarray[5,i] == 0
                @inbounds @fastmath rR = krevs[i]*cs[rarray[4,i]]
            elseif @inbounds rarray[6,i] == 0
                @inbounds @fastmath rR = krevs[i]*cs[rarray[4,i]]*cs[rarray[5,i]]
            else
                @inbounds @fastmath rR = krevs[i]*cs[rarray[4,i]]*cs[rarray[5,i]]*cs[rarray[6,i]]
            end
            @fastmath R = (fR - rR)*V
            @fastmath @inbounds ropvec[i+start] = c*R
        end
    end
end

"""
Calculates sensitivities with respect to `target` at the time point at the end of the simulation
The returned sensitivities are the normalized values

By default uses the InterpolatingAdjoint algorithm with vector Jacobian products calculated with ReverseDiffVJP(true)
this assumes no changes in code branching during simulation, if that were to become no longer true, the Tracker
based alternative algorithm is slower, but avoids this concern.
"""
function getadjointsensitivities(bsol::Q,target::String,solver::W;sensalg::W2=InterpolatingAdjoint(autojacvec=ReverseDiffVJP(false)),abstol::Float64=1e-6,reltol::Float64=1e-3,kwargs...) where {Q,W,W2}
    @assert target in bsol.names || target in ["T","V","P"]

    pethane = 160

    if target in ["T","V","P"]
        if haskey(bsol.domain.thermovariabledict, target)
            ind = bsol.domain.thermovariabledict[target]
        else
            throw(error("$(bsol.domain) doesn't have $target in its thermovariables"))
        end
    else
        ind = findfirst(isequal(target),bsol.names)
        if isempty(bsol.interfaces)
            sensdomain,sensspcnames,senstooriginspcind,senstooriginrxnind = getsensdomain(bsol.domain,ind)
            if :thermovariabledict in fieldnames(typeof(bsol.domain))
                yinds = vcat(senstooriginspcind,collect(values(bsol.domain.thermovariabledict)))
            else
                yinds = vcat(senstooriginspcind)
            end
            pinds = vcat(senstooriginspcind,length(bsol.domain.phase.species).+senstooriginrxnind)
            ind = findfirst(isequal(target),sensspcnames)
        end
    end

    function sensg(y::X,p::Array{Y,1},t::Z) where {Q,V,X,Y<:Float64,Z}
        sensy = y[yinds]
        sensp = p[pinds]
        dy = similar(sensy,length(sensy))
        return dydtreactor!(dy,sensy,t,sensdomain,[],p=sensp)[ind]
    end
    function sensg(y::Array{X,1},p::Y,t::Z) where {Q,V,X<:Float64,Y,Z}
        sensy = y[yinds]
        sensp = p[pinds]
        dy = similar(sensp,length(sensy))
        return dydtreactor!(dy,sensy,t,sensdomain,[],p=sensp)[ind]
    end
    function sensg(y::Array{X,1},p::Array{Y,1},t::Z) where {Q,V,X<:Float64,Y<:Float64,Z}
        sensy = y[yinds]
        sensp = p[pinds]
        dy = similar(sensy,length(sensy))
        return dydtreactor!(dy,sensy,t,sensdomain,[],p=sensp)[ind]
    end
    function sensg(y::Array{X,1},p::Array{Y,1},t::Z) where {Q,V,X<:ForwardDiff.Dual,Y<:ForwardDiff.Dual,Z}
        sensy = y[yinds]
        sensp = p[pinds]
        dy = similar(sensy,length(sensy))
        return dydtreactor!(dy,sensy,t,sensdomain,[],p=sensp)[ind]
    end

    function g(y::X,p::Array{Y,1},t::Z) where {Q,V,X,Y<:Float64,Z}
        dy = similar(y,length(y))
        return dydtreactor!(dy,y,t,bsol.domain,bsol.interfaces,p=p)[ind]
    end
    function g(y::Array{X,1},p::Y,t::Z) where {Q,V,X<:Float64,Y,Z}
        dy = similar(p,length(y))
        return dydtreactor!(dy,y,t,bsol.domain,bsol.interfaces,p=p)[ind]
    end
    function g(y::Array{X,1},p::Array{Y,1},t::Z) where {Q,V,X<:Float64,Y<:Float64,Z}
        dy = zeros(length(y))
        return dydtreactor!(dy,y,t,bsol.domain,bsol.interfaces,p=p)[ind]
    end
    function g(y::Array{X,1},p::Array{Y,1},t::Z) where {Q,V,X<:ForwardDiff.Dual,Y<:ForwardDiff.Dual,Z}
        dy = similar(y,length(y))
        return dydtreactor!(dy,y,t,bsol.domain,bsol.interfaces,p=p)[ind]
    end

    dsensgdu(out, y, p, t) = ForwardDiff.gradient!(out, y -> sensg(y, p, t), y)
    dsensgdp(out, y, p, t) = ForwardDiff.gradient!(out, p -> sensg(y, p, t), p)
    dgdu(out, y, p, t) = ForwardDiff.gradient!(out, y -> g(y, p, t), y)
    dgdp(out, y, p, t) = ForwardDiff.gradient!(out, p -> g(y, p, t), p)
    dsensgdurevdiff(out, y, p, t) = ReverseDiff.gradient!(out, y -> sensg(y, p, t), y)
    dsensgdprevdiff(out, y, p, t) = ReverseDiff.gradient!(out, p -> sensg(y, p, t), p)
    dgdurevdiff(out, y, p, t) = ReverseDiff.gradient!(out, y -> g(y, p, t), y)
    dgdprevdiff(out, y, p, t) = ReverseDiff.gradient!(out, p -> g(y, p, t), p)

    if length(bsol.domain.p)<= pethane
        if target in ["T","V","P"] || !isempty(bsol.interfaces)
            du0,dpadj = adjoint_sensitivities(bsol.sol,solver,g,nothing,(dgdu,dgdp);sensealg=sensalg,abstol=abstol,reltol=reltol,kwargs...)
        else
            du0,dpadj = adjoint_sensitivities(bsol.sol,solver,sensg,nothing,(dsensgdu,dsensgdp);sensealg=sensalg,abstol=abstol,reltol=reltol,kwargs...)
        end
    else
        if target in ["T","V","P"] || !isempty(bsol.interfaces)
            du0,dpadj = adjoint_sensitivities(bsol.sol,solver,g,nothing,(dgdurevdiff,dgdprevdiff);sensealg=sensalg,abstol=abstol,reltol=reltol,kwargs...)
        else
            du0,dpadj = adjoint_sensitivities(bsol.sol,solver,sensg,nothing,(dsensgdurevdiff,dsensgdprevdiff);sensealg=sensalg,abstol=abstol,reltol=reltol,kwargs...)
        end
    end
    dpadj[length(bsol.domain.phase.species)+1:end] .*= bsol.domain.p[length(bsol.domain.phase.species)+1:end]
    if !(target in ["T","V","P"]) && isempty(bsol.interfaces)
        dpadj ./= bsol.sol(bsol.sol.t[end])[senstooriginspcind[ind]]
    elseif !(target in ["T","V","P"]) && !isempty(bsol.interfaces)
        dpadj ./= bsol.sol(bsol.sol.t[end])[ind]
    end
    return dpadj
end

function getadjointsensitivities(syssim::Q,bsol::W3,target::String,solver::W;sensalg::W2=InterpolatingAdjoint(autojacvec=ReverseDiffVJP(false)),abstol::Float64=1e-6,reltol::Float64=1e-3,kwargs...) where {Q,W,W2,W3}
    @assert target in bsol.names || target in ["T","V","P"]
    if target in ["T","V","P"]
        if haskey(bsol.domain.thermovariabledict, target)
            ind = bsol.domain.thermovariabledict[target]
        else
            throw(error("$(bsol.domain) doesn't have $target in its thermovariables"))
        end
    else
        ind = findfirst(isequal(target),bsol.names)+bsol.domain.indexes[1]-1
    end
    domains = Tuple([x.domain for x in syssim.sims])
    function g(y::X,p::Array{Y,1},t::Z) where {Q,V,X,Y<:Float64,Z}
        dy = similar(y,length(y))
        return dydtreactor!(dy,y,t,domains,syssim.interfaces,p=p)[ind]
    end
    function g(y::Array{X,1},p::Y,t::Z) where {Q,V,X<:Float64,Y,Z}
        dy = similar(p,length(y))
        return dydtreactor!(dy,y,t,domains,syssim.interfaces,p=p)[ind]
    end
    function g(y::Array{Float64,1},p::Array{Float64,1},t::Z) where {Q,V,Z}
        dy = similar(p,length(y))
        return dydtreactor!(dy,y,t,domains,syssim.interfaces,p=p)[ind]
    end
    function g(y::Array{X,1},p::Array{Y,1},t::Z) where {Q,V,X<:ForwardDiff.Dual,Y<:ForwardDiff.Dual,Z}
        dy = similar(y,length(y))
        return dydtreactor!(dy,y,t,domains,syssim.interfaces,p=p)[ind]
    end
    dgdu(out, y, p, t) = ForwardDiff.gradient!(out, y -> g(y, p, t), y)
    dgdp(out, y, p, t) = ForwardDiff.gradient!(out, p -> g(y, p, t), p)
    du0,dpadj = adjoint_sensitivities(syssim.sol,solver,g,nothing,(dgdu,dgdp);sensealg=sensalg,abstol=abstol,reltol=reltol,kwargs...)
    for domain in domains
       dpadj[domain.parameterindexes[1]+length(domain.phase.species):domain.parameterindexes[2]] .*= syssim.p[domain.parameterindexes[1]+length(domain.phase.species):domain.parameterindexes[2]]
    end
    if !(target in ["T","V","P"])
       dpadj ./= bsol.sol(bsol.sol.t[end])[ind]
    end
    return dpadj
end
export getadjointsensitivities

function getconcentrationsensitivity(bsol::Simulation{Q,W,L,G}, numerator::String, denominator::String, t::K) where {W<:Union{ConstantVDomain,ConstantTVDomain,ParametrizedTConstantVDomain},K<:Real,Q,G,L}
    @assert numerator in bsol.names
    @assert denominator in bsol.names
    indnum = findfirst(isequal(numerator),bsol.names)+bsol.domain.indexes[1]-1
    inddeno = findfirst(isequal(denominator),bsol.names)+bsol.domain.parameterindexes[1]-1
    Nvars = length(bsol.domain.phase.species)+length(bsol.domain.indexes)-2
    Nrxns = length(bsol.domain.phase.reactions)
    x,dp = extract_local_sensitivities(bsol.sol,t)
    s = dp[inddeno][indnum]
    val = s/bsol.sol(t)[indnum] #constant volume
    if t == 0
        return 0.0
    else
        return val
    end
end

function getconcentrationsensitivity(bsol::Simulation{Q,W,L,G}, numerator::String, denominator::String, t::K) where {W<:Union{ConstantTPDomain,ParametrizedTPDomain,ConstantPDomain,ParametrizedPDomain,ParametrizedVDomain},K<:Real,Q,G,L}
    @assert numerator in bsol.names
    @assert denominator in bsol.names
    indnum = findfirst(isequal(numerator),bsol.names)+bsol.domain.indexes[1]-1
    inddeno = findfirst(isequal(denominator),bsol.names)+bsol.domain.parameterindexes[1]-1
    Nvars = length(bsol.domain.phase.species)+length(bsol.domain.indexes)-2
    Nrxns = length(bsol.domain.phase.reactions)
    x,dp = extract_local_sensitivities(bsol.sol,t)
    svals = dp[inddeno][bsol.domain.indexes[1]:bsol.domain.indexes[2]]
    s = svals[indnum]
    V = getV(bsol,t)
    c = bsol.sol(t)[indnum]/V
    val =  (s-c*sum(svals)*R*getT(bsol,t)/getP(bsol,t))/(c*V) #known T and P
    if t == 0
        return 0.0
    else
        return val
    end
end

function getconcentrationsensitivity(bsol::Simulation{Q,W,L,G}, numerator::String, denominator::Z, t::K) where {W<:Union{ConstantVDomain,ConstantTVDomain,ParametrizedTConstantVDomain},K<:Real,Z<:Integer,Q,G,L}
    @assert numerator in bsol.names
    indnum = findfirst(isequal(numerator),bsol.names)+bsol.domain.indexes[1]-1
    inddeno = denominator+bsol.domain.parameterindexes[1]-1
    Nvars = length(bsol.domain.phase.species)+length(bsol.domain.indexes)-2
    Nrxns = length(bsol.domain.phase.reactions)
    x,dp = extract_local_sensitivities(bsol.sol,t)
    s = dp[inddeno+length(bsol.domain.phase.species)][indnum]
    T = getT(bsol,t)
    P = getP(bsol,t)
    C = getC(bsol,t)
    k = bsol.domain.p[inddeno+length(bsol.domain.phase.species)]
    val = s*k/bsol.sol(t)[indnum] #constant volume
    if t == 0
        return 0.0
    else
        return val
    end
end

function getconcentrationsensitivity(bsol::Simulation{Q,W,L,G}, numerator::String, denominator::Z, t::K) where {W<:Union{ConstantTPDomain,ParametrizedTPDomain,ConstantPDomain,ParametrizedPDomain,ParametrizedVDomain},K<:Real,Z<:Integer,Q,G,L}
    @assert numerator in bsol.names
    indnum = findfirst(isequal(numerator),bsol.names)+bsol.domain.indexes[1]-1
    inddeno = denominator+bsol.domain.parameterindexes[1]-1
    Nvars = length(bsol.domain.phase.species)+length(bsol.domain.indexes)-2
    Nrxns = length(bsol.domain.phase.reactions)
    x,dp = extract_local_sensitivities(bsol.sol,t)
    svals = dp[inddeno+length(bsol.domain.phase.species)][bsol.domain.indexes[1]:bsol.domain.indexes[2]]
    s = svals[indnum]
    V = getV(bsol,t)
    T = getT(bsol,t)
    P = getP(bsol,t)
    C = getC(bsol,t)
    c = bsol.sol(t)[indnum]/V
    k = bsol.domain.p[inddeno+length(bsol.domain.phase.species)]
    val = k*(s-c*sum(svals)*R*T/P)/(c*V) #known T and P
    if t == 0
        return 0.0
    else
        return val
    end
end

export getconcentrationsensitivity

"""
calculate the rates of all reactions at time t
"""
function rates(bsol::Q,t::X) where {Q<:Simulation,X<:Real}
    cs,kfs,krevs = calcthermo(bsol.domain,bsol.sol(t),t)[[2,9,10]]
    V = getdomainsize(bsol,t)
    return [getrate(rxn,cs,kfs,krevs)*V for rxn in bsol.domain.phase.reactions]
end

"""
calculate the rates of all reactions at given times ts
defaults to using bsol.sol.t if ts is not supplied
"""
function rates(bsol::Q;ts::X=Array{Float64,1}()) where {Q<:Simulation,X<:AbstractArray}
    if length(ts) == 0
        ts = bsol.sol.t
    end
    return hcat([rates(bsol,t) for t in ts]...)
end

"""
calculate the rates of all reactions at time t
"""
function rates(ssys::Q,t::X) where {Q<:SystemSimulation,X<:Real}
    rts = zeros(length(ssys.reactions))
    domains = getfield.(ssys.sims,:domain)
    Nrxns = sum([length(sim.domain.phase.reactions) for sim in ssys.sims])+sum([length(inter.reactions) for inter in ssys.interfaces if hasproperty(inter,:reactions)])
    Nspcs = sum([length(sim.domain.phase.species) for sim in ssys.sims])
    cstot = zeros(Nspcs)
    vns = Array{Any,1}(undef,length(domains))
    vcs = Array{Any,1}(undef,length(domains))
    vT = Array{Any,1}(undef,length(domains))
    vP = Array{Any,1}(undef,length(domains))
    vV = Array{Any,1}(undef,length(domains))
    vC = Array{Any,1}(undef,length(domains))
    vN = Array{Any,1}(undef,length(domains))
    vmu = Array{Any,1}(undef,length(domains))
    vkfs = Array{Any,1}(undef,length(domains))
    vkrevs = Array{Any,1}(undef,length(domains))
    vHs = Array{Any,1}(undef,length(domains))
    vUs = Array{Any,1}(undef,length(domains))
    vGs = Array{Any,1}(undef,length(domains))
    vdiffs = Array{Any,1}(undef,length(domains))
    vCvave = Array{Any,1}(undef,length(domains))
    vphi = Array{Any,1}(undef,length(domains))
    index = 1
    for (k,sim) in enumerate(ssys.sims)
        vns[k],vcs[k],vT[k],vP[k],vV[k],vC[k],vN[k],vmu[k],vkfs[k],vkrevs[k],vHs[k],vUs[k],vGs[k],vdiffs[k],vCvave[k],vphi[k] = calcthermo(sim.domain,ssys.sol(t),t)
        cstot[sim.domain.indexes[1]:sim.domain.indexes[2]] = vcs[k]
        rts[index:index+length(vkfs[k])-1] .= getrates(sim.domain.rxnarray,cstot,vkfs[k],vkrevs[k]).*getdomainsize(sim,t)
        index += length(vkfs[k])
    end
    for inter in ssys.interfaces
        if hasproperty(inter,:reactions)
            kfs,krevs=getkfskrevs(inter,vT[inter.domaininds[1]],vT[inter.domaininds[2]],vphi[inter.domaininds[1]],vphi[inter.domaininds[2]],vGs[inter.domaininds[1]],vGs[inter.domaininds[2]],cstot)
            rts[index:index+length(kfs)-1] = getrates(inter.rxnarray,cstot,kfs,krevs).*inter.A
            index += length(kfs)
        end
    end
    return rts
end

"""
calculate the rates of all reactions at given times ts
defaults to using bsol.sol.t if ts is not supplied
"""
function rates(ssys::Q;ts::X=Array{Float64,1}()) where {Q<:SystemSimulation,X<:AbstractArray}
    if length(ts) == 0
        ts = ssys.sol.t
    end
    return hcat([rates(ssys,t) for t in ts]...)
end

export rates


text:
"""
    Revolve

This is a Julia adaptation of the functionality of Revolve; see Alg. 799 published by Griewank et al.
A minor extension is the  optional `bundle` parameter that allows to treat as many loop
iterations in one tape/adjoint sweep. If `bundle` is 1, the default, then the behavior is that of Alg. 799.

"""
mutable struct Revolve <: Scheme
    steps::Int
    bundle::Int
    tail::Int
    acp::Int
    cstart::Int
    cend::Int
    numfwd::Int
    numinv::Int
    numstore::Int
    rwcp::Int
    prevcend::Int
    firstuturned::Bool
    stepof::Vector{Int}
    verbose::Int
    fstore::Union{Function,Nothing}
    frestore::Union{Function,Nothing}
end

function Revolve(
    steps::Int,
    checkpoints::Int,
    fstore::Union{Function,Nothing} = nothing,
    frestore::Union{Function,Nothing} = nothing;
    anActionInstance::Union{Nothing,Action} = nothing,
    bundle_::Union{Nothing,Int} = nothing,
    verbose::Int = 0
)
    if !isa(anActionInstance, Nothing)
        # same as default init above
        anActionInstance.actionflag = 0
        anActionInstance.iteration  = 0
        anActionInstance.cpNum      = 0
    end
    !isa(bundle_, Nothing) ? bundle = bundle_ : bundle = 1
    if bundle < 1 || bundle > steps
       error("Revolve: bundle parameter out of range [1,steps]")
    elseif steps<0
       error("Revolve: negative steps")
    elseif checkpoints < 0
       error("Revolve: negative checkpoints")
    end
    cstart = 0
    tail   = 1
    if bundle > 1
        tail = mod(steps, bundle)
        steps = steps / bundle
        if tail > 0
            step += 1
        else
            tail = bundle
        end
    end
    cend            = steps
    acp             = checkpoints
    numfwd          = 0
    numinv          = 0
    numstore        = 0
    rwcp            = -1
    prevcend        = 0
    firstuturned    = false
    stepof = Vector{Int}(undef, acp+1)

    revolve = Revolve(steps, bundle, tail, acp, cstart, cend, numfwd, numinv, numstore, rwcp, prevcend, firstuturned, stepof, verbose, fstore, frestore)

    if verbose > 0
        predfwdcnt = forwardcount(revolve)
        if predfwdcnt == -1
            error("Revolve: error returned by  revolve::forwardcount")
        else
            @info "prediction:"
            @info " overhead forward steps : $predfwdcnt"
            @info " overhead factor        : $(predfwdcnt/steps)"
        end
    end
    return revolve
end

function adjust(::Revolve)
    error("Not implemented")
end

function next_action!(revolve::Revolve)::Action
    # Default values for next action
    actionflag     = none
    iteration      = 0
    startiteration = 0
    cpnum          = 0
    if revolve.numinv == 0
        # first invocation
        for v in revolve.stepof
            v = 0
        end
        revolve.stepof[1] = revolve.cstart - 1
    end
    prevcstart = revolve.cstart
    revolve.numinv += 1
    rwcptest = (revolve.rwcp == -1)
    if !rwcptest
       rwcptest = revolve.stepof[revolve.rwcp+1] != revolve.cstart
    end
    if (revolve.cend - revolve.cstart) == 0
       # nothing in current subrange
        if (revolve.rwcp == -1) || (revolve.cstart == revolve.stepof[1])
            # we are done
            revolve.rwcp = revolve.rwcp - 1
            if revolve.verbose > 2
                @info "done"
            end
            if revolve.verbose > 0
                @info "summary:"
                @info " overhead forward steps: $(revolve.numfwd)"
                @info " CP stores             : $(revolve.numstore)"
                @info " NextAction calls      : $(revolve.numinv)"
            end
            actionflag = done
        else
           revolve.cstart = revolve.stepof[revolve.rwcp+1]
           revolve.prevcend = revolve.cend
           actionflag = restore
        end
    elseif (revolve.cend - revolve.cstart) == 1
        revolve.cend = revolve.cend - 1
        revolve.prevcend = revolve.cend
        if (revolve.rwcp >= 0) && (revolve.stepof[revolve.rwcp + 1] == revolve.cstart)
            revolve.rwcp -= 1
        end
        if !revolve.firstuturned
            actionflag = firstuturn
            revolve.firstuturned = true
        else
            actionflag = uturn
        end
    elseif rwcptest
        revolve.rwcp += 1
        if revolve.rwcp + 1 > revolve.acp
            error("Revolve: insufficient allowed checkpoints")
        else
            revolve.stepof[revolve.rwcp+1] = revolve.cstart
            revolve.numstore += 1
            revolve.prevcend = revolve.cend
            actionflag = store
        end
    elseif (revolve.prevcend < revolve.cend) && (revolve.acp == revolve.rwcp + 1)
            error("Revolve: insufficient allowed checkpoints")
    else
        availcp = revolve.acp - revolve.rwcp
        if availcp < 1
            error("Revolve: insufficient allowed checkpoints")
        else
            reps = 0
            range = 1
            while range < (revolve.cend - revolve.cstart)
                reps = reps + 1
                range = range * (reps + availcp) / reps
            end
            bino1 = range * reps / (availcp + reps)
            if availcp > 1
                bino2 = bino1 * availcp / (availcp + reps - 1)
            else
                bino2 = 1
            end
            if availcp==1
                bino3 = 0
            elseif availcp>2
                bino3 = bino2 * (availcp - 1) / (availcp + reps - 2)
            else
                bino3 = 1
            end
            bino4 = bino2 * (reps - 1) / availcp
            if availcp < 3
                bino5 = 0
            elseif availcp > 3
                    bino5 = bino3 * (availcp - 1) / reps
            else
                bino5 = 1
            end
            if (revolve.cend - revolve.cstart) <= (bino1 + bino3)
                revolve.cstart = trunc(Int, revolve.cstart + bino4)
            elseif (revolve.cend - revolve.cstart) >= (range-bino5)
                revolve.cstart = trunc(Int, revolve.cstart + bino1)
            else
                revolve.cstart = trunc(Int, revolve.cend - bino2 - bino3)
            end
            if revolve.cstart == prevcstart
                revolve.cstart = prevcstart + 1
            end
            if revolve.cstart == revolve.steps
                revolve.numfwd = (revolve.numfwd + ((revolve.cstart - 1) - prevcstart)
                                    * revolve.bundle + revolve.tail)
            else
                revolve.numfwd = revolve.numfwd + (revolve.cstart - prevcstart) * revolve.bundle
            end
            actionflag = forward
        end
    end
    startiteration = prevcstart * revolve.bundle
    if actionflag == firstuturn
        iteration = revolve.cstart * revolve.bundle + revolve.tail
    elseif actionflag == uturn
        iteration = (revolve.cstart + 1) * revolve.bundle
    else
        iteration = revolve.cstart * revolve.bundle
    end
    if revolve.verbose > 2
        if actionflag == forward
            @info " run forward iterations    [$startiteration, $(iteration - 1)]"
        elseif actionflag == restore
            @info " restore input of iteration $iteration"
        elseif actionflag == firstuturn

            @info " 1st uturn for iterations  [$startiteration, $(iteration - 1)]"
        elseif actionflag == uturn
            @info " uturn for iterations      [$startiteration, $(iteration - 1)]"
        end
    end
    if (revolve.verbose > 1) && (actionflag == store)
        @info " store input of iteration $iteration  "
    end
    cpnum=revolve.rwcp

    return Action(actionflag, iteration, startiteration, cpnum)

end

function guess(revolve::Revolve; bundle::Union{Nothing, Int} = nothing)::Int
    b=1
    bSteps=revolve.steps
    if !isa(bundle, Nothing)
        b=bundle
    end
    if revolve.steps < 1
        error("Revolve: error: steps < 1")
    elseif b<1
        error("Revolve: error: bundle < 1")
    else
        if b > 1
            revolve.tail = mod(bSteps, b)
            bSteps = div(bSteps, b)
            if revolve.tail > 0
                bSteps = bSteps + 1
            end
        end
        if bSteps == 1
            guess = 0
        else
            checkpoints = 1
            reps = 1
            s = 0
            while chkrange(revolve, checkpoints+s, reps+s) > bSteps
                s -= 1
            end
            while chkrange(revolve, checkpoints+s, reps+s) < bSteps
                s += 1
            end
            checkpoints += s
            reps += s
            s = -1
            while chkrange(revolve, checkpoints, reps) >= bSteps
                if checkpoints > reps
                    checkpoints -= 1
                    s = 0
                else
                    reps -= 1
                    s = 1
                end
            end
            if s == 0
                checkpoints += 1
            end
            if s == 1
                reps += 1
            end
            guess = checkpoints
        end
    end
    return guess
end

function factor(revolve::Revolve, steps, checkpoints, bundle::Union{Nothing,Int} = nothing)
    b = 1
    if !isa(bundle, Nothing)
        b = bundle
    end
    f = forwardcount(revolve)
    if f == -1
        error("Revolve: error returned by forwardcount")
    else
        factor = float(f)/steps
    end
    return factor
end

function chkrange(::Revolve, ss, tt)
    ret = Int(0)
    res = 1.
    if tt < 0 || ss < 0
        error("Revolve chkrange: error: negative parameter")
    else
        for i in 1:tt
            res = res * (ss + i)
            res = res / i
            if res > typemax(typeof(ret))
                break
            end
        end
        if res < typemax(typeof(ret))
            ret = trunc(Int, res)
        else
            ret = typemax(typeof(ret))
            @warn "Revolve chkrange: warning: returning maximal integer "
        end
    end
    return ret
end

function forwardcount(revolve::Revolve)
    checkpoints = revolve.acp
    bundle      = revolve.bundle
    steps       = revolve.steps
    if checkpoints < 0
        error("Revolve forwardcount: error: checkpoints < 0")
    elseif steps < 1
        error("Revolve forwardcount: error: steps < 1")
    elseif bundle < 1
        error("Revolve forwardcount: error: bundle < 1")
    else
        s=steps
        if bundle > 1
            tail = mod(s,bundle)
            s = s / bundle
            if tail > 0
            s = s + 1
            end
        end
        if s == 1
            ret = 0
        elseif checkpoints == 0
            error("Revolve forarwdCount: error: given inputs require checkpoints > 0")
        else
            reps = 0
            range = 1
            while range < s
            reps = reps + 1
            range = range*(reps+checkpoints)/reps
            end
            ret = (reps * s - range * reps / (checkpoints + 1)) * bundle
        end
    end
    return ret
end

function checkpoint_struct(body::Function,
        alg::Revolve,
        model_input::MT,
        shadowmodel::MT
    ) where {MT}
    model = deepcopy(model_input)
    storemap = Dict{Int32,Int32}()
    check = 0
    model_check = Array{MT}(undef, alg.acp)
    model_final = []
    while true
        next_action = next_action!(alg)
        if (next_action.actionflag == Checkpointing.store)
            check = check+1
            storemap[next_action.iteration-1]=check
            model_check[check] = deepcopy(model)
        elseif (next_action.actionflag == Checkpointing.forward)
            for j= next_action.startiteration:(next_action.iteration - 1)
                body(model)
            end
        elseif (next_action.actionflag == Checkpointing.firstuturn)
            body(model)
            model_final =  deepcopy(model)
            Enzyme.autodiff(body, Duplicated(model,shadowmodel))
        elseif (next_action.actionflag == Checkpointing.uturn)
            Enzyme.autodiff(body, Duplicated(model,shadowmodel))
            if haskey(storemap,next_action.iteration-1-1)
                delete!(storemap,next_action.iteration-1-1)
                check=check-1
            end
        elseif (next_action.actionflag == Checkpointing.restore)
            model = deepcopy(model_check[storemap[next_action.iteration-1]])
        elseif next_action.actionflag == Checkpointing.done
            if haskey(storemap,next_action.iteration-1-1)
                delete!(storemap,next_action.iteration-1-1)
                check=check-1
            end
            break
        end
    end
    return model_final
end

ext: jl


text:
# Definition, constructors and other operations only to do with the SACtr type

# Composite type for SAC evenly-spaced time series data
@eval mutable struct SACtr
    $([:($(s)::SACFloat) for s in sac_float_hdr]...)
    $([:($(s)::SACInt) for s in sac_int_hdr]...)
    $([:($(s)::SACBool) for s in sac_bool_hdr]...)
    $([:($(s)::SACChar) for s in sac_char_hdr]...)
    # The time series, accessed with .t
    t::Array{SACFloat,1}
    function SACtr(delta_in::Number, npts_in::Integer, b_in=0.)
        delta_in > 0 || error("SACtr: delta must be positive")
        npts_in >= 0 || error("SACtr: npts must be 0 or larger")
        # Variables are by default undefined, or false for bools
        $([:($(s) = sac_rnull) for s in sac_float_hdr]...)
        $([:($(s) = sac_inull) for s in sac_int_hdr]...)
        $([:($(s) = false) for s in sac_bool_hdr]...)
        $([:($(s) = sac_cnull) for s in sac_char_hdr]...)
        # Variables which must be present
        npts = convert(SACInt, npts_in)
        delta = convert(SACFloat, delta_in)
        b = b_in
        e = b + (npts - 1)*delta
        t = zeros(npts)
        depmin = 0.
        depmax = 0.
        nvhdr = sac_ver_num
        iftype = 1
        idep = 5
        iztype = 9
        ievtyp = 5
        leven = true
        lovrok = true
        lcalda = true
        new($([:($(s)) for s in [sac_float_hdr; sac_int_hdr; sac_bool_hdr; sac_char_hdr]]...),
            t)
    end
end
@doc """
    SACtr(delta, npts, b=0.) -> ::SACtr

Construct a composite type holding an evenly-spaced SAC time-series trace, where the trace
is accessed through the field name `t`.  Supply the constant sampling interval `delta`
in seconds, and the number of points in the trace `t`.  Optionally, specify the trace
start time `b` in seconds.

    SACtr(v::AbstractVector, delta, b=0.) -> ::SACtr

Construct a `SACtr` by supplying an array `v`, sampling interval `delta` and optionally
the starting time.

    SACtr(d::Vector{UInt8}, file=""; swap=true, terse=false, check_npts=true) -> ::SACtr

Construct a SACtr from a raw array of bytes representing some data in SAC format.
If `swap` is false, then non-native-endian files are not converted.  If `terse` is
true, then warnings about swapping are not written.  If `check_npts` is false, then
parts of files are read without error.
""" SACtr

function SACtr(data::Vector{UInt8}, file=""; swap::Bool=true, terse::Bool=false,
        check_npts::Bool=true)
    len = sac_byte_len
    clen = 2*sac_byte_len
    # Determine endianness and act accordingly
    nvhdr = reinterpret(SACInt, data[sac_nvhdr_pos*len+1:(sac_nvhdr_pos+1)*len])[1]
    native = if nvhdr == sac_ver_num
        true
    elseif bswap(nvhdr) == sac_ver_num
        false
    else
        error("Array does not appear to be SAC data")
    end
    native && machine_is_little_endian && !swap &&
        error("Data are little-endian but `swap` is `false`.  Not attempting to swap bytes" *
            (file!="" ? " for file '$file'." : "."))
    native && machine_is_little_endian && !terse &&
        @info("Data are little-endian; byteswapping")
    byteswap(x) = native ? x : bswap(x)
    # Create an empty object
    npts_in_file = (length(data) - sac_header_len)÷len
    trace = SACtr(1, npts_in_file)

    ## Read header
    # Float part
    off = 0
    for (i, field) in enumerate(sac_float_hdr)
        setfield!(trace, field,
            byteswap(reinterpret(SACFloat, data[off+1:off+len])[1]))
        off += len
    end
    # Int part
    for (i, field) in enumerate(sac_int_hdr)
        setfield!(trace, field,
            byteswap(reinterpret(SACInt, data[off+1:off+len])[1]))
        off += len
    end
    # Boolean part
    for (i, field) in enumerate(sac_bool_hdr)
        setfield!(trace, field,
            1 == byteswap(reinterpret(SACInt, data[off+1:off+len])[1]))
        off += len
    end
    # Character part
    # kevnm header is double length, so treat separately
    trace.kstnm = String(strip(String(data[off+1:off+clen])))
    off += clen
    trace.kevnm = String(strip(String(data[off+1:off+2clen])))
    trace.kevnm == strip(rpad(sac_cnull, clen)^2) && (trace.kevnm = sac_cnull)
    off += 2clen
    for (i, field) in enumerate(sac_char_hdr)
        i <= 2 && continue
        setfield!(trace, field,
            String(strip(String(data[1+off:clen+off]))))
        off += clen
    end

    # Check length
    @assert off == sac_header_len
    check_npts && npts_in_file < trace.npts &&
        error("Number of points is not as expected: have $npts_in_file versus npts = " *
            "$(trace.npts) in header" * (file!="" ? " for file '$file'." : "."))

    # Now read in the trace
    trace.t .= reinterpret(SACFloat, data[(sac_header_len+1):end])
    native || (trace.t .= bswap.(trace.t))
    update_headers!(trace)
    any(isundefined.([trace.gcarc, trace.az, trace.baz])) && update_great_circle!(trace)
    trace
end

function SACtr(v::AbstractVector, delta::Real, b=0.0)
    s = SACtr(delta, length(v), b)
    s.t .= v
    update_headers!(s)
    s
end

"""
    getindex(A::AbstractArray{SACtr}, s::Symbol) -> Array{typeof(A[:].s)}
    A[:s] -> Array{typeof(A[:].s)}

Return an array of values containing the header with name `s` for the SACtr
traces.  This allows one to get all the headers values by doing A[:kstnm],
for example.
"""
getindex(A::AbstractArray{SACtr}, s::Symbol) = Array{typeof(getfield(A[1], s))}([getfield(a, s) for a in A])
getindex(t::SACtr, s::Symbol) = getfield(t, s) # Also define for single trace for consistency

"""
    setindex!(s::SACtr, value, s::Symbol)
    s[s] = value

Set the header with name `s` (a `Symbol`) to `value`.  E.g.:

    A[:kevnm] = "Blast 1"

    setindex!(A::AbstractArray{SACtr}, value, s::Symbol)
    A[s] = value

Set the header with name `s` for all the SACtr traces in the array `A`.  This
allows one to set all the headers at once for a set of traces by doing e.g.:

    A[:user0] = 1:length(A)

#### Setting trace values

When setting the values of the trace, use `s[:t] = k`.  This will automatically
update headers to reflect the trace and use broadcasting internally so that
the trace cannot be set to the wrong length.  Note however that a
broadcast assignment (`s[:t] .= k`) will not update the headers; this can be
done manually with `SAC.update_headers!(s)`.

    s[:t] = 1:s[:npts] # Sets to trace to be a line with constant slope
    s[:t] = sin.(2π*time(s)) # Sets the trace to a sine function
"""
function setindex!(t::SACtr, v, s::Symbol)
    if s == :t
        t.t .= v
        update_headers!(t)
    else
        setfield!(t, s, convert(typeof(getfield(t, s)), v))
        s in (:evlo, :evla, :stlo, :stla) && update_great_circle!(t)
    end
    t
end
setindex!(A::AbstractArray{SACtr}, V, s::Symbol) = setindex!.(A, V, s)

"""
    (==)(a::SACtr, b::SACtr) -> ::Bool

Return `true` if the traces `a` and `b` are equal (that is, have all fields the same),
and `false` otherwise.
"""
function (==)(a::SACtr, b::SACtr)
    for f in fieldnames(SACtr)
        if getfield(a, f) != getfield(b, f)
            return false
        end
    end
    true
end

"""
    copy(s::SACtr) -> t::SACtr

Return a copy of SAC trace `s`.
"""
function copy(s::SACtr)
    return SACtr(s.delta, s.depmin, s.depmax, s.scale, s.odelta, s.b, s.e, s.o, s.a, s.internal0,
        s.t0, s.t1, s.t2, s.t3, s.t4, s.t5, s.t6, s.t7, s.t8, s.t9, s.f,
        s.resp0, s.resp1, s.resp2, s.resp3, s.resp4, s.resp5, s.resp6, s.resp7, s.resp8, s.resp9,
        s.stla, s.stlo, s.stel, s.stdp, s.evla, s.evlo, s.evel, s.evdp, s.mag,
        s.user0, s.user1, s.user2, s.user3, s.user4, s.user5, s.user6, s.user7, s.user8, s.user9,
        s.dist, s.az, s.baz, s.gcarc, s.internal1, s.internal2, s.depmen, s.cmpaz, s.cmpinc,
        s.xminimum, s.xmaximum, s.yminimum, s.ymaximum,
        s.unused1, s.unused2, s.unused3, s.unused4, s.unused5, s.unused6, s.unused7,
        s.nzyear, s.nzjday, s.nzhour, s.nzmin, s.nzsec, s.nzmsec,
        s.nvhdr, s.norid, s.nevid, s.npts, s.internal3, s.nwfid, s.nxsize, s.nysize, s.unused8,
        s.iftype, s.idep, s.iztype, s.unused9, s.iinst, s.istreg, s.ievreg, s.ievtyp, s.iqual,
        s.isynth, s.imagtyp, s.imagsrc, s.unused10, s.unused11, s.unused12, s.unused13,
        s.unused14, s.unused15, s.unused16, s.unused17,
        s.leven, s.lpspol, s.lovrok, s.lcalda, s.unused18,
        s.kstnm, s.kevnm, s.khole, s.ko, s.ka, s.kt0, s.kt1, s.kt2, s.kt3, s.kt4, s.kt5, s.kt6, s.kt7,
        s.kt8, s.kt9, s.kf, s.kuser0, s.kuser1, s.kuser2, s.kcmpnm, s.knetwk, s.kdatrd, s.kinst,
        Base.copy(s.t))
end

@static if VERSION >= v"0.7-"
    Base.broadcastable(x::SACtr) = Ref(x)
end


text:
function sample_lines(shape::Tuple{Int64},redFac;sampleFunc="random",kargs...)
  error("Not implemented")
end

"""
    sample_lines(shape::Tuple{Int64,Int64},redFac::Float64;sampleFunc="random",kargs...)

generates a pattern to sample complete lines of an Array of size `shape` with a subsampling factor `redFac`.

# Arguments
* `shape::NTuple{N,Int64}` - size of the Array to be sampled
* `redFac::Float64`        - subsampling factor
* `sampleFunc="random"`        - name of the sampling function
                            ("random, "regular", "lines", "poisson" or "vdPoisson")
* `kargs...`               - addional keyword arguments
"""
function sample_lines(shape::Tuple{Int64,Int64},redFac::Float64;sampleFunc="random",kargs...)
  M,N=shape
  A = zeros(Int64,shape)
  yInd = sample((N,),redFac,sampleFunc;kargs...)
  A[1:M,yInd] .= 1;

  return (LinearIndices(A))[findall(x->x!=0, A)]
end

"""
    sample_lines(shape::Tuple{Int64,Int64,Int64},redFac::Float64;sampleFunc="random",kargs...)

generates a pattern to sample complete lines of an Array of size `shape` with a subsampling factor `redFac`.
The arguments are the same as in the 2d case with the exception that shape is of type `NTuple{3,Int64}`
"""
function sample_lines(shape::Tuple{Int64,Int64,Int64},redFac::Float64;sampleFunc="random",kargs...)
  M,N,Z=shape
  A2D = zeros(Int64,N,Z)
  A3D = zeros(Int64,shape)

  idx = sample(size(A2D),redFac,sampleFunc;kargs...)

  A2D[idx] .= 1;

  for i=1:M
    A3D[i,:,:] = A2D[:]
  end

  return (LinearIndices(A3D))[findall(x->x!=0, A3D)]
end


text:
@testset "Unwrap 1D" begin
    @test unwrap([0.1, 0.2, 0.3, 0.4]) ≈ [0.1, 0.2, 0.3, 0.4]
    @test unwrap([0.1, 0.2 + 2pi, 0.3, 0.4]) ≈ [0.1, 0.2, 0.3, 0.4]
    @test unwrap([0.1, 0.2 - 2pi, 0.3, 0.4]) ≈ [0.1, 0.2, 0.3, 0.4]
    @test unwrap([0.1, 0.2 - 2pi, 0.3 - 2pi, 0.4]) ≈ [0.1, 0.2, 0.3, 0.4]
    @test unwrap([0.1 + 2pi, 0.2, 0.3, 0.4]) ≈ [0.1, 0.2, 0.3, 0.4]
    @test unwrap([0.1, 0.2 + 2pi, 0.3, 0.4]) ≈ [0.1, 0.2, 0.3, 0.4]

    test_v = [0.1, 0.2, 0.3 + 2pi, 0.4]
    res_v = unwrap(test_v)
    @test test_v ≈ [0.1, 0.2, 0.3 + 2pi, 0.4]
    res_v .= 0
    unwrap!(test_v)
    @test test_v ≈ [0.1, 0.2, 0.3, 0.4]

    # test unwrapping within multi-dimensional array
    wrapped = [0.1, 0.2 + 2pi, 0.3, 0.4]
    unwrapped = [0.1, 0.2, 0.3, 0.4]
    wrapped = hcat(wrapped, wrapped)
    unwrapped = hcat(unwrapped, unwrapped)

    # test generically typed unwrapping
    types = (Float32, Float64)
    for T in types
        A_unwrapped = collect(range(0, stop=4convert(T, π), length=10))
        A_wrapped = A_unwrapped .% (2convert(T, π))

        test(I, J) = I ≈ J || I .+ 2π ≈ J || I .+ 4π ≈ J

        @test test(unwrap(A_wrapped), A_unwrapped)
        unwrap!(A_wrapped)
        @test test(A_wrapped, A_unwrapped)
    end
end

# tests for multi-dimensional unwrapping
@testset "Unwrap 2D" begin
    types = (Float32, Float64)
    for T in types
        v_unwrapped = collect(range(0, stop=4convert(T, π), length=7))
        A_unwrapped = v_unwrapped .+ v_unwrapped'
        A_wrapped = A_unwrapped .% (2convert(T, π))

        test_unwrapped = unwrap(A_wrapped)
        d = first(A_unwrapped) - first(test_unwrapped)
        @test (test_unwrapped .+ d) ≈ A_unwrapped
        unwrap!(A_wrapped)
        d = first(A_unwrapped) - first(A_wrapped)
        @test (A_wrapped .+ d) ≈ A_unwrapped

        v_unwrapped_range = collect(range(0, stop=4, length=7))
        A_unwrapped_range = v_unwrapped_range .+ v_unwrapped_range'
        test_range = convert(T, 2)
        A_wrapped_range = A_unwrapped_range .% test_range

    end
end

@testset "Unwrap 3D" begin
    types = (Float32, Float64)
    f(x, y, z) = 0.1x^2 - 2y + 2z
    f_wraparound2(x, y, z) = 5*sin(x) + 2*cos(y) + z
    f_wraparound3(x, y, z) = 5*sin(x) + 2*cos(y) - 4*cos(z)
    for T in types
        grid = range(zero(T), stop=2convert(T, π), length=11)
        f_uw = f.(grid, grid', reshape(grid, 1, 1, :))
        f_wr = f_uw .% (2convert(T, π))
        uw_test = unwrap(f_wr)
        offset = first(f_uw) - first(uw_test)
        @test (uw_test.+offset) ≈ f_uw rtol=eps(T) #oop, nowrap
        # test in-place version
        unwrap!(f_wr)
        offset = first(f_uw) - first(f_wr)
        @test (f_wr.+offset) ≈ f_uw rtol=eps(T) #ip, nowrap

        f_uw = f_wraparound2.(grid, grid', reshape(grid, 1, 1, :))
        f_wr = f_uw .% (2convert(T, π))
        uw_test = unwrap(f_wr)
        offset = first(f_uw) - first(uw_test)
        @test (uw_test.+offset) ≈ f_uw #oop, 2wrap
        # test in-place version
        unwrap!(f_wr)
        offset = first(f_uw) - first(f_wr)
        @test (f_wr.+offset) ≈ f_uw #ip, 2wrap

        f_uw = f_wraparound3.(grid, grid', reshape(grid, 1, 1, :))
        f_wr = f_uw .% (2convert(T, π))
        uw_test = unwrap(f_wr)
        offset = first(f_uw) - first(uw_test)
        @test (uw_test.+offset) ≈ f_uw #oop, 3wrap
        # test in-place version
        unwrap!(f_wr)
        offset = first(f_uw) - first(f_wr)
        @test (f_wr.+offset) ≈ f_uw #oop, 3wrap
    end
end


text:
using DataFrames
using Plots
using CSV
using StatsPlots

const FULL_FILE_PROCESSED = joinpath("data","Melbourne_housing_FULL_processed_2.csv")

# Initialise our dataframe
df_full = DataFrame(CSV.File(FULL_FILE_PROCESSED))

"""
    prepare_cumplot(df_full::AbstractDataFrame, column)

Returns a dataframe with two columns, `column` and `cum_rownum`. All `missing` values (in `column` only) are dropped, and the existing values are sorted. `cum_rownum` holds each row's index in reverse - a single variable cumulative sum.
"""
function prepare_cumplot(df_full::AbstractDataFrame, column)
    df_cum = select(df_full, column) |> dropmissing! |> sort!
    df_cum[!, :cum_rownum] = nrow(df_cum):-1:1
    return df_cum
end

###
# Rooms
###

# No missing data, so we can proceed as always

"""
    plot_rooms(df_full::AbstractDataFrame)

Returns Room distribution as a bar graph.

X = num. of rooms, Y = num. of properties that contain that many rooms; Y is in log10 scale.
"""
function plot_rooms(df_full::AbstractDataFrame)
    df_rooms = combine(groupby(df_full, :Rooms), nrow => :count)

    # Discrete data? Bar graph.
    return @df df_rooms plot(:Rooms, :count, seriestype = :bar, xlabel = "Rooms", ylabel = "Property Count", label = "", yscale=:log10, xticks=1:20, plot_title = "Room Distribution")
end

###
# Price
###

# We are operating on a single variable, so we can simply discard missing values
"""
    plot_price(df_full::AbstractDataFrame)

Plots distribution of price as a cumulative sum line.

X = price (\$AUD), Y = number of propetries that have sold for more than that price; Y is in log10 scale.
"""
function plot_price(df_full::AbstractDataFrame)
    # df_price = select(df_full, :Price) |> dropmissing! |> sort!
    # df_price[!, :cum_rownum] = 1:-1:nrow(df_price)
    df_price = prepare_cumplot(df_full, :Price)
    transform!(df_price, :Price => cumsum => :Price_cum)

    # We will use a cumulative plot. It feels appropriate, considering we are using continuous data of a large volume.
    return @df df_price plot(:Price, :cum_rownum, xscale=:log10, xticks = ([100000, 1000000, 10000000], ["\$10,000", "\$100,000", "\$1,000,000"]), yformatter = :plain, label="", xlabel = "Property Sale Price (\$AUD)", ylabel = "Num. of Properties Sold For More", plot_title = "Cumulative Sale Price")#, yticks = 1:10^8:5*10^8)
end

###
# Distance
###

# One special case 'bad value', handled by our cleanup
"""
    plot_distance(df_full::AbstractDataFrame)

Plots distribution of price as a distance sum line.

X = distance (km), Y = number of propetries that are further from the CBD than that distance; Y is in linear scale.
"""
function plot_distance(df_full::AbstractDataFrame)
    df_distance = prepare_cumplot(transform(df_full |> dropmissing, :Distance => ByRow(x->x isa String ? parse(Float16, x) : x) => :Distance), :Distance)

    # We will use a cumulative plot. It feels appropriate, considering we are using continuous data of a large volume.
    return @df df_distance plot(:Distance, :cum_rownum, xticks=0:10:50, yformatter = :plain, label="", xlabel="Distance From CBD (km)", ylabel="Num. of Properties Further From CBD", plot_title="Cumulative CBD Distance")
end

"""
    plot_method(df_full::AbstractDataFrame)

Plots distribution of sale method as a bar chart.

X = sale method (abbreviated), Y = num. of prop. that have been sold via that method; Y is in linear scale.
"""
function plot_method(df_full::AbstractDataFrame)
    df_method = sort!(combine(groupby(df_full, :Method), nrow => :count), :count)

    # Discrete data? Bar graph.
    return @df df_method plot(:Method, :count, seriestype = :bar, xlabel = "Sale Method", ylabel = "Property Count", label = "", yformatter = :plain, plot_title = "Sale Method")
end

"""
    plot_landsize(df_full::AbstractDataFrame)

Plots distribution of land size method as a line graph. Landsizes have been grouped

X = landsize (m²), Y = num. of prop. sold with that landsize; X is in log10 scale.
"""
function plot_landsize(df_full::AbstractDataFrame)
    # Manipulating the data here for convenience
    df_landsize = combine(groupby(transform!(filter!(:Landsize => x->x>5, select(df_full, :Landsize) |> dropmissing! |> sort!), :Landsize => ByRow(x -> convert(Int, x - 5 |> float |> x->(10*log10(x)) |> floor)) => :Landsize_Grouped), :Landsize_Grouped), nrow => :count)

    # Function to convert axes from log10 to what they should be
    adjustment_func(num::Int) = convert(Int, 10^(num / 10.0) |> floor)

    # X-Axis won't be 100% accurate, be difference would be too minute to matter
    return @df df_landsize plot(:Landsize_Grouped, :count, seriestype = :line, xticks=(10:10:56, adjustment_func.(10:10:56)), label="", tickfonthalign=:center, xlabel="Landsize (m²)", ylabel="Number of Properties", plot_title="Property Landsizes")
end

"""
    generate_task1_df(filename=FULL_FILE_PROCESSED)

Generates a `DataFrame` from a CSV file for use throughout this task.
"""
generate_task1_df(filename=FULL_FILE_PROCESSED) = DataFrame(CSV.File(filename))

"""
    save_plots(tasknum::Int, df_creator::Function;plot_funcs::Vector{Function}, filename=FULL_FILE_PROCESSED)

Given a task number, a `DataFrame` generator `Function` that takes a file name as a parameter, and a `Vector` of `Functions` which take this `DataFrame` and output `Plot`s, this function saves all such `Plot`s as appropriately named .png files in the output directory.
"""
function save_plots(tasknum::Int, df_creator::Function;plot_funcs::Vector{Function}, filename=FULL_FILE_PROCESSED)
    df = df_creator(filename)
    for (i, func) in enumerate(plot_funcs)
        savefig(func(df), "output/task$(tasknum).$(i).png")
    end
end;


text:
struct ModifiedState{S<:State, C<:AbstractStateChanges} <: FiniteBasisState
    state::S
    changes::C
end

ModifiedState(state::State) = ModifiedState(state, StateChanges(state))

# accessors
raw_state(s::ModifiedState) = s.state
changes(s::ModifiedState)   = s.changes

# Modifiers
apply_changes!(s::ModifiedState) = begin
    apply!(raw_state(s), changes(s))
    zero!(changes(s))
    return s
end
apply_warn_raw!(s::ModifiedState) = begin
    #@warn "applying changes not optimized"
    apply_changes!(s)
    return raw_state(s)
end
clear_changes!(s::ModifiedState) = zero!(changes(s))

## For DoubleState
apply_chages!(s::DoubleState) = nothing
function apply_chages!(s::DoubleState{<:ModifiedState})
    apply_changes!(row(s))
    apply_changes!(col(s))
    return s
end
clear_changes!(s::DoubleState) = nothing
function clear_changes!(s::DoubleState{<:ModifiedState})
    clear_changes!(row(s))
    clear_changes!(col(s))
    return s
end

raw_config(s::DoubleState) = config(s)
raw_config(s::DoubleState{<:ModifiedState}) = (raw_config(s.σ_row), raw_config(s.σ_col))
## end

## State Interface : custom accessors
@inline config(s::ModifiedState) = config(apply_warn_raw!(s))
@inline raw_config(s::ModifiedState) = config(raw_state(s))

## State Interface : Property Accesors
@inline spacedimension(s::ModifiedState) = spacedimension(raw_state(s))
@inline nsites(s::ModifiedState) = nsites(raw_state(s))
@inline local_dimension(s::ModifiedState{S,C}) where {S,C} = local_dimension(S)
@inline local_dimension(s::Type{ModifiedState{S,C}}) where {S,C} = local_dimension(S)
@inline eltype(s::ModifiedState) = eltype(raw_state(s))

toint(s::ModifiedState) = toint(apply_warn_raw!(s))
index(s::ModifiedState) = index(apply_warn_raw!(s))
index_to_int(s::ModifiedState) = index_to_int(apply_warn_raw!(s))

## State Interface : Checks
same_basis(s1::FiniteBasisState, s2::FiniteBasisState) =
    nsites(s1)==nsites(s2) && local_dimension(s1) == local_dimension(s2) && eltype(s1) == eltype(s2)

## State Interface : Operations on tuples
_toint(left::ModifiedState, right::ModifiedState) = _toint(apply_warn_raw!(left), apply_warn_raw!(right))

## State Interface : flip operations
# This is the standard method, that applies all exhisting transformations
function flipat!(rng::AbstractRNG, state::ModifiedState, i::Int)
    old_val = config(state)[i]
    return _flipat!(rng, state, i, old_val)
end

# This is a faster method, that assumes that you are flipping spins that
# are not flipped in the changes
function flipat_fast!(rng::AbstractRNG, state::ModifiedState, i::Int)
    old_val = config(raw_state(state))[i]
    return _flipat!(rng, state, i, old_val)
end

# Inner method actually performing the flip
function _flipat!(rng::AbstractRNG, state::ModifiedState, i::Int, old_val)
    # Generate new value
    N = local_dimension(raw_state(state))
    T = eltype(raw_state(state))
    new_val = T(rand(rng, 0:(N-1)))
    while new_val == old_val
        new_val = T(rand(rng, 0:(N-1)))
    end

    # Modify
    push!(changes(state), (i, new_val))
    return (old_val, new_val)
end

function setat!(state::ModifiedState, i::Int, val)
    old_val = config(raw_state(state))[i]
    push!(changes(state), (i, val))

    return old_val
end

function setat_fast!(state::ModifiedState, i::Int, val)
    old_val = config(raw_state(state))[i]
    push!(changes(state), (i, val))

    return old_val
end

set_index!(s::ModifiedState, val) = (zero!(changes(s)); set_index!(raw_state(s), val))
set!(s::ModifiedState, val) = (zero!(changes(s)); set!(raw_state(s), val); return s)
add!(s::ModifiedState, val) = add!(apply_warn_raw!(state), val)
rand!(rng, s::ModifiedState) = (zero!(changes(s));
    rand!(rng, raw_state(s)))

Base.show(io::IO, ::MIME"text/plain", s::ModifiedState) =
    print(io, "ModifiedState: $(raw_state(s)) with changes $(changes(s))")


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENCE in the project root.
# ------------------------------------------------------------------

ST.scitype(::Composition{D}, ::ST.DefaultConvention) where {D} = ST.Compositional{D}

function ST.scitype(::CoDaArray{D}, ::ST.DefaultConvention) where {D}
  T = ST.Compositional{D}
  AbstractVector{T}
end


text:
using Distributions;
using NLsolve;
using ForwardDiff;

include(joinpath(@__DIR__, "dipole_response.jl"));

ϕ_dist = Uniform(0.0, 2*π);
θ_dist = Uniform(0.0, π);

# the way that the energy functions are organized is a mess; I'm sorry
# also... not sorry
abstract type Energy end

mutable struct EAPChain
  b::Float64;
  E0::Float64;
  μ::DipoleResponse;
  UFunction::Energy;
  kT::Float64;
  Fz::Float64;
  Fx::Float64;
  ϕs::FdVector;
  cϕs::FdVector;
  sϕs::FdVector;
  θs::FdVector;
  cθs::FdVector;
  sθs::FdVector;
  n̂s::FdMatrix;
  Ω::Float64;
  μs::FdMatrix;
  us::FdVector;
  xs::FdMatrix;
  r::FdVector;
  U::Float64;
end

@inline n(chain::EAPChain) = length(chain.ϕs);

@inline n̂(cϕ, sϕ, cθ, sθ) = [cϕ * sθ; sϕ * sθ; cθ];

@inline n̂j(chain::EAPChain, idx::Int) = n̂(chain.cϕs[idx], chain.sϕs[idx],
                                          chain.cθs[idx], chain.sθs[idx]);

function update_xs!(chain::EAPChain, idx::Int)
  @warn "This implementation of update_xs! is generally slower than the vectorized implementation";
  if idx == 1
    @inbounds chain.xs[:, idx] = chain.b / 2 * n̂j(chain, idx);
    idx += 1;
  end
  for i=idx:n(chain)
    @inbounds chain.xs[:, i] = (chain.xs[:, i-1] +
                                chain.b / 2 * (n̂j(chain, i) + n̂j(chain, i-1)));
  end
end

@inline function update_xs!(chain::EAPChain)
  chain.xs[:, :] = chain.b*(cumsum(chain.n̂s, dims=2) - 0.5*chain.n̂s);
end

#=
function update_xs!(chain::EAPChain)
  @show chain.n̂s
  @show cumsum(chain.n̂s, dims=2);
  @show (cumsum(chain.n̂s, dims=2) - 0.5*chain.n̂s);
  idx = rand(1:n(chain));
  @assert(dot(chain.n̂s[:, idx], chain.n̂s[:, idx]) ≈ 1.0);
  chain.xs[:, :] = chain.b*(cumsum(chain.n̂s, dims=2) - 0.5*chain.n̂s);
end
=#

@inline u(E0::Real, μ::FdVector) = -1/2*E0*μ[3];

function EAPChain(pargs::Dict)
  ϕs = rand(ϕ_dist, pargs["num-monomers"]);
  θs = rand(θ_dist, pargs["num-monomers"]);

  dr = if pargs["chain-type"] == "dielectric"
    DielectricResponse(pargs["K1"], pargs["K2"]);
  elseif pargs["chain-type"] == "polar"
    PolarResponse(pargs["mu"] * [1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 1.0]);
  else
    error("chain-type is not understood.");
  end

  ret =  EAPChain(
                  pargs["mlen"],
                  pargs["E0"],
                  dr,
                  if pargs["energy-type"] == "noninteracting"
                    NonInteractingEnergy();
                  elseif pargs["energy-type"] == "interacting"
                    InteractingEnergy();
                  elseif pargs["energy-type"] == "Ising"
                    IsingEnergy();
                  else
                    error("energy-type is not understood.");
                  end,
                  pargs["kT"],
                  pargs["Fz"],
                  pargs["Fx"],
                  ϕs,
                  map(cos, ϕs),
                  map(sin, ϕs),
                  θs,
                  map(cos, θs),
                  map(sin, θs),
                  zeros(3, pargs["num-monomers"]),
                  log(prod(map(sin, θs))), # store logarithm of solid angle instead
                  zeros(3, pargs["num-monomers"]),
                  zeros(pargs["num-monomers"]),
                  zeros(3, pargs["num-monomers"]),
                  zeros(3),
                  0.0
                 );
  ret.n̂s[:, :] = hcat(map(j -> n̂j(ret, j), 1:n(ret))...);
  @simd for i=1:n(ret)
    @inbounds ret.μs[:, i] = dr(ret.E0, ret.cϕs[i], ret.sϕs[i],
                                ret.cθs[i], ret.sθs[i]);
  end
  ret.us[:] = map(i -> u(ret.E0, view(ret.μs, :, i)), 1:n(ret));
  update_xs!(ret);
  ret.r[:] = end_to_end(ret);
  ret.U = U(ret);
  return ret;
end

function EAPChain(chain::EAPChain)
  return EAPChain(
                  chain.b,
                  chain.E0,
                  chain.μ,
                  chain.UFunction,
                  chain.kT,
                  chain.Fz,
                  chain.Fx,
                  chain.ϕs[:],
                  chain.cϕs[:],
                  chain.sϕs[:],
                  chain.θs[:],
                  chain.cθs[:],
                  chain.sθs[:],
                  chain.n̂s[:, :],
                  chain.Ω,
                  chain.μs[:, :],
                  chain.us[:],
                  chain.xs[:, :],
                  chain.r[:],
                  chain.U
                 );
end

# TODO: consider rewriting this to make better use of past calculations
# this is probably by far the slowest calculation *shrug emoji*
function U_interaction(chain::EAPChain)
  U = 0.0;
  @inbounds for i=1:n(chain)
    @simd for j=i+1:n(chain) # once debugged, use inbounds
      r = chain.xs[:, i] - chain.xs[:, j];
      r2 = dot(r, r);
      rmag = sqrt(r2);
      r̂ = r / rmag;
      r3 = r2*rmag;
      μi = view(chain.μs, :, i);
      μj = view(chain.μs, :, j);
      U += (dot(μi, μj) - 3*dot(μi, r̂)*dot(μj, r̂)) / (4*π*r3);
    end
  end
  return U;
end

# TODO: consider rewriting this to make better use of past calculations
# this is probably by far the slowest calculation *shrug emoji*
function U_Ising(chain::EAPChain)
  U = 0.0;
  @inbounds for i=1:n(chain)-1
    r = chain.xs[:, i] - chain.xs[:, i+1];
    r2 = dot(r, r);
    rmag = sqrt(r2);
    r̂ = r / rmag;
    r3 = r2*rmag;
    μi = view(chain.μs, :, i);
    μj = view(chain.μs, :, i+1);
    U += (dot(μi, μj) - 3*dot(μi, r̂)*dot(μj, r̂)) / (4*π*r3);
  end
  return U;
end

function move!(chain::EAPChain, idx::Int, dϕ::Real, dθ::Real)
  @inbounds begin
    chain.ϕs[idx] += dϕ;
    chain.cϕs[idx] = cos(chain.ϕs[idx]);
    chain.sϕs[idx] = sin(chain.ϕs[idx]);

    chain.θs[idx] = min(π, max(0.0, chain.θs[idx]+dθ));
    sθ = sin(chain.θs[idx]);
    chain.Ω += log(sθ / chain.sθs[idx]); # update solid angle
    chain.cθs[idx] = cos(chain.θs[idx]);
    chain.sθs[idx] = sθ;

    # the order of these updates is important
    chain.n̂s[:, idx] = n̂j(chain, idx);
    chain.μs[:, idx] = chain.μ(chain.E0, chain.cϕs[idx], chain.sϕs[idx],
                               chain.cθs[idx], chain.sθs[idx]);
    chain.us[idx] = u(chain.E0, view(chain.μs, :, idx));
    update_xs!(chain);
    chain.r[:] = end_to_end(chain);
    chain.U = U(chain);
  end
  return true;
end

function flip_n!(chain::EAPChain, idx::Int)
  move!(chain, idx, π, π - 2*chain.θs[idx])
end

pflip_linear(ip::Real) = (1 + ip) / 2;

function cluster_flip!(chain::EAPChain, idx::Int;
                       pflip::Function = pflip_linear,
                       ϵflip::Real = 1.0,
                       ηerr::Real = 1e-10
                      )

  # grow the cluster to the right
  upper_p = NaN;
  upper_idx = idx;
  while true
    if upper_idx >= n(chain)
      upper_p = 0.0;
      break;
    end
    upper_p = pflip_linear(dot(n̂j(chain, upper_idx), n̂j(chain, upper_idx+1)));
    if rand() <= upper_p
      upper_idx += 1
    else
      break;
    end
  end
  @assert(!isnan(upper_p) && upper_idx <= n(chain));

  # grow the cluster to the left
  lower_p = NaN;
  lower_idx = idx;
  while true
    if lower_idx <= 1
      lower_p = 0.0;
      break;
    end
    lower_p = pflip_linear(dot(n̂j(chain, lower_idx), n̂j(chain, lower_idx-1)));
    if rand() <= lower_p
      lower_idx -= 1
    else
      break;
    end
  end
  @assert(!isnan(lower_p) && lower_idx >= 1);

  if rand() <= ϵflip
    prev_n̂s = copy(chain.n̂s[:, lower_idx:upper_idx]);
    for i in lower_idx:upper_idx
      flip_n!(chain, i);
    end

    if norm(prev_n̂s + chain.n̂s[:, lower_idx:upper_idx], Inf) > ηerr
      println("prev");
      display(prev_n̂s)
      println();
      println("curr");
      display(chain.n̂s[:, lower_idx:upper_idx])
      println();
      @show(norm(prev_n̂s + chain.n̂s[:, lower_idx:upper_idx], Inf))
    end
    @assert(norm(prev_n̂s + chain.n̂s[:, lower_idx:upper_idx], Inf) < ηerr);

    new_upper_p = if upper_idx < n(chain)
      pflip_linear(dot(n̂j(chain, upper_idx), n̂j(chain, upper_idx+1)));
    else
      0
    end
    new_lower_p = if lower_idx > 1
      pflip_linear(dot(n̂j(chain, lower_idx), n̂j(chain, lower_idx-1)));
    else
      0
    end
    return (
             ((1 - new_upper_p)*(1 - new_lower_p)) /
             ((1 - upper_p)*(1 - lower_p))
            );
  else
    return 1.0;
  end

end

function move!(chain::EAPChain, idx::Int, dϕ::Real, dθ::Real, r0::AbstractVector;
               frac_mv::Real = 0.15, acc_tol::Real = 1e-1)

  @assert(frac_mv > 0.0 && frac_mv <= 1.0,
          "Fraction of chain to move must be between 0.0 and 1.0");

  L0 = chain.b * n(chain);

  trial_idxs = if (idx / n(chain)) >= (1 - frac_mv)
    max(1, round(Int, (1 - frac_mv)*n(chain))):n(chain);
  elseif (idx / n(chain)) <= frac_mv
    1:(min(n(chain), round(Int, frac_mv*n(chain))));
  else
    (max(1, round(Int, idx - frac_mv*n(chain)/2.0))):round(Int, idx + frac_mv*n(chain)/2.0)
  end

  r1 = chain.b * ([
                   sum(chain.cϕs[trial_idxs] .* chain.sθs[trial_idxs]);
                   sum(chain.sϕs[trial_idxs] .* chain.sθs[trial_idxs]);
                   sum(chain.cθs[trial_idxs])
                  ]);

  # internal lambda function for solving nonlinear system
  f! = (F::Vector, x::Vector) -> begin;
    n = convert(Int, length(x) / 2);
    cϕs = map(cos, x[1:n]);
    sϕs = map(sin, x[1:n]);
    sθs = map(sin, x[(n+1):(2*n)]);
    F[1] = chain.b * sum(cϕs .* sθs) - r1[1];
    F[2] = chain.b * sum(sϕs .* sθs) - r1[2];
    F[3] = chain.b * sum(map(cos, x[(n+1):(2*n)])) - r1[3];
  end

  #=
  j! = (J::Matrix, x::Vector) -> begin;
    n = convert(Int, length(x) / 2);
    cϕs = map(cos, x[1:n]);
    sϕs = map(sin, x[1:n]);
    cθs = map(cos, x[(n+1):(2*n)]);
    sθs = map(sin, x[(n+1):(2*n)]);
    J[:, :] = chain.b * hcat(map(i -> [cϕs[i]*cθs[i] -sϕs[i]*sθs[i];
                                       sϕs[i]*cθs[i]  cϕs[i]*sθs[i];
                                       -sθs[i]        0.0;], 1:n)...);
  end
  =#

  # trial move initial guess
  x0 = vcat(chain.ϕs[trial_idxs], chain.θs[trial_idxs]);
  local_idx = findnext(i -> i==idx, collect(trial_idxs), 1);
  x0[local_idx] += dϕ;
  x0[2*local_idx] += dθ;

  # solve nonlinear system to ensure end to end vector constrain is satisfied
  result = nlsolve(f!, x0);
  if result.residual_norm / L0 > acc_tol
    @warn "residual norm is high", result.residual_norm / L0;
    return false;
  end

  # update chain based on result
  m = length(trial_idxs);
  chain.ϕs[trial_idxs] = result.zero[1:m];
  chain.θs[trial_idxs] = result.zero[(m+1):(2*m)];
  chain.cϕs[trial_idxs] = map(cos, chain.ϕs[trial_idxs]);
  chain.sϕs[trial_idxs] = map(sin, chain.ϕs[trial_idxs]);
  chain.cθs[trial_idxs] = map(cos, chain.θs[trial_idxs]);
  chain.sθs[trial_idxs] = map(sin, chain.θs[trial_idxs]);

  # the order of these updates is important
  chain.n̂s[:, trial_idxs] = hcat(map(idx -> n̂j(chain, idx), trial_idxs)...);
  chain.μs[:, trial_idxs] = hcat(map(idx -> chain.μ(chain.E0, chain.cϕs[idx],
                                                    chain.sϕs[idx], chain.cθs[idx],
                                                    chain.sθs[idx]),
                                     trial_idxs)...);
  chain.us[trial_idxs] = map(idx -> u(chain.E0, view(chain.μs, :, idx)),
                             trial_idxs);
  update_xs!(chain);
  chain.r[:] = end_to_end(chain);
  chain.U = U(chain);

  return true;
end

@inline end_to_end(chain::EAPChain) = @inbounds (chain.xs[:, end] +
                                                 chain.b/2.0*n̂j(chain, n(chain)));

@inline chain_μ(chain::EAPChain) = reshape(sum(chain.μs, dims=2), 3);

include(joinpath(@__DIR__, "energy.jl"));
@inline U(chain::EAPChain) = chain.UFunction(chain); # forward energy call to chain

max_stars_repo_name: grasingerm/polymer-stats

max_issues_repo_name: grasingerm/polymer-stats

max_forks_repo_name: grasingerm/polymer-stats


text:
module Clustering

using LinearAlgebra
using Statistics

export GMM, Kmeans, Xmeans, DBSCAN, HDBSCAN

include("utils.jl")
include("GMM.jl")
include("Kmeans.jl")
include("DBSCAN.jl")
include("HDBSCAN.jl")
include("Xmeans.jl")

end


text:
module SetIntegerParameterTrigger

using ..Ahorn, Maple

@mapdef Trigger "ShaderHelper/SetIntegerParameterTrigger" SetIntegerParameterTrigger(x::Integer, y::Integer, width::Integer=Maple.defaultTriggerWidth, height::Integer=Maple.defaultTriggerHeight, shader::String="shaderhelper/testshader", value::Int=0, key::String="key")

const placements = Ahorn.PlacementDict(
    "Set Integer Parameter Trigger (Shader Helper)" => Ahorn.EntityPlacement(
        SetIntegerParameterTrigger,
        "rectangle"
    )
)

end

max_stars_repo_name: lordseanington/ShaderHelper

max_issues_repo_name: lordseanington/ShaderHelper

max_forks_repo_name: lordseanington/ShaderHelper


text:
function output(t, t_output)
  if t > t_output
    x_array = []
    v_array = []
    rho_array = []
    P_array = []
    outputfilename = @sprintf("%s_%03d.txt", outputfile, output_count)

    open(outputfilename, "w") do io
      println(io, "### x, v, rho, P")
      for i in 1:Npart
        x = sph[i].x
        v = sph[i].p / sph[i].m
        rho = sph[i].rho
        P = sph[i].P
        hsml = grad_h(i)
        push!(x_array, x)
        push!(v_array, v)
        push!(rho_array, rho)
        push!(P_array, P)
        println(io, "$x $v $rho $P $hsml")
      end
    end

    if plot_figure
      step_riemann = lbox / Nsample_riemann
      x_riemann = x_min:step_riemann:x_max
      rho_riemann_array = map(x -> rho_riemann(x, t, v_s, P_s, x_center), x_riemann)
      v_riemann_array = map(x -> v_riemann(x, t, v_s, P_s, x_center), x_riemann)
      P_riemann_array = map(x -> P_riemann(x, t, v_s, P_s, x_center), x_riemann)
      settings = "$formulation, $gradient gradient\n$time_integrator, $kernel, η=$eta_hsml\nvolume element: $volume_element\ntime-dependent art. visc.: $time_dependent_viscosity\n"

      plt_rho = Plots.plot(; ylabel="rho", title=settings, titlefont=font(12),
                           ylims=(rho_min, rho_max), xlims=(x_min_plot, x_max_plot))
      plt_rho = Plots.plot!(x_riemann, rho_riemann_array; label="exact")
      plt_rho = Plots.scatter!(x_array, rho_array; markersize=2, label="sim")

      plt_P = Plots.plot(; ylabel="P", xlabel="x", ylims=(P_min, P_max),
                         xlims=(x_min_plot, x_max_plot))
      plt_P = Plots.plot!(x_riemann, P_riemann_array; label="exact")
      plt_P = Plots.scatter!(x_array, P_array; markersize=2, label="sim")

      plt_v = Plots.plot(; ylabel="v", ylims=(v_min, v_max), xlims=(x_min_plot, x_max_plot))
      plt_v = Plots.plot!(x_riemann, v_riemann_array; label="exact")
      plt_v = Plots.scatter!(x_array, v_array; markersize=2, label="sim")

      Plots.plot(plt_rho, plt_v, plt_P; layout=(3, 1), size=(400, 600))
      figname = @sprintf("%s_%s_%s_%s_%s_%s_%s_fig_%03d.png", outputfile, formulation, gradient, time_integrator, kernel, volume_element, time_dependent_viscosity, output_count)
      Plots.savefig(figname)
    end
    global t_output += output_interval
    global output_count += 1
  end
end


text:
"""
    bson_read([T=Any], src::Union{IO, DenseVector{UInt8}})

Read an object from the BSON document encoded in `src`.
"""
bson_read(::Type{T}, bytes::DenseVector{UInt8}) where T = BSONReader(bytes)[T]

"""
    bson_read([T=Any], path::AbstractString)

Read an object from the BSON document at `path`.
"""
bson_read(::Type{T}, path::AbstractString) where T = bson_read(T, read(path))

function bson_read(::Type{T}, io::IO) where T
    buf = UInt8[]
    readbytes!(io, buf, typemax(Int))
    bson_read(T, buf)
end

bson_read(bytes::DenseVector{UInt8}) = bson_read(Any, bytes)
bson_read(path::AbstractString) = bson_read(Any, path)
bson_read(io::IO) = bson_read(Any, io)

"""
    bson_write(dst::Union{IO, DenseVector{UInt8}}, x)
    bson_write(dst::Union{IO, DenseVector{UInt8}}, xs::Pair...)

Encode `x` or each element of `xs` as a BSON document in `dst` and return `dst`.
"""
function bson_write(buf::DenseVector{UInt8}, x)
    writer = BSONWriter(buf)
    writer[] = x
    close(writer)
    buf
end

function bson_write(buf::DenseVector{UInt8}, xs::Pair...)
    writer = BSONWriter(buf)
    for (k, v) in xs
        writer[k] = v
    end
    close(writer)
    buf
end

function bson_write(io::IO, x...)
    buf = UInt8[]
    bson_write(buf, x...)
    write(io, buf)
    io
end

"""
    bson_write(path::AbstractString, x)
    bson_write(path::AbstractString, xs::Pair...)

Encode `x` or each element of `xs` as fields of a BSON document and write to `path`.
"""
function bson_write(path::AbstractString, x...)
    open(path, "w") do io
        bson_write(io, x...)
    end
    nothing
end


text:
using FFL
using Test

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


text:
#Copyright (c) 2015 Rex Computing and Isaac Yonemoto

#see LICENSE.txt

#this work was supported in part by DARPA Contract D15PC00135

#fraczero-speed.jl

#testing to see how fast we can test to see if the fraction is zero.

#method #1 - compare against the zero array.

const count = 1000000

function constfn(b, i)
  b && (i == zero(UInt64))
end

function forloopfn(b)
  for idx = 1:length(b)
    b[idx] != 0 && return false
  end
  return true
end

function forloopfn2(b)
  for idx = length(b):-1:1
    b[idx] != 0 && return false
  end
  return true
end

function testy()

#make a cell array which will hold four-unit Int64 arrays
arrays = cell(count)
for idx = 1:count
  arrays[idx] = [(rand() < 0.25) ? zero(UInt64) : rand(UInt64) for jdx = 1:4]
end

println("by comparing against a zero array.")
sum = 0
tic()
for (i=1:count)
  sum += (arrays[i] == zeros(UInt64, 4)) ? 1 : 0
end
toc()
#println forces julia to use the result and not compile it away.
println(sum)

println("by comparing against a cached zero array.")
cache = zeros(UInt64, 4)
sum = 0
tic()
for (i=1:count)
  sum += (arrays[i] == cache) ? 1 : 0
end
toc()
#println forces julia to use the result and not compile it away.
println(sum)

println("by mapping an anonymous function")
sum = 0
tic()
for (i=1:count)
  sum += (reduce((b, i) -> b && (i == zero(UInt64)), true, arrays[i])) ? 1 : 0
end
toc()
#println forces julia to use the result and not compile it away.
println(sum)

println("by mapping a cached anonymous function.")
cachefn = (b, i) -> b && (i == zero(UInt64))
sum = 0
tic()
for (i=1:count)
  sum += (reduce(cachefn, true, arrays[i])) ? 1 : 0
end
toc()
#println forces julia to use the result and not compile it away.
println(sum)

println("by mapping a const function.")
sum = 0
tic()
for (i=1:count)
  sum += (reduce(constfn, true, arrays[i])) ? 1 : 0
end
toc()
#println forces julia to use the result and not compile it away.
println(sum)

println("with a breakaway for loop.")
sum = 0
tic()
for (i=1:count)
  sum += (forloopfn(arrays[i])) ? 1 : 0
end
toc()
#println forces julia to use the result and not compile it away.
println(sum)

println("breakaway for loop, downwards.")
sum = 0
tic()
for (i=1:count)
  sum += (forloopfn2(arrays[i])) ? 1 : 0
end
toc()
#println forces julia to use the result and not compile it away.
println(sum)
end

testy()
testy()
testy()

#results:  a breakaway for loop was definitely the fastest.  It's not even a contest.

#by comparing against a zero array.
#elapsed time: 0.234025668 seconds
#by comparing against a cached zero array.
#elapsed time: 0.17063435 seconds
#by mapping an anonymous function
#elapsed time: 0.297591116 seconds
#by mapping a cached anonymous function.
#elapsed time: 0.300173192 seconds
#by mapping a const function.
#elapsed time: 0.290023684 seconds
#with a breakaway for loop.
#elapsed time: 0.031858592 seconds
#breakaway for loop, downwards.
#elapsed time: 0.037607769 seconds

max_stars_repo_name: REX-Computing/unumjl

max_issues_repo_name: REX-Computing/unumjl

max_forks_repo_name: REX-Computing/unumjl


text:
# Dependencies
using Distributions # not required per se
using Plots
using Logging
import Base: show, zero
# Constants
const states = ['U','N','A','F']
const cs = Dict('N'=>1/2, 'A'=>1/4, 'F'=>1/8)
# helper function for indices
inds(i,j) = CartesianIndex.([(i-1, j-1), (i,j-1), (i+1, j-1),
            (i-1, j),             (i+1, j),
            (i-1, j+1), (i,j+1), (i+1, j+1)])

"""
    Cell

Cell in our world. State should be one of 'U', 'N', 'A','F','B' for
respectively Urban, Nature, Agriculture, Forest and Border.
"""
struct Cell
    state::Char
    age::Int
end
zero(::Type{Cell}) = Cell('b',0)
show(io::IO, c::Cell) = print(io, "($(c.state), $(c.age))")

"""
    getneigbors(world::Array{T,2}, i::Int,j::Int)

return the eight neighbors of the world in position (i,j)
"""
function getneigbors(world::Array{T,2}, i::Int,j::Int) where T
    return view(world, inds(i,j))
end

"""
    trfprob(c::Cell)

Return transformation probability for a non-urban to urban cell.
"""
function trfprob(c::Cell)
    return tanh((c.age+1/2)/40)
end

"""
    urbanisation(c::Cell, n)

return the next states for a cell and its neighbors

see also: [getneigbors](@ref)
"""
function urbanisation(c::Cell, n)
    # urban stays urban
    if c.state == 'U'
        @debug "Town gets older"
        return Cell('U', c.age + 1)
    # nature evolves to urban
    else
        # within mutation treshold
        if rand() < sum(map(trfprob, filter(x-> x.state =='U',n))) / 8 * cs[c.state]
            return Cell('U', 0)
            @info "Added new city!!!"
        # otherwise just age
        else
            @debug "$(c.state) gets older"
            return Cell(c.state, c.age + 1)
        end
    end
end

"""
    ruralisation(c::Cell, n)

return the next states for a cell and its neighbors

see also: [getneigbors](@ref)
"""
function ruralisation(c::Cell,n)
    # urban stays urban
    if c.state == 'U'
        @debug "Town gets older"
        return Cell(c.state, c.age + 1)
    # forest stays forest
    elseif c.state =='F'
        @debug "Forest grows older"
        return Cell(c.state, c.age + 1)
    # nature to forest
    elseif c.state == 'N'
        # within mutation treshold
        if rand() < 1/9 * tanh(c.age/20) * (1 + length(filter(x-> x.state =='F',n)) )
            @debug "NEW FOREST!!"
            return Cell('F',0)
        # otherwise just age
        else
            return Cell(c.state, c.age +1)
        end
    # agriculture to nature
    elseif c.state == 'A'
        if rand() < 1/8 * sum( map(x -> tanh(x.age / 10), filter(x-> x.state =='N',n) ) )
            @debug "NEW NATURE!!"
            return Cell('N',0)
        else
            return Cell(c.state, c.age +1)
        end
    end
end

"""
    update(c::Cell, n)

Given a cell and its neighbors, determine the next state
"""
function update(c::Cell, n)
    if rand() < 5/10
        return urbanisation(c, n)
    else
        return ruralisation(c,n)
    end
end

"""
    genworld(N::Int)

Generate an initial world (random composition) with a border.

statedist: default to uniform, but you can also use a custom distribution
"""
function genworld(N::Int;
                 statedist=Categorical(1/length(states).* ones(length(states))))
    @debug statedist
    world = zeros(Cell, N+2, N+2)
    world[2:end-1, 2:end-1] .= reshape([Cell(states[rand(statedist)], 7) for _ in 1:N^2], N, :)
    return world
end

"""
    newworld!(w::Array{Cell,2}, nw::Array{Cell,2})

From a given world, determine the next state and return it.
"""
function newworld!(w::Array{Cell,2}, nw::Array{Cell,2})
    for i in 2:size(w,1)-1
        for j in 2:size(w,2)-1
            nw[i,j] = update(w[i,j], getneigbors(w, i, j))
        end
    end
    w .= nw
    return w
end

# actual program
function main(  N::Int=10,
                d=Categorical(1/length(states).* ones(length(states)));
                nsteps::Int=100)
    w = genworld(N, statedist=d)
    nw = zeros(Cell, size(w))
    record = zeros(nsteps,4)
    for i in 1:nsteps
        newworld!(w, nw)
        record[i,1] = sum(map(x-> x.state == 'U', w)) / N^2
        record[i,2] = sum(map(x-> x.state == 'F', w)) / N^2
        record[i,3] = sum(map(x-> x.state == 'N', w)) / N^2
        record[i,4] = sum(map(x-> x.state == 'A', w)) / N^2
    end
    return record
end

# logger settings
Logging.global_logger(ConsoleLogger(stdout,Logging.Debug))
Logging.disable_logging(Logging.LogLevel(-1000))

# do it
@info "running main"
record = main(100)

# plotting (not was not required)
const cdict = Dict(65=>0, 70=>1, 78=>2, 85=>3, 98=>4)
const cmap = cgrad([    RGBA(222/255,166/255,53/255),       # agriculture brown
                        RGBA(80/255,118/255,40/255),        # forest green (dark)
                        RGBA(105/255,224/255,62/255),       # nature green (light)
                        RGBA(181/255,181/255,181/255),      # urban gray
                        RGBA(240/255,240/255,240/255)]);    # border off-white
const heatmapsettings = Dict(:yflip=>true, :color=>cmap, :axis=>false, :size=>(500,350), :clim=>(0,4));

"""
    plotworld(w::Array{Cell,2}; title::String="")

Make a heatmap of the world
"""
function plotworld(w::Array{Cell,2}; title::String="")
    return heatmap(map(x->cdict[Int(x.state)], w); title=title, heatmapsettings...)
end

"""
    testworld()

Generate a test world and show it
"""
function testworld()
    N = 7
    world = zeros(Cell, N+2, N+2)
    world[2:4, 2:4] = reshape([Cell('A', 0) for _ in 1:3^2], 3, :) # agricultre
    world[2:4, 6:8] = reshape([Cell('N', 0) for _ in 1:3^2], 3, :) # nature
    world[6:8, 2:4] = reshape([Cell('F', 0) for _ in 1:3^2], 3, :) # forest
    world[6:8, 6:8] = reshape([Cell('U', 0) for _ in 1:3^2], 3, :) # urban
    p = plotworld(world)
    title!("Testworld: [agriculture, nature; forest, urban]")
    return p
end

@info "plotting testworld"
display(testworld())
@info "plotting evolution"
plot(1:size(record,1), record, marker=:circle, ylims=(-0.05,1.05),
    label=["urban" "forest" "nature" "agriculture"],
    xlabel="time [years]", ylabel="% of land usage")


text:
module QuerySQLite

import Base: !, &, |, ==, !=, *, +, %, abs, Char, coalesce, collect, eltype, getproperty, length,
lowercase, in, isdone, isequal, isless, ismissing, iterate, IteratorSize, max, min,
occursin, rand, show, showerror, startswith, string, uppercase
using Base: Generator, NamedTuple, RefValue, SizeUnknown, tail
using Base.Meta: quot
import Base.Multimedia: showable
using DataValues: DataValue
import IteratorInterfaceExtensions: getiterator, isiterable
import MacroTools
using MacroTools: @capture
import QueryOperators
import SQLite
import SQLite: getvalue
using SQLite: columns, DB, execute!, generate_namedtuple, juliatype,
SQLITE_DONE, SQLITE_NULL, SQLITE_ROW, sqlite3_column_count, sqlite3_column_name,
sqlite3_column_type, sqlite3_step, sqlitevalue, Stmt, tables
using TableShowUtils: printdataresource, printHTMLtable, printtable
import TableTraits: isiterabletable

export Database

include("utilities.jl")
include("database.jl")
include("code_instead.jl")
include("iterate.jl")
include("model_row.jl")
include("translate.jl")
include("show.jl")

# To add support for a new function
# 1) Use `@code_instead` to specify the argument types
# 2) If the function modifies the model row, define model_row_call(::typeof(function), ...)
# 3a) If the function maps directly onto an SQL function, use @translate_default function SQL_function_symbol
# 3b) Otherwise, define translate_call(::typeof(function), ...)
# 4) If the SQL function uses special syntax, special case the SQLExpression show method

# There are two situtations in which Julia code cannot be translated into SQL
# 1) Patterns based on non-functions or non-overloadable functions, e.g. if, &&. Use if_else and & instead
# 2) Separate operations in SQL that are combined in Julia by dispatch, e.g. * for string concatention. Use `string` instead

end # module