text:
# Julia code
# Copyright (c) 2022 the SLICOTMath.jl developers
# Portions extracted from SLICOT-Reference distribution:
# Copyright (c) 2002-2020 NICONET e.V.

function run_mb03ud(datfile, io=stdout)
    NIN = 5
    NOUT = 6
    NMAX = 10
    LDA = NMAX
    LDQ = NMAX
    LDWORK = max( 1, 5*NMAX )
    A = Array{Float64,2}(undef, LDA,NMAX)
    Q = Array{Float64,2}(undef, LDQ,NMAX)
    SV = Array{Float64,1}(undef, NMAX)
    DWORK = Array{Float64,1}(undef, LDWORK)
    f = open(datfile,"r")
    readline(f)
    vs = split(readline(f))
    N = parse(BlasInt, vs[1])
    JOBQ =  vs[2][1]
    JOBP =  vs[3][1]
if ( N<0 || N>NMAX )
else

    vs = String[]
    _isz,_jsz = (N,N)
    while length(vs) < _isz*_jsz
        append!(vs, replace.(split(readline(f)),'D'=>'E'))
    end
    for i in 1:_isz
       _i0 = (i-1)*_jsz
       A[i,1:_jsz] .= parsex.(Float64, vs[_i0+1:_i0+_jsz])
    end
# interp call 1

    INFO = SLICOT.mb03ud!(JOBQ, JOBP, N, A, Q, SV)
    @test INFO == 0
    INFO == 0 || return

if ( INFO!=0 )
else
# interp output 1
    println(io,"SV:")
    _nr = N
    show(io,"text/plain",SV[1:_nr])
    println(io,)

if ( LSAME( JOBP, 'V' ) )
# interp output 2
    println(io,"A:")
    _nc = N
    _nr = N
    show(io,"text/plain",A[1:_nr,1:_nc])
    println(io,)

end # if
if ( LSAME( JOBQ, 'V' ) )
# interp output 3
    println(io,"Q:")
    _nc = N
    _nr = N
    show(io,"text/plain",Q[1:_nr,1:_nc])
    println(io,)

end # if
end # if
end # if
    close(f)
end # run_X()


text:
function integrate( b, w, P, nstep; fplot=false, fsave=false, nsave=10)

    # Initial guess for h
    h = ones(P.nx,P.ny).*P.p./infilt(b,P)
    In=infilt(b,P)
    steadyh!(h,In, P)

    # Init fields for integral approximation
    σ=initsigma(1.0, 1. +(P.η)*1.1, P.nsigma)
    fg=initfg(σ, P.nx, P.dx)
    bint=zeros(P.nx,P.ny)
    wint=zeros(P.nx,P.ny)
    global iint=0

    u=cat(b,w,dims=3)
    ttot=0.
    for i=1:nstep
        @printf("t=%4.3f <b>=%1.4f (%1.4f, %1.4f) <w>=%1.4f (%1.4f, %1.4f)", ttot, mean(u[:,:,1]), minimum(u[:,:,1]), maximum(u[:,:,1]), mean(u[:,:,2]), minimum(u[:,:,2]), maximum(u[:,:,2]))

        ttot+=P.dt

        # Integrate over chunck P.dt
        prob = ODEProblem(rhs_stat!, u, (0.,P.dt), (P, fg, σ, h, bint, wint))
        @time sol=solve(prob, save_everystep=false, save_start=false);
        u=sol[1]

        # Plot solutions
        if(fplot==true)
            plotbwh(u[:,:,1], u[:,:,2], h, P, ttot)
        end

        # Save solution
        if(fsave & (mod(i, nsave)==0))
            @printf("Saving step %d, time %4.5f\n", i, ttot)
            open(P.outfile, "a") do io
                writedlm(io, reshape(cat(u,h,dims=3), P.nx*P.ny, 3))
            end
        end

    end

    @printf("t=%4.3f <b>=%1.4f (%1.4f, %1.4f) <w>=%1.4f (%1.4f, %1.4f)\n", ttot, mean(u[:,:,1]), minimum(u[:,:,1]), maximum(u[:,:,1]), mean(u[:,:,2]), minimum(u[:,:,2]), maximum(u[:,:,2]))

    return u[:,:,1], u[:,:,2], h
end

ext: jl

max_stars_repo_name: jhardenberg/bwh

max_issues_repo_name: jhardenberg/bwh

max_forks_repo_name: jhardenberg/bwh


text:
"""
    abstract type AbstractChainElement{S<:AbstractCell, F}

Abstract type for representing an element of a chain. Create subtypes of this if there is a
more efficient way to store a simplex and coefficient than as two fields.

An `AbstractChainElement{S, F}` behaves like a `S` when doing comparisons and like a `F`
when doing arithmetic.

Coefficient values must be in a field - they must support `+`, `-`, `*`, `/`, `inv`, `zero`,
`one` and `iszero`.

# Interface

* `AbstractChainElement{S, F}(sx::S)` - default to element with `coefficient` equal to
  `sign(sx)`.
* `AbstractChainElement{S, F}(::S, coefficient)` - convert `coefficient` to `F` if
  necessary.
* `simplex(::ChainElement)`
* `coefficient(::ChainElement)`
"""
abstract type AbstractChainElement{S<:AbstractCell,F<:Number} end

for op in (:+, :-)
    @eval function (Base.$op)(ce1::C, ce2::C) where {C<:AbstractChainElement}
        @boundscheck ce1 == ce2 || throw(ArgumentError("simplices don't match"))
        return C(simplex(ce1), $op(coefficient(ce1), coefficient(ce2)))
    end
end
for op in (:oneunit, :zero, :-, :+)
    @eval function (Base.$op)(elem::C) where {C<:AbstractChainElement}
        return C(simplex(elem), $op(coefficient(elem)))
    end
end
function Base.:*(elem::C, x::F) where {F<:Number,C<:AbstractChainElement{<:Any,F}}
    return C(simplex(elem), coefficient(elem) * x)
end
function Base.:*(x::F, elem::C) where {F<:Number,C<:AbstractChainElement{<:Any,F}}
    return C(simplex(elem), x * coefficient(elem))
end
function Base.:/(elem::C, x::F) where {F<:Number,C<:AbstractChainElement{<:Any,F}}
    return C(simplex(elem), coefficient(elem) / x)
end
Base.one(::Type{<:AbstractChainElement{<:Any,F}}) where {F} = one(F)
Base.one(::AbstractChainElement{<:Any,F}) where {F} = one(F)
Base.iszero(elem::AbstractChainElement) = iszero(coefficient(elem))

Base.hash(elem::AbstractChainElement, u::UInt64) = hash(simplex(elem), u)
function Base.:(==)(elem1::AbstractChainElement, elem2::AbstractChainElement)
    return simplex(elem1) == simplex(elem2)
end
function Base.isless(elem1::AbstractChainElement, elem2::AbstractChainElement)
    return isless(simplex(elem1), simplex(elem2))
end

# Make chain elements useful when they come out as representatives.
birth(elem::AbstractChainElement) = birth(simplex(elem))
index(elem::AbstractChainElement) = index(simplex(elem))
Base.sign(elem::AbstractChainElement) = sign(coefficient(elem))
vertices(elem::AbstractChainElement) = vertices(simplex(elem))

simplex_type(::Type{<:AbstractChainElement{S}}) where {S} = S
coefficient_type(::Type{<:AbstractChainElement{<:Any,F}}) where {F} = F
simplex_type(::C) where {C<:AbstractChainElement} = simplex_type(C)
coefficient_type(::C) where {C<:AbstractChainElement} = coefficient_type(C)

# Make chain elements kinda similar to `Pair`s
function Base.getindex(elem::AbstractChainElement, i)
    if i == 1
        return simplex(elem)
    elseif i == 2
        return coefficient(elem)
    else
        throw(BoundsError(elem, i))
    end
end

Base.eltype(::Type{<:AbstractChainElement{S,F}}) where {S,F} = Union{S,F}
Base.iterate(elem::AbstractChainElement, i=1) = i > 2 ? nothing : (elem[i], i + 1)
Base.firstindex(::AbstractChainElement) = 1
Base.lastindex(::AbstractChainElement) = 2
Base.length(::AbstractChainElement) = 2

function Base.show(io::IO, elem::AbstractChainElement)
    return print(io, simplex(elem), " => ", coefficient(elem))
end

function Base.convert(::Type{C}, simplex::S) where {S,C<:AbstractChainElement{S}}
    return C(simplex)
end
function Base.convert(::Type{S}, elem::AbstractChainElement{S}) where {S<:AbstractCell}
    return simplex(elem)
end
function Base.convert(
    ::Type{C}, elem::AbstractChainElement{S,F}
) where {S,F,C<:AbstractChainElement{S,F}}
    return C(simplex(elem), coefficient(elem))
end

"""
    chain_element_type(simplex, coefficient)

Get the type of `AbstractChainElement` that is to be used for packing `simplex` and
`coefficient`. Overload this function to specialize element types for (simplex, field) type
pairs.
"""
chain_element_type(::S, ::F) where {S,F} = chain_element_type(S, F)

"""
    ChainElement{S<:AbstractCell, F} <: AbstractChainElement{S, F}

The default, basic subtype of [`AbstractChainElement`](@ref).
"""
struct ChainElement{S<:AbstractCell,F} <: AbstractChainElement{S,F}
    simplex::S
    coefficient::F

    function ChainElement{S,F}(simplex::S, coefficient=one(F)) where {S,F}
        return new{S,F}(abs(simplex), sign(simplex) * F(coefficient))
    end
end

simplex(elem::ChainElement) = elem.simplex
coefficient(elem::ChainElement) = elem.coefficient
chain_element_type(::Type{S}, ::Type{F}) where {S,F} = ChainElement{S,F}

"""
    PackedElement{S<:Simplex, F<:Mod} <: AbstractChainElement{S, F}

Like `ChainElement` for `Simplex` and `Mod`. Packs the coefficient value into top `M` bits
of index. Used internally by `Chain`.
"""
struct PackedElement{S<:Simplex,F<:Mod,M,U,T} <: AbstractChainElement{S,F}
    index_coef::U
    birth::T

    function PackedElement{S,F,M,U,T}(
        simplex::S, coefficient=one(F)
    ) where {M,U,T,S<:Simplex{<:Any,T},F<:Mod{M}}
        idx = simplex.index
        coef = F(coefficient) * sign(idx)
        uidx = U(abs(idx))
        index_coef = U(Int(coef)) << (sizeof(U) * 8 - n_bits(Val(M))) | uidx

        return new{S,F,M,U,T}(index_coef, birth(simplex))
    end
end

"""
    n_bits(::Val{M})

Get numer of bits needed to represent number mod `M`.
"""
@inline n_bits(::Val{M}) where {M} = 8 * sizeof(Int) - leading_zeros(M)

function simplex(pe::PackedElement{S,<:Any,M,U}) where {S,M,U}
    mask = typemax(U) << n_bits(Val(M)) >> n_bits(Val(M))
    return S(pe.index_coef & mask, pe.birth)
end
function coefficient(elem::PackedElement{<:Any,F,M,U}) where {F,U,M}
    return F(elem.index_coef >> (sizeof(U) * 8 - n_bits(Val(M))), false)
end

function chain_element_type(
    ::Type{S}, ::Type{F}
) where {M,I,T,S<:Simplex{<:Any,T,I},F<:Mod{M}}
    if n_bits(Val(M)) ≤ 8
        return PackedElement{S,F,M,unsigned(I),T}
    else
        return ChainElement{S,F}
    end
end

function index_overflow_check(
    ::Type{S}, ::Type{F}, nv, message
) where {T,I,F<:Mod,S<:Simplex{<:Any,T,I}}
    vertices = ntuple(i -> I(nv - i + 1), length(S))
    element = chain_element_type(S, F)(S(index(I.(vertices)), zero(T)))

    index_I = index(element)
    index_big = index(BigInt.(vertices))

    if index_I ≠ index_big
        throw(OverflowError(message))
    end
end

function index_overflow_check(::Type{<:AbstractCell}, _, _, _)
    return nothing
end


text:
# This file was generated, do not modify it. # hide
get_mse(lr) = rms(MLJ.predict(lr, rows=test), y[test])^2

@show get_mse(lr1)
@show get_mse(lr2)
@show get_mse(lr3)

max_issues_repo_name: alan-turing-institute/MLJTutorials

max_forks_repo_name: alan-turing-institute/MLJTutorials


text:
using BinaryBuilder

# Collection of sources required to build libRmath
sources = [
    "https://github.com/JuliaLang/Rmath-julia/archive/v0.2.0.tar.gz" =>
    "087ada2913c5401c5772cde1606f9924dcb159f1c9d755630dcce350ef8036ac",
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
cd Rmath-julia-*/

# The Rmath-julia makefile is kind of silly
MAKEVARS="fPIC=-fPIC prefix=$prefix"
libdir=$prefix/lib

# It needs to be told it's on Windows
if [[ ${target} == *mingw* ]]; then
    MAKEVARS="${MAKEVARS} OS=Windows_NT"
    libdir=$prefix/bin
fi

# Build it
make ${MAKEVARS} CC="$CC" CXX="$CXX" -j${nproc}

# Install it
mkdir -p $libdir $prefix/include
mv src/libRmath-julia*.${dlext} $libdir/
mv include/* $prefix/include/
"""

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

# The products that we will ensure are always built
products(prefix) = [
    LibraryProduct(prefix, "libRmath", :libRmath)
]

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

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, "libRmath", sources, script, platforms, products, dependencies)


text:
struct UniformDiscrete <: Distribution{Int} end

"""
    uniform_discrete(low::Integer, high::Integer)

Sample an `Int` from the uniform distribution on the set {low, low + 1, ..., high-1, high}.
"""
const uniform_discrete = UniformDiscrete()

function logpdf(::UniformDiscrete, x::Int, low::Integer, high::Integer)
    d = Distributions.DiscreteUniform(low, high)
    Distributions.logpdf(d, x)
end

function logpdf_grad(::UniformDiscrete, x::Int, lower::Integer, high::Integer)
    (nothing, nothing, nothing)
end

function random(::UniformDiscrete, low::Integer, high::Integer)
    rand(Distributions.DiscreteUniform(low, high))
end
is_discrete(::UniformDiscrete) = true

(::UniformDiscrete)(low, high) = random(UniformDiscrete(), low, high)

has_output_grad(::UniformDiscrete) = false
has_argument_grads(::UniformDiscrete) = (false, false)

export uniform_discrete


text:
#Collection of functions relating to finding fixed points of reaction networks, as well as the creation of bifurcation diagrams.

### Functions used during the construction of the reaction network. ###

#Equilibrate content structure, with all the stuff required for homotopy continuation base equilibration calculations.
mutable struct EquilibrateContent
    make_polynomial::Function
    constraints::Vector{Polynomial{true,Float64}}
    homotopy_continuation_templates::Vector{NamedTuple{(:p, :sol),Tuple{Vector{Complex{Float64}},Vector{Vector{Complex{Float64}}}}}}
    equilibrium_polynomial::Union{Vector{Polynomial{true,Float64}},Nothing}
    is_polynomial_system::Bool
    polyvars::NamedTuple{(:u, :t, :p),Tuple{Vector{PolyVar{true}},PolyVar{true},Vector{PolyVar{true}}}}
end

function EquilibrateContent(make_polynomial,n_vars,n_params)
    pvus = (HomotopyContinuation.@polyvar internal___polyvar___u[1:n_vars])[1]
    pvps = (HomotopyContinuation.@polyvar internal___polyvar___p[1:n_params])[1]
    pvt = (HomotopyContinuation.@polyvar internal___polyvar___t)[1]
    is_polynomial_system = try
        make_polynomial(pvus,pvt,pvps,internal___var___paramvals=fill(1,length(pvps)))
        true
    catch; false; end
    equilibrium_polynomial = try
        tmp_pol = make_polynomial(pvus,pvt,pvps);
        (typeof(tmp_pol[1])<:Polynomial) ? tmp_pol : map(pol->pol.num,tmp_pol)
    catch; nothing; end

    return EquilibrateContent(
        make_polynomial,
        Vector{Polynomial{true,Float64}}(),
        Vector{NamedTuple{(:p, :sol),Tuple{Vector{Complex{Float64}},Vector{Vector{Complex{Float64}}}}}}(),
        equilibrium_polynomial,
        is_polynomial_system,
        (u=pvus,t=pvt,p=pvps),
    )
end

abstract type AbstractBifurcationSolver end
struct SimpleHCBifurcationSolver <: AbstractBifurcationSolver end
struct HCBifurcationSolver <: AbstractBifurcationSolver end

abstract type AbstractSteadyStateSolver end
struct HCSteadyStateSolver <: AbstractSteadyStateSolver end

#Various Functions for accessing and updating an EquilibrateContent structure.
get_equi_poly(rn::DiffEqBase.AbstractReactionNetwork) = (return rn.equilibrate_content.equilibrium_polynomial)
has_equi_poly(rn::DiffEqBase.AbstractReactionNetwork) = (return rn.equilibrate_content.equilibrium_polynomial != nothing)
is_poly(rn::DiffEqBase.AbstractReactionNetwork) = (return rn.equilibrate_content.is_polynomial_system)
get_polyvars(rn::DiffEqBase.AbstractReactionNetwork) = (return rn.equilibrate_content.polyvars)
get_hc_templates(rn::DiffEqBase.AbstractReactionNetwork) = (return rn.equilibrate_content.homotopy_continuation_templates)
has_hc_templates(rn::DiffEqBase.AbstractReactionNetwork) = (return length(rn.equilibrate_content.homotopy_continuation_templates)!=0)

push_hc_template!(rn::DiffEqBase.AbstractReactionNetwork,p_template,sol_template) = (push!(rn.equilibrate_content.homotopy_continuation_templates,(p=p_template,sol=sol_template)))
push_constraint!(rn::DiffEqBase.AbstractReactionNetwork,constraint) = (push!(rn.equilibrate_content.constraints,constraint); has_equi_poly(rn) && push!(rn.equilibrate_content.equilibrium_polynomial,constraint);)
set_equi_poly!(rn::DiffEqBase.AbstractReactionNetwork,t;kwargs...) = (rn.equilibrate_content.equilibrium_polynomial = derationalise_polys(rn.equilibrate_content.make_polynomial(rn.equilibrate_content.polyvars.u,t,rn.equilibrate_content.polyvars.p;kwargs...)))
derationalise_polys(polys) = (typeof(polys[1])<:Polynomial) ? (return polys) : (return map(pol->pol.num,polys))
add_constraints!(rn::DiffEqBase.AbstractReactionNetwork) = (has_equi_poly(rn) && foreach(constraint -> push!(rn.equilibrate_content.equilibrium_polynomial,constraint), rn.equilibrate_content.constraints))
number_of_constraints(rn::DiffEqBase.AbstractReactionNetwork) = return length(rn.equilibrate_content.constraints)

reset_equi_poly!(rn::DiffEqBase.AbstractReactionNetwork) = (rn.equilibrate_content.equilibrium_polynomial = nothing)
reset_hc_templates!(rn::DiffEqBase.AbstractReactionNetwork) = (rn.equilibrate_content.homotopy_continuation_templates = Vector{NamedTuple{(:p, :sol),Tuple{Vector{Complex{Float64}},Vector{Vector{Complex{Float64}}}}}}())

check_is_polynomial(rn::DiffEqBase.AbstractReactionNetwork) = (!is_poly(rn)) && (error("This reaction network does not correspond to a polynomial system. Some of the reaction rate must contain non polynomial terms."))
check_has_polynomial(rn::DiffEqBase.AbstractReactionNetwork) = (!has_equi_poly(rn)) && (error("No equilibrium polynomial have been created. Please use the @fix_concentration macro to do so."))

#Generates the expression which is then converted into a function which generates polynomials for given parameters. Parameters not given can be given at a later stage, but parameters in the exponential must be given here. Also checks whenever the system is polynomial.
function get_equilibration(params::Vector{Symbol}, reactants::OrderedDict{Symbol,Int}, f_expr::Vector{ExprValues})
    func_body = Expr(:block)
    push!(func_body.args,Expr(:if,:(length(internal___var___paramvals) != 0),Expr(:block)))
    for i = 1:length(params)
        push!(func_body.args[1].args[2].args, :((internal___var___paramvals___exemption!=$i) && ($(params[i])=(isinteger(internal___var___paramvals[$i]) ? Int64(internal___var___paramvals[$i]) : internal___var___paramvals[$i]))))
    end
    push!(func_body.args,:([]))
    foreach(poly -> push!(func_body.args[2].args, recursive_replace!(poly,(reactants,:internal___polyvar___u),(OrderedDict(:t=>1),:internal___polyvar___t))), deepcopy(f_expr))
    func_expr = :((;TO___BE___REMOVED=to___be___removed) -> $(deepcopy(func_body)))
    foreach(i -> push!(func_expr.args[1].args[1].args,Expr(:kw,params[i],:(internal___polyvar___p[$i]))), 1:length(params))
    push!(func_expr.args[1].args[1].args,Expr(:kw,:internal___var___paramvals,:([])))
    push!(func_expr.args[1].args[1].args,Expr(:kw,:internal___var___paramvals___exemption,-1))
    deleteat!(func_expr.args[1].args[1].args,1)
    push!(func_expr.args[1].args,:internal___polyvar___u)
    push!(func_expr.args[1].args,:internal___polyvar___t)
    push!(func_expr.args[1].args,:internal___polyvar___p)
    return func_expr
end

### Functions required for the preparation to use homotopy continuation to find fixed points. ###

#--- Fixes parameters, required when some exponents contain parameters. --#
#Calls the internal make_polynomial function to create a polynomial. Used to model creation (will succced unless there is parameters in the exponents). Input include list of parameters to fix. All parameters occuring as exponent must be here (none other is required at this stage). Will insert fixed concentrations if such exists.
"""
    fix_parameters(reaction_network, [parameter_values]; t=get_polyvars(rn).t, kwargs...)

Fixes certain parameters of a reaction network, preparing it for homotopy continuation.

In the polynomial needed to be solved by homotopy continuation to find steady states solutions, both varriables and parameters are initially represented as polynomial varriables. If a parameter exists in an exponent (e.g. X^n), that parameter might be needed to be fixed to a desired value before solving. Fixing is not required, but if you solve a polynomial several times, it might increase performance.

## args
-  reaction_network: a reaction network.
-  parameter_values: An entire vector can be given, fixing all parameters. Generaly you do not want to do this.
-  t=get_polyvars(rn).t: If you want to fix the time value.

## kwargs
-  How parameters generally are fixed. Write e.g. n=4 to fix the parameter value of n.
"""
function fix_parameters(rn::DiffEqBase.AbstractReactionNetwork, params=Vector{Number}()::Vector{Number}; t=get_polyvars(rn).t, full_vector_exemption=-1, kwargs...)
    check_is_polynomial(rn)
    set_equi_poly!(rn,t;internal___var___paramvals=params,internal___var___paramvals___exemption=full_vector_exemption,kwargs...)
    add_constraints!(rn)
    return nothing
end

#--- Adds additional constraints to the system, typically a fixed concentration such as in the system X <--> Y.
#In some networks some linear combinations of concentrations remain fixed. This supplies this information directly to the network.
"""
    @add_contraint(reaction_network, constraints)

Some networks (like X ↔ Y) have an infinite number of potential steady states, depending on initial conditions. For these some additional contraint(s) is required. Typically in the form of fixed concentrations.

## args
-  reaction_network: a reaction network.
-  constraints: Constraints, typically given as X + Y = C. Both parameters, varriables, and numbers can be given. However, parameters much be declared previously in the reaction network.
"""
macro add_constraint(rn::Symbol, constraints::Expr...)
    func_expr = Expr(:escape,:(internal___add___constraint!($rn)))
    foreach(constraint -> push!(func_expr.args[1].args,recursive_replace_vars!(balance_poly(constraint), rn)),constraints)
    return func_expr
end
#In some networks some linear combinations of concentrations remain fixed. This supplies this information directly to the network.
"""
    @add_contraints(reaction_network, constraints)

Same as @add_constraint, but a more natural way of adding several contraints such as:
@add_contraints rn begin
    X1 + Y1 = C1
    X2 + Y2 = C2
end

## args
-  reaction_network: a reaction network.
-  constraints: Constraints, typically given as X + Y = C. Both parameters, varriables, and numbers can be given. However, parameters much be declared previously in the reaction network.
"""
macro add_constraints(rn::Symbol, constraints::Expr...)
    func_expr = Expr(:escape,:(internal___add___constraint!($rn)))
    foreach(constraint -> push!(func_expr.args[1].args,recursive_replace_vars!(balance_poly(constraint), rn)),MacroTools.striplines(constraints...).args)
    return func_expr
    output = Expr(:block)
    for constraint in MacroTools.striplines(constraints...).args
        push!(output.args,Expr(:escape,:(internal___add___constraint!($rn))))
        push!(output.args[end].args[1].args,recursive_replace_vars!(balance_poly(constraint), rn))
    end
    return output
end
#Changes a polynomial expression of the form p(x) = q(x) to p(x) - q(x). Used by add_constraint(s).
function balance_poly(poly::Expr)
    (poly.head != :call) && (return :($(poly.args[1])-$(poly.args[2])))
    return poly
end
#Complicated function which replaces the named varriables in the expression supplied to @fix_concentration with the corresponding polyvar stored in the reaction network.
function recursive_replace_vars!(expr::Union{ExprValues,LineNumberNode}, rn::Symbol)
    if (typeof(expr) == Symbol)&&(expr!=:+)&&(expr!=:-)&&(expr!=:*)&&(expr!=:^)&&(expr!=:/)
        return :(in($(QuoteNode(expr)),$(rn).syms) ? $(rn).equilibrate_content.polyvars.u[$(rn).syms_to_ints[$(QuoteNode(expr))]] : $(rn).equilibrate_content.polyvars.p[$(rn).params_to_ints[$(QuoteNode(expr))]])
    elseif typeof(expr) == Expr
        foreach(i -> expr.args[i] = recursive_replace_vars!(expr.args[i], rn), 1:length(expr.args))
    end
    return expr
end
#Function which does the actual work of adding the fixed concentration information. Not meant to be called directly, but is called through the fixed_concentration macro.
function internal___add___constraint!(rn::DiffEqBase.AbstractReactionNetwork, constraints::Polynomial...)
    check_is_polynomial(rn)
    foreach(constraint -> push_constraint!(rn,constraint),constraints)
    remove_replicate_polys!(rn)
end
#Replicate polynomials might be bad for HomotopyContinuation, this function removes them when they might occur.
function remove_replicate_polys!(rn::DiffEqBase.AbstractReactionNetwork)
    for i = length(get_equi_poly(rn)):-1:2, j = (i-1):-1:1
        (get_equi_poly(rn)[i]==get_equi_poly(rn)[j])&&deleteat!(get_equi_poly(rn),i)
        (get_equi_poly(rn)[i]==-get_equi_poly(rn)[j])&&deleteat!(get_equi_poly(rn),i)
    end
end

#--- Makes homotopy continuation templates to be used at a later stage when solving for fixed points ---#
#Can be called separatly, but will otherwise be called first time a steady state is to be found. Solves the system once using random parameters. Saves the solution as a template to be used for further solving.
"""
    make_hc_template(reaction_network)

    Makes a template for solving for the given reaction networks steady states using homotopy continuation. Generally not needed as this is genrated automatically when netowkr is first solved *if one does not already exists). However, it can be sued to replace the current template, for whatever reason.

## args
-  reaction_network: a reaction network.
"""
function make_hc_template!(rn::DiffEqBase.AbstractReactionNetwork)
    check_is_polynomial(rn)
    check_has_polynomial(rn)
    p_template = randn(ComplexF64, length(rn.params))
    f_template = DynamicPolynomials.subs.(get_equi_poly(rn), Ref(get_polyvars(rn).p => p_template))
    solution_template = HomotopyContinuation.solutions(HomotopyContinuation.solve(f_template, show_progress=false))
    reset_hc_templates!(rn)
    push_hc_template!(rn,p_template,solution_template)
end
#Similar to make_hc_template, but if a template is already presence this will add another one (as opposed to make which creates one and discard previous ones). Having several templates might add some robustness.
"""
    add_hc_template(reaction_network)

    Adds another hc template to the network. Having more HC templates might increase run time a little bit, but might give some additional certainity that all solutions are found (this should generally not be a problem anyway)

## args
-  reaction_network: a reaction network.
"""
function add_hc_template!(rn::DiffEqBase.AbstractReactionNetwork)
    p_template = randn(ComplexF64, length(rn.params))
    f_template = DynamicPolynomials.subs.(get_equi_poly(rn), Ref(get_polyvars(rn).p => p_template))
    solution_template = HomotopyContinuation.solutions(HomotopyContinuation.solve(f_template, show_progress=false))
    if (isempty(get_hc_templates(rn))) || (length(solution_template) == length(get_hc_templates(rn)[1].sol))
        push_hc_template!(rn,p_template,solution_template)
    elseif length(solution_template) > length(get_hc_templates(rn)[1].sol)
        reset_hc_templates!(rn)
        push_hc_template!(rn,p_template,solution_template)
        @warn "Tried to add another homotopy continuation template to the set. However the new template contained more solutions than the previous one(s). Previous template(s) discarded."
    else
        @warn "Tried to add another homotopy continuation template to the set. However the new template contained less solutions than the previous one(s). No new template added."
    end
end
#Macro running the HC template function.
"""
    @make_hc_template(reaction_network)

    Calls the "make_hc_template function on the given reaction network."

## args
-  reaction_network: a reaction network.
"""
macro make_hc_template(rn::Symbol)
    return Expr(:escape,:(make_hc_template!($rn)))
end
#Macro running the HC template function.
"""
    @make_hc_template(reaction_network)

    Calls the "add_hc_template function on the given reaction network."

## args
-  reaction_network: a reaction network.
"""
macro add_hc_template(rn::Symbol)
    return Expr(:escape,:(add_hc_template!($rn)))
end

### Functions for finding single steady states of fixed points, and for analysing their stability. ###

#Finds the steady states of a reaction network
"""
    steady_states(reaction_network, [solver], [parameter_values]; kwargs...)

Finds the steady states of a given reaction network, at given parameter values.

## args
-  reaction_network: a reaction network.
-  solver (optional): a subtype of AbstractSteadyStateSolver that specifies how to solve the bifurcation diagram. Default is HCSteadyStateSolver (currently only solver avaiable).
-  parameter_values: a vector which specifies the parameter values for which the steady states are to be found. Not required in case the reaction network lacks steady states.

## kwargs
 - potential arguments for steady state solvers (however, only current solver does not have any such arguments).
"""
function steady_states(rn::DiffEqBase.AbstractReactionNetwork, args...)
    return steady_states(rn, HCSteadyStateSolver(), args...)
end
#Finds steady states of a system using homotopy continuation.
function steady_states(rn::DiffEqBase.AbstractReactionNetwork, solver::HCSteadyStateSolver, p=Vector{Float64}())
    using_temp_poly  = initialise_solver!(rn,p)
    (length(p)==0) && (return positive_real_solutions(HomotopyContinuation.solutions(HomotopyContinuation.solve(get_equi_poly(rn),show_progress=false))))
    result = hc_solve_at(rn,p)
    using_temp_poly && finalise_solver!(rn)
    return positive_real_solutions(result)
end
# Solves the reaction networks at the specific parameter values using Homotopy Continuation.
function hc_solve_at(rn::DiffEqBase.AbstractReactionNetwork,p::Vector{Float64})
    for hc_template in get_hc_templates(rn)
        result = HomotopyContinuation.solutions(HomotopyContinuation.solve(get_equi_poly(rn), hc_template.sol, parameters=get_polyvars(rn).p, p₁=hc_template.p, p₀=p, show_progress=false))
        (length(result) == length(hc_template.sol)) && (return result)
    end
    #@warn "While solving the system using homotopy continuation some solutions were lost."
    best_length = 0; best_solution = [];
    for hc_template in get_hc_templates(rn)
        result = HomotopyContinuation.solutions(HomotopyContinuation.solve(get_equi_poly(rn), hc_template.sol, parameters=get_polyvars(rn).p, p₁=hc_template.p, p₀=p, show_progress=false))
        (length(result) > best_length) && (best_length = length(result); best_solution = result;)
    end
    return best_solution
end

#Prepares the reaction network for being solved with Homotopy Continuation, in case some initialising operations have not been done yet.
function initialise_solver!(rn::DiffEqBase.AbstractReactionNetwork, p::Vector{Float64}, bifurcation_exception_parameter::Int=-1)
    check_is_polynomial(rn)
    using_temp_poly = !has_equi_poly(rn)
    using_temp_poly && fix_parameters(rn, p, full_vector_exemption=bifurcation_exception_parameter)
    !has_hc_templates(rn) && make_hc_template!(rn)
    return using_temp_poly
end
# In case a temporary equilibrium polynomial were used, this one resets it.
function finalise_solver!(rn::DiffEqBase.AbstractReactionNetwork)
    reset_equi_poly!(rn); reset_hc_templates!(rn);
end
# Extracts the biological plausible solutions (real and non negative).
function positive_real_solutions(result)
    tmp_sol = filter(result) do y
        all(yᵢ -> abs(imag(yᵢ)) < 0.001, y)&&all(yᵢ -> real(yᵢ) ≥ -0.001, y)
    end
    return map(u->real.(u),tmp_sol)
end

#Returns a boolean which will be true if the given fixed point is stable.
"""
    stability(solution,reaction_network, parameter_values,  t=0.)

    Gives the steady stability of the steady state of a certain network. Returns true if the point is stable, false if it is not.

## args
- solution: the values of the reactants in the steady state.
- parameter_values: the parameter values in for which the steady state was calculated.
- reaction_network: a reaction network.
- t (optional): the value of t, for which the stability is calculated.
"""
function stability(solution::Vector{Float64}, rn::DiffEqBase.AbstractReactionNetwork, p::Vector{Float64}, t=0.::Float64)
    jac = zeros(length(rn.syms),length(rn.syms))
    rn.jac(jac,solution,p,t)
    return sort(real.(eigen(jac).values))[end-number_of_constraints(rn)]<0.
end
#Runs stability on an vector of steady states.
"""
    stability(solutions,reaction_network, parameter_values,  t=0.)

    Gives the steady stability of the steady state of a certain network.

## args
- solutions: An array of steady states (which each in turn, is an array of reactant concentrations in the steady state).
- parameter_values: the parameter values in for which the steady state was calculated.
- reaction_network: a reaction network.
- t (optional): the value of t, for which the stability is calculated.
"""
function stability(solutions::Vector{Vector{Float64}}, rn::DiffEqBase.AbstractReactionNetwork, p::Vector{Float64}, t=0.::Float64)
    return map(sol->stability(sol,rn,p,t), solutions)
end

### Structures holding information about bifurcation diagrams. ###

### Bifurcation Diagrams Structures###

#Contains information from a single bifurcation line. A bifurcation plot would typically contain several such lines. These are the results from a method which tracks a solution through parameter space.
struct BifurcationPath
    p_vals::Vector{Float64}
    vals::Vector{Vector{Float64}}
    jac_eigenvals::Vector{Vector{ComplexF64}}
    stability_types::Vector{Int8}
    length::Int
end
#A bifurcation diagram, contains a set of bifurcation paths.
struct BifurcationDiagram
    param::Symbol
    range::Tuple{Float64,Float64}
    paths::Vector{BifurcationPath}
end

#This is a single point in a grid like bifurcation scheme. It contains a parameter value and corresponding fixed points, as well as stability information.
struct BifurcationPoint
    vals::Vector{Vector{Float64}}
    jac_eigenvals::Vector{Vector{ComplexF64}}
    stability_types::Vector{Int8}
end
#Contains a grid like bifurgation diagrams with a vector of bifurcation points.
struct BifurcationGrid
    param::Symbol
    range::AbstractRange
    grid_points::Vector{BifurcationPoint}
    length::Int
end
#Contains a grid like bifurgation diagrams over 2 parameters. The bifurcation points are contained in a matrix.
struct BifurcationGrid2D
    param1::Symbol
    range1::AbstractRange
    param2::Symbol
    range2::AbstractRange
    grid_points::Matrix{BifurcationPoint}
    size::Tuple{Int64,Int64}
end

#A 2d hybrid of grid and diagram. Contains a 1d grid of bifurcation diargam.
struct BifurcationGridDiagram
    param1::Symbol
    range1::AbstractRange
    param2::Symbol
    range2::Tuple{Float64,Float64}
    bifurcation_diagrams::Vector{BifurcationDiagram}
    length::Int
end

### Functions called to generate bifurcation diagrams. ###

#Generates a bifurcation diagram for a given reaction network and between two given parameter values. Does so by tracking paths

#Generates a bifurcation diagram.
"""
    bifurcations(reaction_network, [solver], parameter_values, parameter_symbol, parameter_range; kwargs...)

Get a bifurcation diagram of the specified system.

## args
-  reaction_network: a reaction network.
-  solver (optional): a subtype of AbstractBifurcationSolver that specifies how to solve the bifurcation diagram. Default is HCBifurcationSolver.
-  parameter_values: a vector which specifies the point in parameter space around which the bifurcation diagram is evaluated.
-  parameter_symbol: the parameter which is varied, given as a symbol.
-  parameter_range: the range over which the specified parameter is varied.

## kwargs
- stepsize = (parameter_range[2] - parameter_range[1])/200: The maximum distance between two parameter values for which steady states are calculated. Smaller values will give a smoother diagram.
- Δp = (parameter_range[2] - parameter_range[1])/200: The distance to jump after finding a bifurcation. After discovering a bifurcation, the solver jumps ahead a distance `Δp` and starts back-tracking. If this distance is too small, the method may error or cause visible artifiacts in the bifurcation diagram. If it is too large, then you might jump over another bifurcation and it will be missed (this should also be farily obvious in a plot). Only applicable for the HCBifurcationSolver.
- Δu = 1e-7: (alternatively reduced to a tenth of the minimum value of any initial solutions) used by algorithm to determine if two points are identical. Set this to less than the expected minimum distance between two different solutions.
"""
function bifurcations(rn::DiffEqBase.AbstractReactionNetwork, args...; kwargs...)
    return bifurcations(rn, HCBifurcationSolver(), args...; kwargs...)
end
function bifurcations(rn::DiffEqBase.AbstractReactionNetwork, solver::AbstractBifurcationSolver, p, param::Symbol, range; kwargs...)
    return BifurcationDiagram(param,range,solve_bifurcation(solver,rn,p,param,range; kwargs...))
end

#Generates a grid of bifurcation points, using a given steady state method.
"""
    bifurcation_grid(reaction_network, [solver], parameter_values, parameter_symbol, parameter_range; kwargs...)

    Solves the steady states of a system for every parameter values in a given range of values. Returns a bifurcation grid type object (which can be plotted).

## args
 - reaction_network: a reaction network.
 - solver (optional): a subtype of AbstractSteadyStateSolver that specifies how to solve the bifurcation diagram. Default is HCSteadyStateSolver (currently only solver avaiable).
 - parameter_values: a vector which specifies the parameter values for which the steady states are to be found. Not required in case the reaction network lacks steady states.
 - parameter_symbol: the parameter which is to be varried.
 - parameter_range: a range with all the values of the parameter for which steady states are to be solved.

## kwargs
 - potential arguments for steady state solvers (however, only current solver does not have any such arguments).
"""
function bifurcation_grid(rn::DiffEqBase.AbstractReactionNetwork, args...)
    return bifurcation_grid(rn, HCSteadyStateSolver(), args...)
end
function bifurcation_grid(rn::DiffEqBase.AbstractReactionNetwork, solver::AbstractSteadyStateSolver, p::Vector{Float64}, param::Symbol, range::AbstractRange)
    grid_points = Vector{Union{BifurcationPoint,Nothing}}(fill(nothing,length(range)))
    for i = 1:length(range)
        p_i=copy(p); p_i[rn.params_to_ints[param]]=range[i];
        sol = steady_states(rn, solver, p_i)
        jac_eigenvals = get_jac_eigenvals(map(s->ComplexF64.(s),sol),param,fill(range[i],length(sol)),rn,p)
        grid_points[i] = BifurcationPoint(sol,jac_eigenvals,stability_type.(jac_eigenvals,number_of_constraints(rn)))
    end
    return BifurcationGrid(param,range,Vector{BifurcationPoint}(grid_points),length(range))
end

#Generates a 2d grid of bifurcation points, using a given steady state method.
"""
    bifurcation_grid_2d(reaction_network, [solver], parameter_values, parameter_symbol1, parameter_range1, parameter_symbol2, parameter_range2; kwargs...)

    Ad bifurcation grid, but varries two different parameters. Returning a 2d grid with steady states in each grid point.

## args
 - reaction_network: a reaction network.
 - solver (optional): a subtype of AbstractSteadyStateSolver that specifies how to solve the bifurcation diagram. Default is HCSteadyStateSolver (currently only solver avaiable).
 - parameter_values: a vector which specifies the parameter values for which the steady states are to be found. Not required in case the reaction network lacks steady states.
 - parameter_symbo1l: the first parameter which is to be varried.
 - parameter_range1: a range with all the values of the first parameter for which steady states are to be solved.
 - parameter_symbol2: the second parameter which is to be varried.
 - parameter_range2: a range with all the values of the second parameter for which steady states are to be solved.

## kwargs
 - potential arguments for steady state solvers (however, only current solver does not have any such arguments).
"""
function bifurcation_grid_2d(rn::DiffEqBase.AbstractReactionNetwork, args...)
    return bifurcation_grid_2d(rn, HCSteadyStateSolver(), args...)
end
function bifurcation_grid_2d(rn::DiffEqBase.AbstractReactionNetwork, solver::AbstractSteadyStateSolver, p::Vector{Float64}, param1::Symbol, range1::AbstractRange, param2::Symbol, range2::AbstractRange)
    grid_points = Matrix{Union{BifurcationPoint,Nothing}}(fill(nothing,length(range1),length(range2)))
    for i = 1:length(range1), j = 1:length(range2)
        p_ij = copy(p);
        p_ij[rn.params_to_ints[param1]] = range1[i];
        p_ij[rn.params_to_ints[param2]] = range2[j];
        sol = steady_states(rn, solver, p_ij)
        jac_eigenvals = get_jac_eigenvals(map(s->ComplexF64.(s),sol),param1,fill(range1[i],length(sol)),rn,p)
        grid_points[i,j] = BifurcationPoint(sol,jac_eigenvals,stability_type.(jac_eigenvals,number_of_constraints(rn)))
    end
    return BifurcationGrid2D(param1,range1,param2,range2,Matrix{BifurcationPoint}(grid_points),(length(range1),length(range2)))
end

#Generates a grid of bifurcation diagram.
"""
    bifurcation_grid_diagram(reaction_network, [solver], parameter_values, parameter_symbol1, parameter_range1, parameter_symbol2, parameter_range2; kwargs...)

    Varries a first parameter over a range of discrete values, for each values makes a bifurcation diagram over a second, continious range. Allows to vsiualise changes ins steady states over two parameters.

## args
- reaction_network: a reaction network.
- solver (optional): a subtype of AbstractBifurcationSolver that specifies how to solve the bifurcation diagram. Default is HCBifurcationSolver.
- parameter_values: a vector which specifies the point in parameter space around which the bifurcation diagram is evaluated.
- parameter_symbol1: the first parameter which is varied, given as a symbol.
- parameter_range1: the range over which the first parameter is varied (this is a range of discrete numbers).
- parameter_symbol2: the second parameter which is varied, given as a symbol.
- parameter_rang21: the range over which the second parameter is varied (this is a tuple of two numbers).

## kwargs
- stepsize = (parameter_range[2] - parameter_range[1])/200: The maximum distance between two parameter values for which steady states are calculated. Smaller values will give a smoother diagram.
- Δp = (parameter_range[2] - parameter_range[1])/200: The distance to jump after finding a bifurcation. After discovering a bifurcation, the solver jumps ahead a distance `Δp` and starts back-tracking. If this distance is too small, the method may error or cause visible artifiacts in the bifurcation diagram. If it is too large, then you might jump over another bifurcation and it will be missed (this should also be farily obvious in a plot). Only applicable for the HCBifurcationSolver.
- Δu = 1e-7: (alternatively reduced to a tenth of the minimum value of any initial solutions) used by algorithm to determine if two points are identical. Set this to less than the expected minimum distance between two different solutions.
"""
function bifurcation_grid_diagram(rn::DiffEqBase.AbstractReactionNetwork, args...)
    bifurcation_grid_diagram(rn, HCBifurcationSolver(), args...)
end
function bifurcation_grid_diagram(rn::DiffEqBase.AbstractReactionNetwork, solver::AbstractBifurcationSolver, p::Vector{Float64}, param1::Symbol, range1::AbstractRange, param2::Symbol, range2::Tuple{Float64,Float64}; stepsize=(range2[2]-range2[1])/200.::Float64)
    diagram_grid = Vector{Union{BifurcationDiagram,Nothing}}(fill(nothing,length(range1)))
    for i = 1:length(range1)
        p_i=copy(p); p_i[rn.params_to_ints[param1]]=range1[i];
        diagram_grid[i] = bifurcations(rn, solver, p_i, param2, range2; stepsize=stepsize)
    end
    return BifurcationGridDiagram(param1,range1,param2,range2,Vector{BifurcationDiagram}(diagram_grid),length(range1))
end

### Functions required for generating bifurcation diagrams. These should return a vector of bifurcation paths. ###

### Homotopy continuation based method for making bifurcation diagrams. Only works for polynomial systems and contains some heuretics for dealing with certain bifurcation points. ###

#--- Function for the simple HC bifurcation solver ---#
#Simple bifurcation solver which tracks the path from both sides. Might have problem with some diagrams with several bifurcation points.

function solve_bifurcation(
        rn::DiffEqBase.AbstractReactionNetwork,
        args...;
        kwargs...,
    )
    return solve_bifurcation(BifurcationSolver(), rn, args...; kwargs...)
end

function solve_bifurcation(
        ::SimpleHCBifurcationSolver,
        rn::DiffEqBase.AbstractReactionNetwork,
        p::Vector{Float64},
        param::Symbol,
        range::Tuple{Float64,Float64};
        stepsize=(range[2]-range[1])/200.::Float64,
        Δu=1e-4,
    )
    using_temp_poly = initialise_solver!(rn,p,rn.params_to_ints[param])
    p1 = copy(p); p1[rn.params_to_ints[param]] = range[1];
    p2 = copy(p); p2[rn.params_to_ints[param]] = range[2];
    sol1 = hc_solve_at(rn, p1)
    sol2 = hc_solve_at(rn, p2)
    Δu = min(Δu,sort(norm.([sol1...,sol2...]))[2]/10.)
    tracker1 = make_coretracker(rn,sol1,p1,p2,stepsize/(range[2]-range[1]))
    tracker2 = make_coretracker(rn,sol2,p2,p1,stepsize/(range[2]-range[1]))
    paths_complete   = Vector{NamedTuple{(:p, :u),Tuple{Vector{Float64},Vector{Vector{Complex{Float64}}}}}}()
    paths_incomplete = Vector{NamedTuple{(:p, :u),Tuple{Vector{Float64},Vector{Vector{Complex{Float64}}}}}}()
    for sol in sol1
        path = track_path(sol,tracker1,range...)
        if (tracker1.state.status == HomotopyContinuation.CoreTrackerStatus.success)
            remove_sol!(sol2,path.u[end],Δu)
            push!(paths_complete,path)
        else
            push!(paths_incomplete,path)
        end
    end
    for sol in sol2
        path = track_path(sol,tracker2,reverse(range)...)
        (tracker2.state.status == HomotopyContinuation.CoreTrackerStatus.success) && remove_path!(paths_incomplete,path.u[end],Δu)
        push!(paths_complete,path)
    end
    append!(paths_complete,paths_incomplete)
    using_temp_poly && finalise_solver!(rn)
    return bifurcation_paths(positive_real_projection.(paths_complete),param,range[1],range[2],rn,p)
end
#Used in the homotopy continuation bifurcation traces heuristics. If the path was traced succesfully to its end, checks if it corresponds to a solution in that end and removes that solution (it does not need to be attempted in the other direction).
function remove_sol!(results::Vector{Vector{ComplexF64}},path_fin::Vector{ComplexF64},Δu::Float64)
    for i = length(results):-1:1
        if maximum(abs.([imag.(path_fin.-results[i])..., real.(path_fin.-results[i])...]))<Δu
            deleteat!(results,i)
            return
        end
    end
end
#Used in the homotopy continuation bifurcation traces heuristics. Similar to the previous but removes an unfinished path.
function remove_path!(paths::Vector{NamedTuple{(:p, :u),Tuple{Vector{Float64},Vector{Vector{Complex{Float64}}}}}},path_fin::Vector{ComplexF64},Δu::Float64)
    for i = length(paths):-1:1
        if maximum(abs.([imag.(path_fin.-paths[i][2][1])..., real.(path_fin.-paths[i][2][1])...]))<Δu
            deleteat!(paths,i)
            return
        end
    end
end

#--- Function for the main HC bifurcation solver ---#
#Main bifurcation solver. Track paths from the first to the last parameter values. Detects whenever a solution is lost due to a complicated bifurcation. In this case it resumes the path tracking just after that bifurcation.
#Δp = how long step in parameter value to take after a bifrucation is encountered (should be smaller than the minimum expected distance between two bifurcations).
#Δu = distance between two solutions for them to be considered identical. Should larger than the distance between two solutions.

function solve_bifurcation(
        ::HCBifurcationSolver,
        rn::DiffEqBase.AbstractReactionNetwork,
        p::Vector{Float64},
        param::Symbol,
        range::Tuple{Float64,Float64};
        stepsize=(range[2]-range[1])/200.::Float64,
        Δp=(range[2]-range[1])/200.::Float64,
        Δu=1e-4::Float64,
    )
#= function HcBifurcationSolver(rn::DiffEqBase.AbstractReactionNetwork,p::Vector{Float64},param::Symbol,range::Tuple{Float64,Float64};stepsize=(range[2]-range[1])/200.::Float64,Δp=0.01::Float64,Δx=0.01::Float64) =#
    using_temp_poly = initialise_solver!(rn,p,rn.params_to_ints[param])
    parameters(p_val) = (p = copy(p); p[rn.params_to_ints[param]] = p_val; return p;)
    p_cur = range[1]; paths = Vector{NamedTuple{(:p, :u),Tuple{Vector{Float64},Vector{Vector{Complex{Float64}}}}}}();
    while p_cur < range[2]
        sol = hc_solve_at(rn, parameters(p_cur))
        Δu = min(Δu,sort(norm.(sol))[min(2,length(sol))]/10.)
        substract_sols!(sol,p_cur,paths,Δu)
        ct1 = make_coretracker(rn,sol,parameters(p_cur),parameters(range[1]),stepsize/(range[2]-range[1]))
        ct2 = make_coretracker(rn,sol,parameters(p_cur),parameters(range[2]),stepsize/(range[2]-range[1]))
        new_paths = track_path_two_ways(sol,ct1,ct2,p_cur,range)
        paths = combine_paths(paths,new_paths,p_cur,Δu,Δp)
        p_cur = minimum(map(path->path.p[end],new_paths)) + Δp
    end
    using_temp_poly && finalise_solver!(rn)
    return bifurcation_paths(positive_real_projection.(paths),param,range[1],range[2],rn,p)
end
#Tracks a given solution in both directions, and combines the paths together.
function track_path_two_ways(start_points,coretracker1,coretracker2,p_cur,range)
    paths = Vector{NamedTuple{(:p, :u),Tuple{Array{Float64,1},Array{Array{Complex{Float64},1},1}}}}();
    for sp in start_points
        (P1,U1) = track_path(sp,coretracker1,p_cur,range[1])
        (P2,U2) = track_path(sp,coretracker2,p_cur,range[2])
        push!(paths,(p=[reverse(P1)...,P2...],u=[reverse(U1)...,U2...]))
    end
    return paths
end
#For a given solution at a certain parameter value, remove those solutions already present in the given paths.
function substract_sols!(solutions,p,paths,Δu)
    for path in filter(path -> path.p[end]≥p, paths)
        idx = findfirst(path.p .> p)
        p_pre = path.p[idx-1]; p_post = path.p[idx];
        f1 = (p_post-p)/(p_post-p_pre); f2 = (p_post-p)/(p_post-p_pre);
        path_value_at_p = f1*path.u[idx-1]+f2*path.u[idx]
        closest_sol = argmin(map(sol->norm(sol-path_value_at_p),solutions))
        (norm(solutions[closest_sol]-path_value_at_p) < Δu) && deleteat!(solutions,closest_sol)
    end
end
#For two sets of paths, removes the paths in the old path vector which seems to be similar to those in the new (starts in the same point).
function combine_paths(paths_old,paths_new,p,Δu,Δp)
    output_paths = copy(paths_new)
    first_vals = map(path -> path.u[1],paths_new)
    last_vals = map(path -> path.u[end],paths_new)
    midpoint_vals = map(path -> path.u[findfirst(path.p .> p/2)],paths_new)
    for path in paths_old
        same_start_p = (path.p[1] .- map(p -> p.p[1], paths_new)) .< Δp
        same_end_p = (path.p[end] .- map(p -> p.p[end], paths_new)) .< Δp
        same_first_vals = norm.(map(fv->fv-path.u[1],first_vals)) .< Δu
        same_last_vals = norm.(map(lv->lv-path.u[end],last_vals)) .< Δu
        any(same_start_p .& same_end_p .& same_first_vals .& same_last_vals) && continue
        #any(norm.(map(fv->fv-path.u[1],first_vals)) .< Δu) && continue
        #(findfirst(path.p .> p/2)==nothing) || any(norm.(map(mv->mv-path.u[findfirst(path.p .> p/2)],midpoint_vals)) .< Δu) && continue
        push!(output_paths,path)
    end
    return output_paths
end

#--- Bifurcation functions used by both solvers ---#
#Makes a coretracker from the given information
function make_coretracker(rn,sol,p1,p2,Δt)
    return HomotopyContinuation.coretracker(get_equi_poly(rn), sol, parameters=get_polyvars(rn).p, p₁=p1, p₀=p2, max_step_size=Δt)
end
#Tracks a path for the given coretracker, returns the path.
function track_path(solution,coretracker,p_start,p_end)
    P = Vector{Float64}(); U = Vector{Vector{ComplexF64}}();
    for (u,t) in HomotopyContinuation.iterator(coretracker, solution)
        push!(P,t2p(t,p_start,p_end)); push!(U,u);
    end
    return (p=P,u=U)
end
#Coverts a given t value (as used by homotopy continuation, 0<=t<=1), to the corresponding parameter value.
function t2p(t,p_start,p_end)
    return p_start + (1-t)*(p_end-p_start)
end
#For a given path, filters away all biologically unplausible values (negative and imaginary).
function positive_real_projection(track_result::NamedTuple{(:p, :u),Tuple{Array{Float64,1},Array{Array{Complex{Float64},1},1}}})
    P = Vector{Float64}(); U = Vector{Vector{ComplexF64}}();
    for i = 1:length(track_result.p)
        if (minimum(real.(track_result.u[i]))>-0.0001)&&(maximum(abs.(imag.(track_result.u[i])))<0.0001)
            push!(P,track_result.p[i]); push!(U,real.(track_result.u[i]))
        end
    end
    return (p=P,u=U)
end
#Takes a set of paths as tracked by homotopy continuation and turns them into a vector of bifurcation paths.
function bifurcation_paths(paths::Vector{NamedTuple{(:p, :u),Tuple{Array{Float64,1},Array{Array{Complex{Float64},1},1}}}},param::Symbol,r1::Number,r2::Number,rn::DiffEqBase.AbstractReactionNetwork,p::Vector{Float64})
    bps = Vector{BifurcationPath}()
    for path in paths
        (length(path.p)==0) && continue
        (abs(1-path.p[1]/path.p[end])<0.0001) && continue
        jac_eigenvals = get_jac_eigenvals(path.u,param,path.p,rn,p)
        push!(bps,BifurcationPath(path.p, path.u, jac_eigenvals, stability_type.(jac_eigenvals,number_of_constraints(rn)), length(path.p)))
    end
    return bps
end

#---Functions relating to the detection of stability in bifurcation diagrams -###

#Generates a number between 0 and 3 corresponing to some stability type (0=unstable, 1=stable, 2=unstable with imaginary eigenvalues, 3=stable with imaginary eigenvalues).
function stability_type(eigenvalues::Vector{Any},nbr_of_constraints)
    stab_type = Int8(0)
    (sort(real(eigenvalues))[end-nbr_of_constraints]<1e-6)&&(stab_type+=1)
    any(imag(eigenvalues).>1e-6)&&(stab_type+=2)
    return stab_type
end
#For a set of fixed points, return a corresponding set with the eigen values of the jacobian at the specified points.
function get_jac_eigenvals(vals::Vector{Vector{ComplexF64}},param::Symbol,param_vals::Vector{Float64},rn::DiffEqBase.AbstractReactionNetwork,p::Vector{Float64})
    stabs = Vector{Vector{Any}}()
    (length(vals)==0)&&(return stabs)
    Jac_temp = zeros(length(vals[1]),length(vals[1]))
    for i = 1:length(vals)
        params_i = copy(p)
        params_i[rn.params_to_ints[param]] = param_vals[i]
        rn.jac(Jac_temp,vals[i],params_i,0.)
        push!(stabs,eigen(Jac_temp).values)
    end
    return stabs
end


text:
using Documenter, ImageUtils

makedocs(
    modules = [ImageUtils],
    format = :html,
    checkdocs = :exports,
    sitename = "ImageUtils.jl",
    pages = Any["index.md"]
)

deploydocs(
    repo = "github.com/tknopp/ImageUtils.jl.git",
)


text:
#=

StatsBase.coeftable(model::LMM) = error("coeftable is not defined for $(typeof(model)).")

StatsBase.deviance(model::LMM) = error("deviance is not defined for $(typeof(model)).")

StatsBase.islinear(model::LMM) = error("islinear is not defined for $(typeof(model)).")

StatsBase.nulldeviance(model::LMM) =
    error("nulldeviance is not defined for $(typeof(model)).")

StatsBase.nullloglikelihood(model::LMM) =
        error("nullloglikelihood is not defined for $(typeof(model)).")

StatsBase.score(model::LMM) = error("score is not defined for $(typeof(model)).")

=#
"""
    StatsBase.confint(lmm::LMM{T}; level::Real=0.95, ddf::Symbol = :satter) where T

Confidece interval for coefficients.

ddf = :satter/:residual

```math
CI_{U/L} = β ± SE * t_{ddf, 1-α/2}
```

See also: [`dof_satter`](@ref), [`dof_residual`](@ref)
"""
function StatsBase.confint(lmm::LMM{T}; level::Real=0.95, ddf::Symbol = :satter) where T
    isfitted(lmm) || error("Model not fitted")
    alpha = 1.0 - level
    if ddf == :satter
        ddfv = dof_satter(lmm)
    elseif ddf == :contain
        ddfv = dof_contain(lmm)
    elseif ddf == :residual
        ddfv = fill!(Vector{Float64}(undef, coefn(lmm)), dof_residual(lmm))
    end
    cis = Vector{Tuple{T, T}}(undef, coefn(lmm))
    for i = 1:coefn(lmm)
        #ERROR: ArgumentError: TDist: the condition ν > zero(ν) is not satisfied
        d = lmm.result.se[i] * quantile(TDist(ddfv[i]), 1.0 - alpha / 2)
        cis[i] = (lmm.result.beta[i] - d, lmm.result.beta[i] + d)
    end
    cis
end

#=
REML: n = total number of observation - number fixed effect parameters; d = number of covariance parameters
ML:, n = total number of observation; d = number of fixed effect parameters + number of covariance parameters.
=#
"""
    StatsBase.coef(lmm::LMM) = copy(lmm.result.beta)

Model coefficients (β).
"""
StatsBase.coef(lmm::LMM) = copy(lmm.result.beta)

"""
    StatsBase.coefnames(lmm::LMM) = StatsBase.coefnames(lmm.mf)

Coefficients names.
"""
StatsBase.coefnames(lmm::LMM) = StatsBase.coefnames(lmm.mf)

"""
    StatsBase.nobs(lmm::LMM)

Number of observations.
"""
function StatsBase.nobs(lmm::LMM)
    return length(lmm.data.yv)
end

"""
    StatsBase.dof_residual(lmm::LMM)

DOF residuals: N - rank(X).
"""
function StatsBase.dof_residual(lmm::LMM)
    nobs(lmm) - lmm.rankx
end

"""
    StatsBase.dof(lmm::LMM)

DOF.
"""
function StatsBase.dof(lmm::LMM)
    lmm.nfixed + lmm.covstr.tl
end

"""
    StatsBase.loglikelihood(lmm::LMM)

Return loglikelihood value.
"""
function StatsBase.loglikelihood(lmm::LMM)
    -lmm.result.reml/2
end

"""
    StatsBase.aic(lmm::LMM)

Akaike Information Criterion.
"""
function StatsBase.aic(lmm::LMM)
    l = loglikelihood(lmm)
    d = lmm.covstr.tl
    -2l + 2d
end

"""
    StatsBase.bic(lmm::LMM)

Bayesian information criterion.
"""
function StatsBase.bic(lmm::LMM)
    l = loglikelihood(lmm)
    d = lmm.covstr.tl
    n = nobs(lmm) - lmm.nfixed
    -2l + d * log(n)
end

"""
    StatsBase.aicc(lmm::LMM)

Corrected Akaike Information Criterion.
"""
function StatsBase.aicc(lmm::LMM)
    l = loglikelihood(lmm)
    d = lmm.covstr.tl
    n = nobs(lmm) - lmm.nfixed
    -2l + (2d * n) / (n - d - 1.0)
end

"""
    caic(lmm::LMM)

Conditional Akaike Information Criterion.
"""
function caic(lmm::LMM)
    l = loglikelihood(lmm)
    d = lmm.covstr.tl
    n = nobs(lmm) - lmm.nfixed
    -2l + d * (log(n) + 1.0)
end

"""
    StatsBase.isfitted(lmm::LMM)
"""
function StatsBase.isfitted(lmm::LMM)
    lmm.result.fit
end
"""
    StatsBase.vcov(lmm::LMM)

Variance-covariance matrix of β.
"""
StatsBase.vcov(lmm::LMM) = copy(lmm.result.c)
"""
    StatsBase.stderror(lmm::LMM)

Standard error
"""
StatsBase.stderror(lmm::LMM) = sqrt.(diag(vcov(lmm)))
"""
    StatsBase.modelmatrix(lmm::LMM)

Fixed effects matrix.
"""
StatsBase.modelmatrix(lmm::LMM) = lmm.data.xv
"""
    StatsBase.response(lmm::LMM)

Response vector.
"""
StatsBase.response(lmm::LMM) = lmm.data.yv

StatsBase.crossmodelmatrix(lmm::LMM) = (x = modelmatrix(lmm); Symmetric(x' * x))

#=
StatsBase.mss(model::LMM) = error("mss is not defined for $(typeof(model)).")

StatsBase.rss(model::LMM) = error("rss is not defined for $(typeof(model)).")

StatsBase.informationmatrix(model::LMM; expected::Bool = true) =
    error("informationmatrix is not defined for $(typeof(model)).")

StatsBase.weights(model::LMM) = error("weights is not defined for $(typeof(model)).")

function StatsBase.r2(model::LMM)
    Base.depwarn("The default r² method for linear models is deprecated. " *
                 "Packages should define their own methods.", :r2)
    mss(model) / deviance(model)
end

function StatsBase.r2(model::LMM, variant::Symbol)
    loglikbased = (:McFadden, :CoxSnell, :Nagelkerke)
    if variant in loglikbased
        ll = loglikelihood(model)
        ll0 = nullloglikelihood(model)
        if variant == :McFadden
            1 - ll/ll0
        elseif variant == :CoxSnell
            1 - exp(2 * (ll0 - ll) / nobs(model))
        elseif variant == :Nagelkerke
            (1 - exp(2 * (ll0 - ll) / nobs(model))) / (1 - exp(2 * ll0 / nobs(model)))
        end
    elseif variant == :devianceratio
        dev  = deviance(model)
        dev0 = nulldeviance(model)
        1 - dev/dev0
    else
        error("variant must be one of $(join(loglikbased, ", ")) or :devianceratio")
    end
end

const r² = r2

StatsBase.adjr2(model::LMM) = error("adjr2 is not defined for $(typeof(model)).")

function StatsBase.adjr2(model::LMM, variant::Symbol)
    k = dof(model)
    if variant == :McFadden
        ll = loglikelihood(model)
        ll0 = nullloglikelihood(model)
        1 - (ll - k)/ll0
    elseif variant == :devianceratio
        n = nobs(model)
        dev  = deviance(model)
        dev0 = nulldeviance(model)
        1 - (dev*(n-1))/(dev0*(n-k))
    else
        error("variant must be one of :McFadden or :devianceratio")
    end
end

const adjr² = adjr2

StatsBase.responsename(model::LMM) = error("responsename is not defined for $(typeof(model)).")

StatsBase.meanresponse(model::LMM) = error("meanresponse is not defined for $(typeof(model)).")

StatsBase.leverage(model::LMM) = error("leverage is not defined for $(typeof(model)).")

StatsBase.residuals(model::LMM) = error("residuals is not defined for $(typeof(model)).")

StatsBase.predict(model::LMM) = error("predict is not defined for $(typeof(model)).")

StatsBase.predict!(model::LMM) = error("predict! is not defined for $(typeof(model)).")

=#


text:
function conv1x1_grad_v(
    X::AbstractArray{T,N},
    ΔY::AbstractArray{T,N},
    C::Conv1x1;
    adjoint = false,
) where {T,N}

    # Reshape input
    n_in, batchsize = size(X)[N-1:N]
    v1 = C.v1.data
    v2 = C.v2.data
    v3 = C.v3.data
    k = length(v1)

    dv1 = InvertibleNetworks.cuzeros(X, k)
    dv2 = InvertibleNetworks.cuzeros(X, k)
    dv3 = InvertibleNetworks.cuzeros(X, k)

    V1 = v1 * v1' / (v1' * v1)
    V2 = v2 * v2' / (v2' * v2)
    V3 = v3 * v3' / (v3' * v3)

    dV1 = InvertibleNetworks.partial_derivative_outer(v1)
    dV2 = InvertibleNetworks.partial_derivative_outer(v2)
    dV3 = InvertibleNetworks.partial_derivative_outer(v3)

    ∂V1 = deepcopy(dV1)
    ∂V2 = deepcopy(dV2)
    ∂V3 = deepcopy(dV3)

    M1 = (I - 2 * (V2 + V3) + 4 * V2 * V3)
    M3 = (I - 2 * (V1 + V2) + 4 * V1 * V2)
    tmp = InvertibleNetworks.cuzeros(X, k, k)
    for i = 1:k
        # ∂V1
        mul!(tmp, ∂V1[i, :, :], M1)
        adjoint ? adjoint!(∂V1[i, :, :], tmp) : copyto!(∂V1[i, :, :], tmp)
        # ∂V2
        v2 = ∂V2[i, :, :]
        broadcast!(+, tmp, v2, 4 * V1 * v2 * V3 - 2 * (V1 * v2 + v2 * V3))
        adjoint ? adjoint!(∂V2[i, :, :], tmp) : copyto!(∂V2[i, :, :], tmp)
        # ∂V3
        mul!(tmp, M3, ∂V3[i, :, :])
        adjoint ? adjoint!(∂V3[i, :, :], tmp) : copyto!(∂V3[i, :, :], tmp)
    end

    prod_res = InvertibleNetworks.cuzeros(X, size(∂V1, 1), prod(size(X)[1:N-2]), n_in)
    inds = [i < N ? (:) : 1 for i = 1:N]
    for i = 1:batchsize
        inds[end] = i
        Xi = -2.0f0 * reshape(view(X, inds...), :, n_in)
        ΔYi = reshape(view(ΔY, inds...), :, n_in)
        broadcast!(
            +,
            dv1,
            dv1,
            InvertibleNetworks.custom_sum(
                InvertibleNetworks.mat_tens_i(prod_res, Xi, ∂V1, ΔYi),
                (3, 2),
            ),
        )
        broadcast!(
            +,
            dv2,
            dv2,
            InvertibleNetworks.custom_sum(
                InvertibleNetworks.mat_tens_i(prod_res, Xi, ∂V2, ΔYi),
                (3, 2),
            ),
        )
        broadcast!(
            +,
            dv3,
            dv3,
            InvertibleNetworks.custom_sum(
                InvertibleNetworks.mat_tens_i(prod_res, Xi, ∂V3, ΔYi),
                (3, 2),
            ),
        )
    end
    return dv1, dv2, dv3
end


text:
module LightGraphsIO

using LightGraphs
using Parsers

abstract type AbstractGraphFormat end
struct DOTFormat <: AbstractGraphFormat end
struct EdgeListFormat <: AbstractGraphFormat end
struct GMLFormat <: AbstractGraphFormat end
struct NETFormat <: AbstractGraphFormat end
struct GEXFFormat <: AbstractGraphFormat end
struct Graph6Format <: AbstractGraphFormat end
struct GraphMLFormat <: AbstractGraphFormat end
struct CDFFormat <: AbstractGraphFormat end
struct LGCompressedFormat <: AbstractGraphFormat end

function write_file(io, file_name)
    open(file_name, "w") do f
        write(f, String(take!(io)))
    end
end

################ DOT

"""
    load_graph(file, ::Type{DOTFormat}) -> Union{SimpleGraph, SimpleDiGraph}

    Parses and produces a `file` written in the `DOTFormat`.
    CAVEAT:
    * Parses simple directed and undirected graphs.
    * Assumes 1 graph per file (for now)
    * Please sponsor/PR/file issues for more parsing support.

    # Arguments
    # - `file`: the file were the graph is written
    # - `DOTFormat`: DOT graph format
    #
    # # Returns
    # - `g`: a simple un/directed LightGraphs.jl graph
"""
function load_graph(file, ::Type{DOTFormat})
    if occursin("digraph", readline(file))
        _load_graph(SimpleDiGraph(0), DOTFormat, file)
    else
        _load_graph(SimpleGraph(0), DOTFormat, file)
    end
    return nothing
end

### Jacob Says:
# I'd be happy to talk through things. It sounds like Parsers.jl (https://github.com/JuliaData/Parsers.jl)
# could be helpful, though I know it can be a bit intimidating to dive into.
#
# Basically you can create a Parsers.Options with the config you want (e.g. delim='\t', comment="#" )
# then you would call Parsers.xparse , which is the internal (yet part of the public API)
# method that Parsers.parse and Parsers.tryparse call.
#
# For best performance, I'd recommend Mmap.mmap(file) you get a Vector{UInt8}, then pass that to xparse .
# You can see an example of how CSV.jl uses xparse directly for parsing column names from a file:
# https://github.com/JuliaData/CSV.jl/blob/master/src/detection.jl#L231.
#
# Anyway, that's probably enough to get started, but happy to chat if something doesn't make sense.

### Jacob also says
# One last thing that's probably helpful; xparse won't throw any errors (by design),
# but you get back a code return value.
#
# It's a bitmask of basically everything that happened while parsing and whether it ran into problems or whatnot.
#
# The documented codes are here:
# https://github.com/JuliaData/Parsers.jl/blob/ab5ef1bbdc81fe8ee979a5b287ea065d991ba0ce/src/utils.jl#L44.
#
# In particular, you can check Parsers.newline(code) to see if a newline was encountered while parsing
# to know if you hit the end of the line.

function _load_graph(g::SimpleGraph, ::Type{DOTFormat}, file)
    # 1. Create a Parsers.Options with the right presets
    # 2. Mmap the file (Mmap.mmap(file)) to get a Vector{UInt8}
    # 3. Pass that to Parsers.xparse
    # 4. Throw errors appropriately with the returned bytecode value
    return nothing
end

function _load_graph(g::SimpleDiGraph, ::Type{DOTFormat}, file)
    return nothing
end

function write_graph(g::SimpleGraph, ::Type{DOTFormat}, file_name="graph.dot", graph_name="")
    io = IOBuffer()
    println(io, "graph ", graph_name == "" ? '{' : string(graph_name, " {"))
    for e in edges(g)
        source = string(src(e))
        dest = string(dst(e))
        println(io, '\t', source, " -- ", dest)
    end
    write(io, '}')
    write_file(io, file_name)
end

function write_graph(g::SimpleDiGraph, ::Type{DOTFormat}, file_name="digraph.dot", graph_name="")
    io = IOBuffer()
    println(io, "digraph ", graph_name == "" ? '{' : string(graph_name, " {"))
    for v in vertices(g)
        out_nbrs = outneighbors(g, v)
        length(out_nbrs) == 0 && continue
        println(io, '\t', string(v), " -> {", join(out_nbrs, ','), '}')
    end
    write(io, '}')
    write_file(io, file_name)
    return nothing
end

################ EdgeList
function load_graph(file, ::Type{EdgeListFormat})
    2
end

function write_graph(g, ::Type{EdgeListFormat}, file_name="graph.edgelist")
    io = IOBuffer()
    for e in edges(g)
        println(io, join((src(e), dst(e)), ','))
    end
    write_file(io, filename)
    return nothing
end

################ GML
function load_graph(file, ::Type{GMLFormat})
    3
end

function write_graph(g, ::Type{GMLFormat}, file_name="graph.gml", graph_name="")
    io = IOBuffer()
    println(io, "graph\n[")
    length(graph_name) > 0 && println(io, "label ", graph_name)
    is_directed(g) && println(io, "directed 1")
    for i in 1:nv(g)
        println(io, "\tnode")
        println(io, "\t[")
        println(io, "\t\tid ",i)
        println(io, "\t]")
    end
    for e in LightGraphs.edges(g)
        s, t = Tuple(e)
        println(io, "\tedge")
        println(io, "\t[")
        println(io, "\t\tsource ", s)
        println(io, "\t\ttarget ", t)
        println(io, "\t]")
    end
    write(io, ']')
    write_file(io, file_name)
    return nothing
end

################ NET
function load_graph(file, ::Type{NETFormat})
    4
end

function write_graph(g, ::Type{NETFormat}, file_name="graph.net")
    io = IOBuffer()
    println(io, "*Vertices ", nv(g))
    println(io, is_directed(g) ? "*Arcs" : "*Edges")
    for e in edges(g)
        println(io, src(e), ' ', dst(e))
    end
    write_file(io, (is_directed(g) ? "di" : "") * file_name)
    return nothing
end

############### GEXF
function load_graph(file, ::Type{GEXFFormat})
    5
end

function write_graph(g, ::Type{GEXFFormat})

    5 #LOL not rewriting EzXML soon
end

############### Graph6
function load_graph(file, ::Type{Graph6Format})
    6
end

function write_graph(g, ::Type{Graph6Format}, file_name="graph.g6")
    A = adjacency_matrix(g, Bool)
    n = nv(g)
    nbits = div(n * (n - 1), 2)
    x = BitVector(undef, nbits)
    ind = 0
    for col in 2:n, row in 1:(col - 1)
        x[ind] = A[row, col]
        ind += 1
    end
    #join(">>graph6<<",string(
    return nothing
end

############### GraphML
function load_graph(file, ::Type{GraphMLFormat})
    7
end

function write_graph(g, ::Type{GraphMLFormat})

end

############### CDF
function load_graph(file, ::Type{CDFFormat})
    8
end

function write_graph(g, ::Type{CDFFormat})
    8
end

############### LGCompressed
function load_graph(file, ::Type{LGCompressedFormat})
    9
end

function write_graph(g, ::Type{LGCompressedFormat}, file_name="graph.jgz")

end

# Exports!
export AbstractGraphFormat
export DOTFormat
export EdgeListFormat, GMLFormat
export NETFormat,      GEXFFormat
export Graph6Format,   GraphMLFormat
export CDFFormat,      LGCompressedF
export load_graph, write_graph

end # module


text:
#=

  Murder mystery

  From
  Andy Gordpn:
  "Reverend Bayes, meet Countess Lovelace: Probabilistic Programming for Machine Learning"
  https://channel9.msdn.com/Events/Lang-NEXT/Lang-NEXT-2012/Reverend-Bayes-meet-Countess-Lovelace-Probabilistic-Programming-for-Machine-Learning
  Around @14:00

  """
  Miss Scarlett dunnit 30%. Col Mustard dunnit 70%.
  Scarlett uses gun 3%, uses pipe 97%.
  Mustard uses gun 80%, uses pipe 20%.
  ...
  We found a gun at the Scene.
  What is the probability that Scarlett dunnit?
  """

  * Original question: we found a gun, what is the
    probability that Scarlett did it?

   It's not very likely that Scarlett did it!

   Summary Statistics
   parameters      mean       std   naive_se      mcse         ess      rhat   ess_per_sec
      Symbol   Float64   Float64    Float64   Float64     Float64   Float64       Float64

    scarlett    0.0142    0.1183     0.0012    0.0013   7957.8849    1.0001      605.3925
     mustard    0.7024    0.4572     0.0046    0.0057   7373.2272    1.0002      560.9150
     withGun    1.0000    0.0000     0.0000    0.0000         NaN       NaN           NaN
    withPipe    0.2031    0.4023     0.0040    0.0038   8323.5860    0.9999      633.2131

  * If we had found a pipe instead:

    It's a little more probable that Mustard did it. But Scarlett is still under
    investigation for it.

   Summary Statistics
   parameters      mean       std   naive_se      mcse         ess      rhat   ess_per_sec
      Symbol   Float64   Float64    Float64   Float64     Float64   Float64       Float64

    scarlett    0.6620    0.4731     0.0047    0.0059   6800.0219    1.0003      441.2161
     mustard    0.6979    0.4592     0.0046    0.0054   6727.5710    0.9999      436.5151
     withGun    0.2908    0.4542     0.0045    0.0060   6720.3716    1.0002      436.0480
    withPipe    1.0000    0.0000     0.0000    0.0000         NaN       NaN           NaN

  * Here's the priors (after we just heard about the murder)
    Mustard probably did it and it's more likely that it was with a Gun than a knife.

   Summary Statistics
   parameters      mean       std   naive_se      mcse         ess      rhat   ess_per_sec
      Symbol   Float64   Float64    Float64   Float64     Float64   Float64       Float64

    scarlett    0.2923    0.4548     0.0045    0.0049   9712.4171    1.0002    33035.4322
     mustard    0.7020    0.4574     0.0046    0.0052   8739.3445    1.0002    29725.6616
     withGun    0.5762    0.4942     0.0049    0.0052   9747.4247    1.0000    33154.5059
    withPipe    0.4247    0.4943     0.0049    0.0046   9917.5952    0.9999    33733.3169

  cf ~/blog/murder_mystery.blog
     ~/psi/murder_mystery.psi
     ~/webppl/murder_mystery.wppl

=#

using Turing, StatsPlots
include("jl_utils.jl")

@model function murder_mystery()
    scarlett ~ flip(0.30)
    mustard ~ flip(0.70)

    withGun ~ scarlett ? flip(0.03) : flip(0.80)
    withPipe ~ scarlett ? flip(0.97) : flip(0.20)

    true ~ Dirac(withGun == true)
    # true ~ Dirac(withPipe == true)
    # true ~ Dirac(withGun < withPipe)
    # true ~ Dirac(withGun > withPipe)

end

model = murder_mystery()

# chns = sample(model, Prior(), 10_000)
chns = sample(model, MH(), 10_000)
# chns = sample(model, PG(15), 10_000)
# chns = sample(model, SMC(), 10_000)
# chns = sample(model, IS(), 10_000)

display(chns)
# display(plot(chns))

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
make_nondeterministicthreading(; basesize = 1, ntasks = Threads.nthreads(), kwargs...) =
    (NondeterministicThreading(; basesize = basesize, ntasks = ntasks), (; kwargs...))

function Transducers.transduce(xf, rf, init, xs, ex::NondeterministicEx)
    ndt, kwargs = make_nondeterministicthreading(; ex.kwargs...)
    xf0, xs0 = extract_transducer(xs)
    # TODO: don't assume all transducers are parallelizable
    return transduce(xf ∘ xf0 ∘ ndt, rf, init, xs0; kwargs...)
end

# TODO: move ThreadedNondeterministicEx to this file


text:
function get_nbins(alg, data)

    n = length(data)

    if alg == :sqrt
        nbins = ceil(Int, sqrt(n))
    elseif alg == :sturges
        nbins = ceil(Int, log(2,n)) + 1
    elseif alg == :rice
        nbins = ceil(Int, 2*cbrt(n))
    elseif alg == :doane
        g₁ = moment(data, 3)
        σ = sqrt((6*(n-2))/((n+1)*(n+3)))
        nbins = ceil(Int, 1 + log(2,n) + log(2, 1+abs(g₁)/σ))
    elseif alg == :scott
        σ = std(data)
        binwidth = 3.5*σ/cbrt(n)
        lo, hi = extrema(data)
        nbins = ceil(Int, (hi - lo)/binwidth)
    elseif alg == :fd
        binwidth = 2*iqr(data)/cbrt(n)
        lo, hi = extrema(data)
        nbins = ceil((hi - lo)/binwidth)
    else
        binwidth = 2*iqr(data)/cbrt(n)
        lo, hi = extrema(data)
        nbins_fd = ceil(Int, (hi - lo)/binwidth)
        nbins_sturges = ceil(Int, log(2,n)) + 1
        nbins = max(nbins_fd, nbins_sturges)
    end

    return(round(nbins))
end

# M: data
# ndim: number of dimensions
# k: number of eigenvalues to keep
# d: discretizer used per principal components: :bb for bayesian blocks :uw for uniform width :km for kmeans
# bintype: type of binning, only used if disc=:uw, select :manual for manual setting of number of bins
# bins: number of bins, only used if bintype=:manual
function CPD(M, ndim::Int64, npoints::Int64, k::Int64, d=:km, bintype=:manual, bins=Int(round(ndim)),processing=:none)
    results = DataFrame()

    #calculate correlation matrix
    C=cor(M');

    #select principal components
    K = size(C, 1)
    ef=eigen(C)
    ef.values[K-k+1:K]
    S=ef.vectors[:, K-k+1:K]
    S=Matrix(S)

    #project onto principal components
    D=S'*M

    #determine cutpoints
    B = []
    for m in 1:k
        if d == :km
            R = kmeans(copy(convert(Array,D[m,:])'),bins)
            cutpoints=zeros(bins+1)
            cutpoints[1]=minimum(D[m,:])
            cutpoints[bins+1]=maximum(D[m,:])
            centers = sort(R.centers,dims=2)
            [cutpoints[i] = (centers[i-1]+centers[i])/2 for i in 2:bins]
            push!(B,cutpoints)
        else
            if d == :bb
                disc = LinearDiscretizer(binedges(DiscretizeBayesianBlocks(), D[m,:]))
            else
                if bintype==:manual
                    nbins=bins
                else
                    nbins=get_nbins(bintype, D[m,:])
                end

                disc = LinearDiscretizer(binedges(DiscretizeUniformWidth(nbins), D[m,:]))
            end
            cutpoints = zeros(nlabels(disc)+1)
            cutpoints[1] = extrema(disc,1)[1]
            #cutpoints = zeros(nlabels(disc))
            [cutpoints[i+1]=extrema(disc,i)[2] for i in 1:nlabels(disc)]
            push!(B,cutpoints)
        end
    end

    #reprojection the cutpoints
    cp_proj = []
    for i in B
        for j in i
            for l in 1:k
                push!(cp_proj, j*S[:,l])
            end
        end
    end

    projections = zeros(ndim,size(cp_proj)[1])

    for i in 1:ndim
        for j in 1:size(cp_proj)[1]
            projections[i,j]=cp_proj[j][i]
        end
        lindisc = LinearDiscretizer(sort(projections[i,:]))
        results[!,i] = [encode(lindisc,one) for one in M[i,:]]
    end

    if processing==:Integer
        results = postprocessing_integer(results,ndim,npoints)
    elseif processing==:Cutpoints
        results = postprocessing_cutpoints(results,ndim,npoints,projections)
    end

    return(results)
end

function CPD_clustered(M, ndim::Int64, npoints::Int64, c::Int64, p::Float64, d=:km)
    results = zeros(size(M))
    SPC = 1/2*(I-corspearman(M'))
    clust = hclust(SPC,linkage=:ward_presquared)
    clustered_data = cutree(clust;k=c)

    for i in 1:c
        selector = clustered_data .== i
        data = M[selector,:]
        k = Int(ceil(size(data)[1]/p))
        bins = get_nbins(:doane, data)
        res = CPD(data,size(data)[1],npoints,k,d,:manual,bins)
        #print("\n")
        #print(typeof(res))
        results[selector,:] .= Array(res)'
    end

    return(DataFrame(results))
end


text:
# Preamble ########################################################################################
using Pkg
Pkg.activate(Base.current_project())

cd("..") # Go to the parent directory
###################################################################################################
# Packages needed #################################################################################
using PolyGibbs
using DelimitedFiles
using ImageCore
using ImageIO
using ImageMagick
using FileIO
using PyPlot
###################################################################################################
# Angle of rotation ###############################################################################

const θ = π/4

###################################################################################################
# Retrieve of data from image #####################################################################

@info "Retrieving the image and extracting the two-dimensional field"
input = Float64.(channelview(Gray.(load("data/origin/symm/symm50.png"))))

###################################################################################################
# Identification of the grayscale field dimensions ################################################

j = (size(input)[1] - 1)/2;

N = Integer(2*j + 1);

@info "The size of the field is $N"

###################################################################################################
# Rotation of the image ###########################################################################

@info "Performing rotation at θ = $θ"

output = Emodes(input, θ)

###################################################################################################
# Saving results ##################################################################################

@info "Writting the results"
writedlm("data/target/symm/symm50_rpi4.dat", real(output))

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

begin
	fig = figure()
	ax = gca()
	ax.imshow(real(output), cmap = "gray")
	tight_layout()
	savefig("data/target/symm/symm50_rpi4.png", dpi = 300)
end

max_stars_repo_name: rurz/PolyGibbs

max_issues_repo_name: rurz/PolyGibbs

max_forks_repo_name: rurz/PolyGibbs


text:
export optimal_control

"""
	optimal_control(CP::ControlProcess, μ0::Dict, d::Int, trange::AbstractVector{<:Real}, solver)

returns a **lower** bound on the objective value of the (stochastic) optimal
control problem specified by CP. `μ0` encodes information about the distribution of
the initial state of the process; specifically, `μ0` maps a given monomial to the
corresponding moment of the initial distribution. trange refers to an *ordered*
set of time points discretizing the control horizon. trange[end] should coincide
with the end of the control horizon, i.e., `trange[end] = Inf` in case of an infinite
horizon problem. The bound is computed via a SOS program of degree d solved with
an appropriate method given by solver.

The bound can be tightened by populating trange or increasing `d`.
"""
function optimal_control(CP::ControlProcess, μ0::Dict, d::Int, trange::AbstractVector{<:Real}, solver, P::Partition = trivial_partition(CP.MP.X))
	if trange[1] == 0
		trange = trange[2:end]
	end
	if isinf(trange[end])
		trange = trange[1:end-1]
	end
	if isinf(CP.T)
		model, w = infinite_horizon_control(CP, μ0, d, trange, solver, P)
	else
		model, w = finite_horizon_control(CP, μ0, d, trange, solver, P)
	end
	optimize!(model)
	return Bound(objective_value(model), model, P, dual_poly(w, CP.MP.iv, trange))
end

## Finite horizon control problems
function finite_horizon_control(CP::ControlProcess, μ0::Dict, d::Int, trange::AbstractVector{<:Real}, solver, P::Partition = trivial_partition(CP.MP.X))
	if isa(CP.Objective, LagrangeMayer)
		model, w = finite_horizon_LM(CP, μ0, d, trange, solver, P)
	elseif isa(CP.Objective, ExitProbability)
		model, w = finite_horizon_EP(CP, μ0, d, trange, solver, P)
	elseif isa(CP.Objective, TerminalSetProbability)
		model, w = finite_horizon_TP(CP, μ0, d, trange, solver, P)
	else
		error("Objective function type not supported")
	end
	return model, w
end

## Lagrange Mayer
function finite_horizon_LM(CP::ControlProcess, μ0::Dict, d::Int, trange::AbstractVector{<:Real}, solver, P::Partition = trivial_partition(CP.MP.X))
    MP = CP.MP
    t = MP.iv
    nT = length(trange)
    Δt = vcat(trange[1], [trange[i] - trange[i-1] for i in 2:nT])
    T = @set(t >= 0 && t <= 1)

    model = SOSModel(solver)
    w = Dict((k, v) => (props(P.graph, v)[:cell] isa Singleton ?
                    @variable(model, [1], Poly(monomials(t, 0:d)))[1] :
                    @variable(model, [1], Poly(monomials(sort(vcat(MP.x, t), rev = true), 0:d)))[1]) for v in vertices(P.graph), k in 1:nT)

    for v in vertices(P.graph)
        add_dynamics_constraints!(model, MP, v, P, props(P.graph, v)[:cell], intersect(T, CP.U), Δt, w, - CP.Objective.l)
    end

    for v in vertices(P.graph)
        add_transversality_constraints!(model, MP, props(P.graph, v)[:cell], w[nT, v], CP.Objective.m, v)
    end

    for e in edges(P.graph)
        add_coupling_constraints!(model, MP, e, P, T, Δt, w)
    end
	obj = expectation(μ0[first(keys(μ0))] isa Dict ? [subs(w[1, v], t => 0) for v in vertices(P.graph)] : subs(w[1, 1], t => 0), μ0)

	if !isempty(CP.PathChanceConstraints)
		# TODO: improve formulation to take advantage of partition
		for PC in CP.PathChanceConstraints
			s_pc = @variable(model)
			w_pc = @variable(model, [1:nT], Poly(monomials(sort(sort(vcat(MP.x, t), rev = true),rev=true), 0:d)))
			add_path_chance_constraint!(model, MP, PC, intersect(T, CP.U), Δt, w_pc, s_pc)
			obj += s_pc*(1-PC.α)
		end
	end
	if !isempty(CP.TerminalChanceConstraints)
		for TC in CP.TerminalChanceConstraints
			s_pc = @variable(model)
			@constraint(model, s_p >= 0)
			for v in vertices(P.graph)
				add_terminal_chance_constraint!(model, MP, TC, v, P, props(P.graph, v)[:cell], w, - s_pc + CP.Objective.m)
			end
			obj += s_pc*(1-TC.α)
		end
	end
	@objective(model, Max, obj)
    return model, w
end

## Exit probability
function finite_horizon_EP(CP::ControlProcess, μ0::Dict, d::Int, trange::AbstractVector{<:Real}, solver, P::Partition = trivial_partition(CP.MP.X))
	nT = length(trange)
	MP = CP.MP
	t = MP.iv
	Δt = [(i == 1 ? trange[i] : trange[i] - trange[i-1]) for i in 1:nT]
	T = @set(t >= 0 && t <= 1)

	model = SOSModel(solver)
	w = Dict((k, v) => (props(P.graph, v)[:cell] isa Singleton ?
                    @variable(model, [1], Poly(monomials(t, 0:d)))[1] :
                    @variable(model, [1], Poly(monomials(sort(vcat(MP.x, t), rev = true), 0:d)))[1]) for v in vertices(P.graph), k in 1:nT)

	∂X = ∂(CP.Objective.X)
	for v in vertices(P.graph)
		add_dynamics_constraints!(model, MP, v, P, props(P.graph, v)[:cell], intersect(T, CP.U), Δt, w, 0)
		add_boundary_constraints!(model, MP, v, P, props(P.graph, v)[:cell], T, w, CP.Objective.X, ∂X)
	end

	for v in vertices(P.graph)
		add_transversality_constraints!(model, MP, props(P.graph, v)[:cell], w[nT,v], 0)
	end

	for e in edges(P.graph)
		add_coupling_constraints!(model, MP, e, P, T, Δt, w)
	end

	obj = expectation([subs(w[1,v], t => 0) for v in vertices(P.graph)], μ0)
	if !isempty(CP.PathChanceConstraints)
		# TODO: improve formulation to take advantage of partition
		for PC in CP.PathChanceConstraints
			s_pc = @variable(model)
			w_pc = @variable(model, [1:nT], Poly(monomials(sort(vcat(MP.x, t),rev=true), 0:d)))
			add_path_chance_constraint!(model, MP, PC, intersect(T, CP.U), Δt, w_pc, s_pc)
			obj += s_pc*(1-PC.α)
		end
	end
	if !isempty(CP.TerminalChanceConstraints)
		error("Terminal chance constraints currently not supported for Exit Probability problems!")
	end
	@objective(model, Max, obj)
	return model, w
end

# terminal probability
function finite_horizon_TP(CP::ControlProcess, μ0::Dict, d::Int, trange::AbstractVector{<:Real}, solver)
	nT = length(trange)
	MP = CP.MP
	t = MP.iv
	Δt = [(i == 1 ? trange[i] : trange[i] - trange[i-1]) for i in 1:nT]
	T = @set(t >= 0 && t <= 1)

	model = SOSModel(solver)
	w = Dict((k, v) => (props(P.graph, v)[:cell] isa Singleton ?
                    @variable(model, [1], Poly(monomials(t, 0:d)))[1] :
                    @variable(model, [1], Poly(monomials(sort(vcat(MP.x, t), rev = true), 0:d)))[1]) for v in vertices(P.graph), k in 1:nT)

	for v in vertices(P.graph)
		add_dynamics_constraints!(model, MP, v, P, props(P.graph, v)[:cell], intersect(T, CP.U), Δt, w, 0)
	end

	for v in vertices(P.graph)
		add_transversality_constraints!(model, MP, intersect(props(P.graph, v)[:cell], CP.Objective.X), w[nT, v], -1)
		add_transversality_constraints!(model, MP, props(P.graph, v)[:cell], w[nT, v], 0)
	end

	for e in edges(P.graph)
		add_coupling_constraints!(model, MP, e, P, T, Δt, w)
	end
	obj = expectation(μ0[first(keys(μ0))] isa Dict ? [subs(w[1, v], t => 0) for v in vertices(P.graph)] : subs(w[1, 1], t => 0), μ0)

	if !isempty(CP.PathChanceConstraints)
		# TODO: improve formulation to take advantage of partition
		for PC in CP.PathChanceConstraints
			s_pc = @variable(model)
			w_pc = @variable(model, [1:nT], Poly(monomials(sort(vcat(MP.x, t),rev=true), 0:d)))
			add_path_chance_constraint!(model, MP, PC, intersect(T, CP.U), Δt, w_pc, s_pc)
			obj += s_pc*(1-PC.α)
		end
	end
	if !isempty(CP.TerminalChanceConstraints)
		for TC in CP.TerminalChanceConstraints
			s_pc = @variable(model)
			@constraint(model, s_p >= 0)
			for v in vertices(P.graph)
				add_terminal_chance_constraint!(model, MP, TC, v, P, props(P.graph, v)[:cell], w, - s_pc + CP.Objective.m)
			end
			obj += s_pc*(1-TC.α)
		end
	end
	@objective(model, Max, obj)
	return model, w
end

## Discounted inf horizon control problems
function infinite_horizon_control(CP::ControlProcess, μ0::Dict, d::Int, trange::AbstractVector{<:Real}, solver)
	@assert isa(CP.Objective, LagrangeMayer) "Objective function type not supported"
	@assert CP.Objective.m == 0 "Mayer term must be 0 for infinite horizon problem"
	@assert isempty(CP.ChancePathConstraints) "Chance path constraints are not supported in infinite horizon problems"
	@assert isempty(CP.ChanceTerminalConstraints) "Chance terminal constraints are not supported in infinite horizon problems"
	return infinite_horizon_LM(CP, μ0, d, trange, solver)
end

function infinite_horizon_control(CP::ControlProcess, μ0::Dict, d::Int, trange::AbstractVector, solver, P::Partition = trivial_partition(CP.MP.X))
    MP = CP.MP
    t = MP.iv
    nT = length(trange)
    Δt = vcat(trange[1], [trange[i] - trange[i-1] for i in 2:nT])
    ρ = CP.discount_factor
	T = @set(t >= 0 && t <= 1)

    model = SOSModel(solver)
	w = Dict((k, v) => (props(P.graph, v)[:cell] isa Singleton ?
                    @variable(model, [1], Poly(monomials(t, 0:d)))[1] :
                    @variable(model, [1], Poly(monomials(sort(vcat(MP.x, t), rev = true), 0:d)))[1]) for v in vertices(P.graph), k in 1:nT)
	w∞ = Dict(v => (props(P.graph, v)[:cell] isa Singleton ?
					@variable(model) :
					@variable(model, [1], Poly(monomials(MP.x, 0:d)))[1]) for v in vertices(P.graph))

	for v in vertices(P.graph)
        add_dynamics_constraints!(model, MP, v, P, props(P.graph, v)[:cell], intersect(T, CP.U), Δt, w, - CP.Objective.l; ρ = ρ)
		add_dynamics_constraints!(model, MP, v, P, props(P.graph, v)[:cell], intersect(CP.U), [1], w∞, - CP.Objective.l; ρ = ρ)
		if props(P.graph, v)[:cell] isa Singleton
			@constraint(model, w∞[v] - w[nT, v](t => 1) >= 0)
		else
			@constraint(model, w∞[v] - subs(w[nT, v], t => 1) >= 0, domain = props(P.graph, v)[:cell])
		end
		add_transversality_constraints!(model, MP, props(P.graph, v)[:cell], w∞, 0)
    end

    for e in edges(P.graph)
        add_coupling_constraints!(model, MP, e, P, T, Δt, w)
		add_coupling_constraints!(model, MP, e, P, T, [1], w∞)
    end

    @objective(model, Max, expectation([subs(w[1,v], t => 0) for v in vertices(P.graph)], μ0))
    return model, merge(w, w∞)
end

#=
function optimal_control(CP::ControlProcess, μ0::Dict, d::Int, trange::AbstractVector{<:Real}, P::Partition, solver)
	if trange[1] == 0
		trange = trange[2:end]
	end
	if trange[end] == Inf
		trange = trange[1:end-1]
	end
	if isinf(CP.T)
        model = infinite_horizon_control(CP, μ0, d, trange, P, solver)
    else
        model = finite_horizon_control(CP, μ0, d, trange, P, solver)
    end
    optimize!(model)
    return objective_value(model), termination_status(model), MOI.get(model, MOI.SolveTime())
end

function finite_horizon_control(CP::ControlProcess, μ0::Dict, d::Int, trange::AbstractVector{<:Real}, P::Partition, solver)
    MP = CP.MP
    t = MP.iv
    nT = length(trange)
    Δt = vcat(trange[1], [trange[i] - trange[i-1] for i in 2:nT])
    T = @set(t >= 0 && t <= 1)

    model = SOSModel(solver)
    w = [Dict(v => (props(P.graph, v)[:cell] isa Singleton ?
                    @variable(model, [1], Poly(monomials(t, 0:d)))[1] :
                    @variable(model, [1], Poly(monomials(vcat(MP.x, t), 0:d)))[1]) for v in vertices(P.graph))
         for k in 1:nT]

    for v in vertices(P.graph)
        add_dynamics_constraints!(model, MP, v, P, props(P.graph, v)[:cell], intersect(T, CP.U), Δt, w, - CP.Objective.l)
    end

    for v in vertices(P.graph)
        add_transversality_constraints!(model, MP, props(P.graph, v)[:cell], w[nT][v], CP.Objective.m)
    end

    for e in edges(P.graph)
        add_coupling_constraints!(model, MP, e, P, T, Δt, w)
    end

    @objective(model, Max, expectation([subs(w[1][v], t => 0) for v in vertices(P.graph)], μ0))
    return model
end

function infinite_horizon_control(CP::ControlProcess, μ0::Dict, d::Int, trange::AbstractVector, p::Partition, solver)
	if trange[1] == 0
        trange = trange[2:end]
    end
    MP = CP.MP
    t = MP.iv
    nT = length(trange)
    Δt = vcat(trange[1], [trange[i] - trange[i-1] for i in 2:nT])
    ρ = CP.discount_factor
	T = @set(t >= 0 && t <= 1)

    model = SOSModel(solver)
	w = [Dict(v => (props(P.graph, v)[:cell] isa Singleton ?
					@variable(model, [1], Poly(monomials(t, 0:d)))[1] :
					@variable(model, [1], Poly(monomials(vcat(MP.x, t), 0:d)))[1]) for v in vertices(P.graph))
		 for k in 1:nT]
	w∞ = Dict(v => (props(P.graph, v)[:cell] isa Singleton ?
					@variable(model) :
					@variable(model, [1], Poly(monomials(MP.x, 0:d)))[1]) for v in vertices(P.graph))

	for v in vertices(P.graph)
        add_dynamics_constraints!(model, MP, v, P, props(P.graph, v)[:cell], intersect(T, CP.U), Δt, w, - CP.Objective.l; ρ = ρ)
		add_dynamics_constraints!(model, MP, v, P, props(P.graph, v)[:cell], intersect(CP.U), [1], w∞, - CP.Objective.l; ρ = ρ)
		if props(P.graph, v)[:cell] isa Singleton
			@constraint(model, w∞[v] - w[nT][v](t => 1) >= 0)
		else
			@constraint(model, w∞[v] - subs(w[nT][v], t => 1) >= 0, domain = props(P.graph, v)[:cell])
		end
		add_transversality_constraints!(model, MP, props(P.graph, v)[:cell], w∞, 0)
    end

    for e in edges(P.graph)
        add_coupling_constraints!(model, MP, e, P, T, Δt, w)
		add_coupling_constraints!(model, MP, e, P, T, [1], w∞)
    end

    @objective(model, Max, expectation(w[1], μ0, P))
    return model
end
=#


text:
# script for setting the PyJulia Python path to the venv folder

ENV["PYTHON"] = "/home/calvang/Projects/pyjulia-experiments/env/bin/python3"
# packages to rebuild
Pkg.build("PyCall")

max_stars_repo_name: calvang/pyjulia-experiments

max_issues_repo_name: calvang/pyjulia-experiments

max_forks_repo_name: calvang/pyjulia-experiments


text:
using Distributed
@everywhere using Mamba, LinearAlgebra

## Data
dyes = Dict{Symbol, Any}(
  :y =>
    [1545, 1440, 1440, 1520, 1580,
     1540, 1555, 1490, 1560, 1495,
     1595, 1550, 1605, 1510, 1560,
     1445, 1440, 1595, 1465, 1545,
     1595, 1630, 1515, 1635, 1625,
     1520, 1455, 1450, 1480, 1445],
  :batches => 6,
  :samples => 5
)

dyes[:batch] = vcat([fill(i, dyes[:samples]) for i in 1:dyes[:batches]]...)
dyes[:sample] = vcat(fill(collect(1:dyes[:samples]), dyes[:batches])...)

## Model Specification

model = Model(

  y = Stochastic(1,
    (mu, batch, s2_within) -> MvNormal(mu[batch], sqrt(s2_within)),
    false
  ),

  mu = Stochastic(1,
    (theta, batches, s2_between) -> Normal(theta, sqrt(s2_between))
  ),

  theta = Stochastic(
    () -> Normal(0, 1000)
  ),

  s2_within = Stochastic(
    () -> InverseGamma(0.001, 0.001)
  ),

  s2_between = Stochastic(
    () -> InverseGamma(0.001, 0.001)
  )

)

## Initial Values
inits = [
  Dict(:y => dyes[:y], :theta => 1500, :s2_within => 1, :s2_between => 1,
       :mu => fill(1500, dyes[:batches])),
  Dict(:y => dyes[:y], :theta => 3000, :s2_within => 10, :s2_between => 10,
       :mu => fill(3000, dyes[:batches]))
]

## Sampling Schemes
scheme = [NUTS([:mu, :theta]),
          Slice([:s2_within, :s2_between], 1000.0)]

scheme2 = [MALA(:theta, 50.0),
           MALA(:mu, 50.0, Matrix{Float64}(I, dyes[:batches], dyes[:batches])),
           Slice([:s2_within, :s2_between], 1000.0)]

scheme3 = [HMC(:theta, 10.0, 5),
           HMC(:mu, 10.0, 5, Matrix{Float64}(I, dyes[:batches], dyes[:batches])),
           Slice([:s2_within, :s2_between], 1000.0)]

scheme4 = [RWM(:theta, 50.0, proposal=Cosine),
           RWM(:mu, 50.0),
           Slice([:s2_within, :s2_between], 1000.0)]

## MCMC Simulations
setsamplers!(model, scheme)
sim = mcmc(model, dyes, inits, 10000, burnin=2500, thin=2, chains=2)
describe(sim)

setsamplers!(model, scheme2)
sim2 = mcmc(model, dyes, inits, 10000, burnin=2500, thin=2, chains=2)
describe(sim2)

setsamplers!(model, scheme3)
sim3 = mcmc(model, dyes, inits, 10000, burnin=2500, thin=2, chains=2)
describe(sim3)

setsamplers!(model, scheme4)
sim4 = mcmc(model, dyes, inits, 10000, burnin=2500, thin=2, chains=2)
describe(sim4)


text:
using Random, Distributions
using ProgressBars

# Does this work?
ENV["JULIA_DEBUG"] = "all"

const τ = 0.6 # time correlation coefficient
const t = 3 # number of frames wrong
const expectedrates = [
    0.585, # delay-critical user u
    2 # delay-tolerant user v
]
const ζ = 2 # path loss exponent
const a = [
    0.7, # power for user u
    0.3  # power for user v
]
const d = [
    2, # distance of user u
    1  # distance of user v
]

const SNR = 10 # SNR in decibels
const total_frames = 20000;

channel_distribution = Normal(0, √2 / 2)

successful_frames = [0, 0]

frame = 0
for frame ∈ ProgressBar(1:total_frames)
    # Step 1: Calculate g of channel for each user
    h = rand(channel_distribution, (2,2)) * [1, 1im]
    g = h ./ sqrt.(1 .+ d.^ζ)

    # Step 1.1: Sort channel gains
    if abs2(h[1]) > abs2(h[2])
        h[2], h[1] = h[1], h[2]
        g = h ./ sqrt.(1 .+ d.^ζ)
    end

    # Step 1.5: Apply time-correlated small fading coefficient
    ω = rand(channel_distribution, (2,2)) * [1, 1im]
    h = τ^t * h + √(1 - τ^(2t)) * ω
    g = h ./ sqrt.(1 .+ d.^ζ)

    # Step 2: Calculate AWGN
    σ = 10^(- SNR / 20)
    n = rand(Normal(0, σ), 2)

    # Step 3: Calculate reception SNR
    # Note: We are using the actual noise value n, instead of just the SNR-based σ
    SNRuu = a[1] * abs2(g[1]) / (a[2] * abs2(g[1]) + σ^2)
    SNRvv = a[2] * abs2(g[2]) / σ^2
    SNRvu = a[1] * abs2(g[2]) / (a[2] * abs2(g[2]) + σ^2)

    receiver_rates = [
        log2(1 + SNRuu),
        log2(1 + SNRvv),
        log2(1 + SNRvu),
    ]

    global successful_frames += [
        (receiver_rates[1] >= expectedrates[1]) ? 1 : 0, #uu
        (receiver_rates[2] >= expectedrates[2] && receiver_rates[3] >= expectedrates[1]) ? 1 : 0
        #vv, vu
    ]
end

@debug "Completed after $total_frames iterations"
@info "Packet success: $successful_frames / $total_frames"
@info "Outage rate u: $(1 - successful_frames[1] / total_frames)"
@info "Outage rate v: $(1 - successful_frames[2] / total_frames)"

# Calculate theoretical values
σ = 10^(- SNR / 20)
ρ = 1 / σ^2
λ = 1 ./ (1 .+ d.^ζ)
r = 2 .^ expectedrates .- 1
φ = max(r[1] / (a[1] - a[2] * r[1]) , r[2] / a[2])

# NOMA - Outdated CSI
PoutTheoretical = [
    1 - exp(- 2 * r[1] * σ^2 / (2 - τ^(2t)) / (a[1] - a[2] * r[1]) / λ[1]),
    1 - 2 * exp(- φ * σ^2 / λ[2]) + exp(- 2 * φ * σ^2 / (2 - τ^(2t)) / λ[2])
]

@info "Theor. Outage Rate u: $(PoutTheoretical[1])"
@info "Theor. Outage Rate v: $(PoutTheoretical[2])"