text:
using Distributed, DistributedArrays,
     PyPlot, LinearAlgebra, Formatting, Dates

import .AbstractOperator.Operator
import .AbstractOperator.Sounding
import .AbstractOperator.Operator
import .AbstractOperator.get_misfit
import .AbstractOperator.Sounding
import .AbstractOperator.makeoperator
import .AbstractOperator.make_tdgp_opt
import .AbstractOperator.getresidual

struct Tpointer
    fp   :: IOStream
    fstr :: String
end

mutable struct Chain
    pid           :: Int
    npidsperchain :: Int
    T             :: Float64
    misfit        :: Float64
end

function Chain(chainprocs::Array{Int, 1};
               Tmax          = 2.5,
               nchainsatone  = 1)

    @assert Tmax > 1
    nchains = length(chainprocs)
    npidsperchain = floor(Int, length(chainprocs)/nchains)
    @info "npidsperchain = $npidsperchain"
    T = 10.0.^range(0, stop = log10(Tmax), length = nchains-nchainsatone+1)
    append!(T, ones(nchainsatone-1))
    chains = Array{Chain, 1}(undef, nchains)

    pid_end = 0
    for ichain in 1:nchains
        pid_start      = pid_end + 1
        pid_end        = pid_start + npidsperchain - 1
        pids           = chainprocs[pid_start:pid_end]

        chains[ichain] = Chain(pids[1], npidsperchain, T[ichain], 0.0)
    end
    chains
end

function accept(current_misfit, new_misfit, stat, Temp, movetype)
    logalpha = (current_misfit[1] - new_misfit)/Temp
    accepted = true
    if log(rand()) < logalpha
        current_misfit[1] = new_misfit
        stat.accepted_moves[movetype] += 1
    else
        accepted = false
    end
    accepted
end

function mh_step!(m::ModelStat, mn::ModelNuisance,
    F::Operator, opt::OptionsStat, stat::Stats,
    Temp::Float64, movetype::Int, current_misfit::Array{Float64, 1})
    # for purely stat move and nuisance
    new_misfit = get_misfit(m, mn, opt, F)
    accepted = accept(current_misfit, new_misfit, stat, Temp, movetype)
    !accepted && undo_move!(movetype, m, opt)
end

function mh_step!(mns::ModelNonstat, m::ModelStat, mn::ModelNuisance,
    F::Operator, optns::OptionsNonstat, stat::Stats,
    Temp::Float64, movetype::Int, current_misfit::Array{Float64, 1})
    # for nonstat move and nuisance
    new_misfit = get_misfit(mns, mn, optns, F)
    accepted = accept(current_misfit, new_misfit, stat, Temp, movetype)
    !accepted && undo_move!(movetype, mns, optns, m)
end

function mh_step!(m::ModelStat, mns::ModelNonstat, mn::ModelNuisance,
    F::Operator, opt::OptionsStat, optns::OptionsNonstat,
    stat::Stats, Temp::Float64, movetype::Int, current_misfit::Array{Float64, 1})
    # for stat move updating nonstat and nuisance
    new_misfit = get_misfit(mns, mn, opt, F)
    accepted = accept(current_misfit, new_misfit, stat, Temp, movetype)
    !accepted && undo_move!(movetype, m, opt, mns, optns)
end

function mh_step!(mn::ModelNuisance, m::Model, F::Operator,
    optn::OptionsNuisance, statn::Stats,
    Temp::Float64, movetype::Int, current_misfit::Array{Float64,1})
    # this only gets called for moves on the nuisance chain
    # can be either nonstat or stat
    new_misfit = get_misfit(m, mn, optn, F)
    accepted = accept(current_misfit, new_misfit, statn, Temp, movetype)
    !accepted && undo_move!(mn)
end

function mh_step!(m::ModelStat, F::Operator, opt::OptionsStat, stat::Stats,
    Temp::Float64, movetype::Int, current_misfit::Array{Float64, 1})
    # for purely stat move
    new_misfit = get_misfit(m, opt, F)
    accepted = accept(current_misfit, new_misfit, stat, Temp, movetype)
    !accepted && undo_move!(movetype, m, opt)
end

function mh_step!(mns::ModelNonstat, m::ModelStat,
    F::Operator, optns::OptionsNonstat, statns::Stats,
    Temp::Float64, movetype::Int, current_misfit::Array{Float64, 1})
    # for purely nonstat move
    new_misfit = get_misfit(mns, optns, F)
    accepted = accept(current_misfit, new_misfit, statns, Temp, movetype)
    !accepted && undo_move!(movetype, mns, optns, m)
end

function mh_step!(m::ModelStat, mns::ModelNonstat,
    F::Operator, opt::OptionsStat, optns::OptionsNonstat,
    stat::Stats, Temp::Float64, movetype::Int, current_misfit::Array{Float64, 1})
    # for stat move updating nonstat
    new_misfit = get_misfit(mns, optns, F)
    accepted = accept(current_misfit, new_misfit, stat, Temp, movetype)
    !accepted && undo_move!(movetype, m, opt, mns, optns)
end

function do_mcmc_step(m::ModelStat, mns::ModelNonstat, mn::ModelNuisance,
    opt::OptionsStat, optns::OptionsNonstat, stat::Stats,
    current_misfit::Array{Float64, 1}, F::Operator,
    Temp::Float64, isample::Int, wp::Writepointers)
    # Stationary GP changes which update nonstationary GP + nuisance
    movetype, priorviolate = do_move!(m, opt, stat, mns, optns)
    if !priorviolate
        mh_step!(m, mns, mn, F, opt, optns, stat, Temp, movetype, current_misfit)
    end
    get_acceptance_stats!(isample, opt, stat)
    writemodel = false
    abs(Temp-1.0) < 1e-12 && (writemodel = true)
    write_history(isample, opt, m, current_misfit[1], stat, wp, Temp, writemodel)
    return current_misfit[1]
end

function do_mcmc_step(m::DArray{ModelStat}, mns::DArray{ModelNonstat},
    mn::DArray{ModelNuisance}, opt::DArray{OptionsStat},
    optns::DArray{OptionsNonstat}, stat::DArray{Stats},
    current_misfit::DArray{Array{Float64, 1}}, F::DArray{x}, Temp::Float64,
    isample::Int, wp::DArray{Writepointers}) where x<:Operator
    misfit = do_mcmc_step(localpart(m)[1], localpart(mns)[1], localpart(mn)[1],
        localpart(opt)[1], localpart(optns)[1], localpart(stat)[1],
        localpart(current_misfit)[1], localpart(F)[1],
        Temp, isample, localpart(wp)[1])
end

function do_mcmc_step(mns::ModelNonstat, m::ModelStat, mn::ModelNuisance,
    optns::OptionsNonstat, statns::Stats, current_misfit::Array{Float64, 1},
    F::Operator, Temp::Float64, isample::Int, wp::Writepointers)
    # purely nonstationary GP moves + nuisance
    movetype, priorviolate = do_move!(mns, m, optns, statns)
    if !priorviolate
        mh_step!(mns, m, mn, F, optns, statns, Temp, movetype, current_misfit)
    end
    get_acceptance_stats!(isample, optns, statns)
    writemodel = false
    abs(Temp-1.0) < 1e-12 && (writemodel = true)
    write_history(isample, optns, mns, current_misfit[1], statns, wp, Temp, writemodel)
    return current_misfit[1]
end

function do_mcmc_step(mns::DArray{ModelNonstat}, m::DArray{ModelStat},
    mn::DArray{ModelNuisance}, optns::DArray{OptionsNonstat}, statns::DArray{Stats},
    current_misfit::DArray{Array{Float64, 1}}, F::DArray{x}, Temp::Float64,
    isample::Int, wpns::DArray{Writepointers}) where x<:Operator
    misfit = do_mcmc_step(localpart(mns)[1], localpart(m)[1], localpart(mn)[1],
        localpart(optns)[1], localpart(statns)[1],localpart(current_misfit)[1],
        localpart(F)[1], Temp, isample, localpart(wpns)[1])
end

function do_mcmc_step(m::ModelStat, mn::ModelNuisance,
    opt::OptionsStat, stat::Stats, current_misfit::Array{Float64, 1},
    F::Operator, Temp::Float64, isample::Int, wp::Writepointers)
    # purely stationary GP moves + nuisance
    movetype, priorviolate = do_move!(m, opt, stat)
    if !priorviolate
        mh_step!(m, mn, F, opt, stat, Temp, movetype, current_misfit)
    end
    get_acceptance_stats!(isample, opt, stat)
    writemodel = false
    abs(Temp-1.0) < 1e-12 && (writemodel = true)
    write_history(isample, opt, m, current_misfit[1], stat, wp, Temp, writemodel)
    return current_misfit[1]
end

function do_mcmc_step(m::DArray{ModelStat}, mn::DArray{ModelNuisance},
    opt::DArray{OptionsStat}, stat::DArray{Stats},
    current_misfit::DArray{Array{Float64, 1}}, F::DArray{x}, Temp::Float64,
    isample::Int, wp::DArray{Writepointers}) where x<:Operator
    misfit = do_mcmc_step(localpart(m)[1], localpart(mn)[1],
        localpart(opt)[1], localpart(stat)[1],localpart(current_misfit)[1],
        localpart(F)[1], Temp, isample, localpart(wp)[1])
end

function do_mcmc_step(mn::ModelNuisance, m::Model,
    optn::OptionsNuisance, statn::Stats, current_misfit::Array{Float64,1},
    F::Operator, Temp::Float64, isample::Int, wpn::Writepointers_nuisance)
    # this only gets called for moves on the nuisance chain
    # can be either nonstat or stat model
    movetype = do_move!(mn, optn, statn)
    mh_step!(mn, m, F, optn, statn, Temp, movetype, current_misfit)
    get_acceptance_stats!(isample, optn, statn)
    writemodel = false
    abs(Temp - 1.0) < 1e-12 && (writemodel = true)
    write_history(isample, optn, mn, current_misfit[1], statn, wpn, Temp, writemodel)
    return current_misfit[1]
end

function do_mcmc_step(mn::DArray{ModelNuisance}, m::DArray{S},
    optn::DArray{OptionsNuisance}, statn::DArray{Stats},
    current_misfit::DArray{Array{Float64, 1}},
    F::DArray{x}, Temp::Float64, isample::Int,
    wpn::DArray{Writepointers_nuisance}) where {x<:Operator, S<:Model}
    misfit = do_mcmc_step(localpart(mn)[1], localpart(m)[1],
            localpart(optn)[1], localpart(statn)[1], localpart(current_misfit)[1],
            localpart(F)[1], localpart(Temp)[1], localpart(isample)[1], localpart(wpn)[1])
end

function do_mcmc_step(m::ModelStat, mns::ModelNonstat,
    opt::OptionsStat, optns::OptionsNonstat, stat::Stats,
    current_misfit::Array{Float64, 1}, F::Operator,
    Temp::Float64, isample::Int, wp::Writepointers)
    # Stationary GP changes which update nonstationary GP
    movetype, priorviolate = do_move!(m, opt, stat, mns, optns)
    if !priorviolate
        mh_step!(m, mns, F, opt, optns, stat, Temp, movetype, current_misfit)
    end
    get_acceptance_stats!(isample, opt, stat)
    writemodel = false
    abs(Temp-1.0) < 1e-12 && (writemodel = true)
    write_history(isample, opt, m, current_misfit[1], stat, wp, Temp, writemodel)
    return current_misfit[1]
end

function do_mcmc_step(m::DArray{ModelStat}, mns::DArray{ModelNonstat}, opt::DArray{OptionsStat},
    optns::DArray{OptionsNonstat}, stat::DArray{Stats},
    current_misfit::DArray{Array{Float64, 1}}, F::DArray{x}, Temp::Float64,
    isample::Int, wp::DArray{Writepointers}) where x<:Operator
    misfit = do_mcmc_step(localpart(m)[1], localpart(mns)[1],
        localpart(opt)[1], localpart(optns)[1], localpart(stat)[1],
        localpart(current_misfit)[1], localpart(F)[1],
        Temp, isample, localpart(wp)[1])
end

function do_mcmc_step(mns::ModelNonstat, m::ModelStat,
    optns::OptionsNonstat, statns::Stats, current_misfit::Array{Float64, 1},
    F::Operator, Temp::Float64, isample::Int, wpns::Writepointers)
    # purely nonstationary GP moves
    movetype, priorviolate = do_move!(mns, m, optns, statns)
    if !priorviolate
        mh_step!(mns, m, F, optns, statns, Temp, movetype, current_misfit)
    end
    get_acceptance_stats!(isample, optns, statns)
    writemodel = false
    abs(Temp-1.0) < 1e-12 && (writemodel = true)
    write_history(isample, optns, mns, current_misfit[1], statns, wpns, Temp, writemodel)
    return current_misfit[1]
end

function do_mcmc_step(mns::DArray{ModelNonstat}, m::DArray{ModelStat},
    optns::DArray{OptionsNonstat}, statns::DArray{Stats},
    current_misfit::DArray{Array{Float64, 1}}, F::DArray{x}, Temp::Float64,
    isample::Int, wpns::DArray{Writepointers}) where x<:Operator
    misfit = do_mcmc_step(localpart(mns)[1], localpart(m)[1],
        localpart(optns)[1], localpart(statns)[1],localpart(current_misfit)[1],
        localpart(F)[1], Temp, isample, localpart(wpns)[1])
end

function do_mcmc_step(m::ModelStat, opt::OptionsStat, stat::Stats, current_misfit::Array{Float64, 1},
    F::Operator, Temp::Float64, isample::Int, wp::Writepointers)
    # purely stationary GP moves
    movetype, priorviolate = do_move!(m, opt, stat)
    if !priorviolate
        mh_step!(m, F, opt, stat, Temp, movetype, current_misfit)
    end
    get_acceptance_stats!(isample, opt, stat)
    writemodel = false
    abs(Temp-1.0) < 1e-12 && (writemodel = true)
    write_history(isample, opt, m, current_misfit[1], stat, wp, Temp, writemodel)
    return current_misfit[1]
end

function do_mcmc_step(m::DArray{ModelStat}, opt::DArray{OptionsStat}, stat::DArray{Stats},
    current_misfit::DArray{Array{Float64, 1}}, F::DArray{x}, Temp::Float64,
    isample::Int, wp::DArray{Writepointers}) where x<:Operator
    misfit = do_mcmc_step(localpart(m)[1], localpart(opt)[1],
        localpart(stat)[1],localpart(current_misfit)[1],
        localpart(F)[1], Temp, isample, localpart(wp)[1])
end

function close_history(wp::DArray)
    for (idx, pid) in enumerate(procs(wp))
        @sync @spawnat pid close_history(wp[idx])
    end
end

function getlastiter(opt_in::Options, chains, idx)
    iterlast = 0
    if opt_in.history_mode=="a"
        # only the value of iterlast at last chain is used
        # as we want to avoid history reads
        if idx == length(chains)
            iterlast = history(opt_in, stat=:iter)[end]
        end
        chains[idx].T = history(opt_in, stat=:T)[end]
    end
    iterlast
end

function close_temperature_file(fp::IOStream)
    close(fp)
end

function makewritefilenames(opt_in)
    costs_filename = "misfits_"*opt_in.fdataname
    fstar_filename = "models_"*opt_in.fdataname
    x_ftrain_filename = "points_"*opt_in.fdataname
    nu_filename = "values_nuisance_"*opt_in.fdataname
    costs_filename, fstar_filename, x_ftrain_filename, nu_filename
end

function makestatfilenames(opt_in::OptionsStat, idx)
    costs_filename, fstar_filename, x_ftrain_filename, = makewritefilenames(opt_in)
    opt_in.costs_filename      = costs_filename*"s_$idx.bin"
    opt_in.fstar_filename      = fstar_filename*"s_$idx.bin"
    opt_in.x_ftrain_filename   = x_ftrain_filename*"s_$idx.bin"
end

function makenonstatfilenames(optns_in::OptionsNonstat, idx)
    costs_filename, fstar_filename, x_ftrain_filename, = makewritefilenames(optns_in)
    optns_in.costs_filename    = costs_filename*"ns_$idx.bin"
    optns_in.fstar_filename    = fstar_filename*"ns_$idx.bin"
    optns_in.x_ftrain_filename = x_ftrain_filename*"ns_$idx.bin"
end

function makenuisancefilenames(optn_in::OptionsNuisance, idx)
    costs_filename, fstar_filename, x_ftrain_filename, nu_filename = makewritefilenames(optn_in)
    optn_in.costs_filename     = costs_filename*"nuisance_$idx.bin"
    optn_in.vals_filename      = nu_filename*"$idx.bin"
end

function init_chain_darrays(opt_in::OptionsStat,
                            optns_in::OptionsNonstat,
                            optn_in::OptionsNuisance,
                            F_in::Operator, chains::Array{Chain, 1})
    # for nonstat, stat, and nuisances all together

    m_, mns_, mn_, opt_, optns_, optn_, F_in_, stat_, statns_, statn_,
    current_misfit_, wp_, wpns_, wpn_  = map(x -> Array{Future, 1}(undef, length(chains)), 1:14)

    opt_in.history_mode == "a" && setrestartflag.([opt_in, optns_in, optn_in])

    iterlast = 0
    @sync for(idx, chain) in enumerate(chains)

        makestatfilenames(opt_in, idx)
        makenonstatfilenames(optns_in, idx)
        makenuisancefilenames(optn_in, idx)

        opt_[idx]            = @spawnat chain.pid [opt_in]
        optns_[idx]          = @spawnat chain.pid [optns_in]
        optn_[idx]           = @spawnat chain.pid [optn_in]

        m_[idx]              = @spawnat chain.pid [init(opt_in)]
        mns_[idx]            = @spawnat chain.pid [init(optns_in, fetch(m_[idx])[1])]
        mn_[idx]             = @spawnat chain.pid [init(optn_in)]

        @sync wp_[idx]       = @spawnat chain.pid [open_history(opt_in)]
        @sync wpns_[idx]     = @spawnat chain.pid [open_history(optns_in)]
        @sync wpn_[idx]      = @spawnat chain.pid [open_history(optn_in)]

        stat_[idx]           = @spawnat chain.pid [Stats()]
        statns_[idx]         = @spawnat chain.pid [Stats()]
        statn_[idx]          = @spawnat chain.pid [Stats(nmoves=optn_in.nnu)]

        F_in_[idx]           = @spawnat chain.pid [F_in]

        current_misfit_[idx] = @spawnat chain.pid [[ get_misfit(fetch(mns_[idx])[1],
                                fetch(mn_[idx])[1],
                                fetch(optns_[idx])[1],
                                fetch(F_in_[idx])[1]) ]]

        iterlast = getlastiter(opt_in, chains, idx)
    end

    m, mns, mn, opt, optns, optn, stat, statns, statn, F,
    current_misfit, wp, wpns, wpn = map(x -> DArray(x), (m_, mns_, mn_, opt_, optns_, optn_,
                                    stat_, statns_, statn_, F_in_, current_misfit_,
                                    wp_, wpns_, wpn_))
    @info "initialisation complete"
    return m, mns, mn, opt, optns, optn, stat, statns, statn, F, current_misfit,
            wp, wpns, wpn, iterlast
end

function domcmciters(iterlast, nsamples, chains, mns::DArray{ModelNonstat}, m::DArray{ModelStat},
            mn::DArray{ModelNuisance}, optns::DArray{OptionsNonstat}, opt::DArray{OptionsStat},
            optn::DArray{OptionsNuisance}, statns, stat, statn, current_misfit, F, wpns, wp, wpn)
    # for nonstat, stat, and nuisances all together

    t2 = time()
    for isample = iterlast+1:iterlast+nsamples
        # we do need each remotecall to finish before
        # moving on to the next kind of move
        swap_temps(chains)
        @sync for chain in chains
            # purely nonstationary GP moves + nuisance
            @async chain.misfit = remotecall_fetch(do_mcmc_step, chain.pid,
                                            mns, m, mn, optns, statns,
                                            current_misfit, F,
                                            chain.T, isample, wpns)
        end
        @sync for chain in chains
            # purely nuisance move
            @async chain.misfit = remotecall_fetch(do_mcmc_step, chain.pid,
                                            mn, m, optn, statn,
                                            current_misfit, F,
                                            chain.T, isample, wpn)
        end
        @sync for chain in chains
            # Stationary GP changes which update nonstationary GP + nuisance
            @async chain.misfit = remotecall_fetch(do_mcmc_step, chain.pid,
                                            m, mns, mn, opt, optns, stat,
                                            current_misfit, F,
                                            chain.T, isample, wp)
        end
        t2 = disptime(isample, t2, iterlast, nsamples)
    end
end

function main(opt_in     ::OptionsStat,
            optns_in     ::OptionsNonstat,
            optn_in      ::OptionsNuisance,
            F_in         ::Operator,
            chainprocs   ::Array{Int, 1};
            nsamples     = 4001,
            nchainsatone = 1,
            Tmax         = 2.5)
    # for nonstat, stat, and nuisances all together

    chains = Chain(chainprocs, nchainsatone=nchainsatone, Tmax=Tmax)

    m, mns, mn, opt, optns, optn, stat,
    statns, statn, F, current_misfit, wp, wpns, wpn,
    iterlast = init_chain_darrays(opt_in, optns_in, optn_in, F_in, chains)

    domcmciters(iterlast, nsamples, chains, mns, m, mn, optns, opt,
        optn, statns, stat, statn, current_misfit, F, wpns, wp, wpn)

    close_history.([wp, wpns, wpn])
    nothing
end

function main(opt_in ::OptionsStat,
    optns_in     ::OptionsNonstat,
    optn         ::OptionsNuisance,
    F_in         ::Operator;
    nchains      = 4,
    nsamples     = 4001,
    nchainsatone = 1,
    Tmax         = 2.5)
    # unnecessary but for backwards compat
    # for nonstat and stat together
    pids = workers()
    @assert length(pids) == nchains
    main(opt_in, optns_in, optn, F_in, workers(), nsamples=nsamples, Tmax=Tmax, nchainsatone=nchainsatone)
end

function init_chain_darrays(opt_in::OptionsStat,
                            optns_in::OptionsNonstat,
                            F_in::Operator, chains::Array{Chain, 1})
    # for nonstat and stat together

    m_, mns_, opt_, optns_, F_in_, stat_, statns_,
    current_misfit_, wp_, wpns_,  = map(x -> Array{Future, 1}(undef, length(chains)), 1:10)

    opt_in.history_mode == "a" && setrestartflag.([opt_in, optns_in])

    iterlast = 0
    @sync for(idx, chain) in enumerate(chains)

        makestatfilenames(opt_in, idx)
        makenonstatfilenames(optns_in, idx)

        opt_[idx]            = @spawnat chain.pid [opt_in]
        optns_[idx]          = @spawnat chain.pid [optns_in]

        m_[idx]              = @spawnat chain.pid [init(opt_in)]
        mns_[idx]            = @spawnat chain.pid [init(optns_in, fetch(m_[idx])[1])]

        @sync wp_[idx]       = @spawnat chain.pid [open_history(opt_in)]
        @sync wpns_[idx]     = @spawnat chain.pid [open_history(optns_in)]

        stat_[idx]           = @spawnat chain.pid [Stats()]
        statns_[idx]         = @spawnat chain.pid [Stats()]

        F_in_[idx]           = @spawnat chain.pid [F_in]

        current_misfit_[idx] = @spawnat chain.pid [[ get_misfit(fetch(mns_[idx])[1],
                                fetch(optns_[idx])[1],
                                fetch(F_in_[idx])[1]) ]]

        iterlast = getlastiter(opt_in::Options, chains, idx)
    end

    m, mns, opt, optns, stat, statns, F,
    current_misfit, wp, wpns = map(x -> DArray(x), (m_, mns_, opt_, optns_,
                                    stat_, statns_, F_in_, current_misfit_,
                                    wp_, wpns_))
    @info "initialisation complete"
    return m, mns, opt, optns, stat, statns, F, current_misfit,
            wp, wpns, iterlast
end

function domcmciters(iterlast, nsamples, chains, mns::DArray{ModelNonstat}, m::DArray{ModelStat},
            optns::DArray{OptionsNonstat}, opt::DArray{OptionsStat},
            statns, stat, current_misfit, F, wpns, wp)
    # for nonstat and stat together

    t2 = time()
    for isample = iterlast+1:iterlast+nsamples
        # we do need each remotecall to finish before
        # moving on to the next kind of move
        swap_temps(chains)
        @sync for chain in chains
            # purely nonstationary GP moves
            @async chain.misfit = remotecall_fetch(do_mcmc_step, chain.pid,
                                            mns, m, optns, statns,
                                            current_misfit, F,
                                            chain.T, isample, wpns)
        end
        @sync for chain in chains
            # Stationary GP changes which update nonstationary GP
            @async chain.misfit = remotecall_fetch(do_mcmc_step, chain.pid,
                                            m, mns, opt, optns, stat,
                                            current_misfit, F,
                                            chain.T, isample, wp)
        end
        t2 = disptime(isample, t2, iterlast, nsamples)
    end
end

function main(opt_in     ::OptionsStat,
            optns_in     ::OptionsNonstat,
            F_in         ::Operator,
            chainprocs   ::Array{Int, 1};
            nsamples     = 4001,
            nchainsatone = 1,
            Tmax         = 2.5)
    # for nonstat and stat together

    chains = Chain(chainprocs, nchainsatone=nchainsatone, Tmax=Tmax)

    m, mns, opt, optns, stat,
    statns, F, current_misfit, wp, wpns,
    iterlast = init_chain_darrays(opt_in, optns_in, F_in, chains)

    domcmciters(iterlast, nsamples, chains, mns, m, optns, opt,
        statns, stat, current_misfit, F, wpns, wp)

    close_history.([wp, wpns])
    nothing
end

function main(opt_in ::OptionsStat,
    optns_in     ::OptionsNonstat,
    F_in         ::Operator;
    nchains      = 4,
    nsamples     = 4001,
    nchainsatone = 1,
    Tmax         = 2.5)
    # unnecessary but for backwards compat
    # for nonstat and stat together
    pids = workers()
    @assert length(pids) == nchains
    main(opt_in, optns_in, F_in, workers(), nsamples=nsamples, Tmax=Tmax, nchainsatone=nchainsatone)
end

function init_chain_darrays(opt_in::OptionsStat, optn_in::OptionsNuisance,
                            F_in::Operator, chains::Array{Chain, 1})
    # purely stationary GP moves + nuisance

    m_, mn_, opt_, optn_, F_in_, stat_, statn_,
    current_misfit_, wp_, wpn_  = map(x -> Array{Future, 1}(undef, length(chains)), 1:10)

    opt_in.history_mode == "a" && setrestartflag.([opt_in, optn_in])

    iterlast = 0
    @sync for(idx, chain) in enumerate(chains)

        makestatfilenames(opt_in, idx)
        makenuisancefilenames(optn_in, idx)

        opt_[idx]            = @spawnat chain.pid [opt_in]
        optn_[idx]           = @spawnat chain.pid [optn_in]

        m_[idx]              = @spawnat chain.pid [init(opt_in)]
        mn_[idx]             = @spawnat chain.pid [init(optn_in)]

        @sync wp_[idx]       = @spawnat chain.pid [open_history(opt_in)]
        @sync wpn_[idx]      = @spawnat chain.pid [open_history(optn_in)]

        stat_[idx]           = @spawnat chain.pid [Stats()]
        statn_[idx]          = @spawnat chain.pid [Stats(nmoves=optn_in.nnu)]

        F_in_[idx]           = @spawnat chain.pid [F_in]

        current_misfit_[idx] = @spawnat chain.pid [[ get_misfit(fetch(m_[idx])[1],
                                fetch(mn_[idx])[1],
                                fetch(opt_[idx])[1],
                                fetch(F_in_[idx])[1]) ]]

        iterlast = getlastiter(opt_in, chains, idx)
    end

    m, mn, opt, optn, stat, statn, F,
    current_misfit, wp, wpn = map(x -> DArray(x), (m_, mn_, opt_, optn_,
                                    stat_, statn_, F_in_, current_misfit_,
                                    wp_, wpn_))
    @info "initialisation complete"
    return m, mn, opt, optn, stat, statn, F, current_misfit,
        wp, wpn, iterlast
end

function domcmciters(iterlast, nsamples, chains, m::DArray{ModelStat}, mn::DArray{ModelNuisance},
            opt::DArray{OptionsStat}, optn::DArray{OptionsNuisance}, stat, statn,
            current_misfit, F, wp, wpn)
    # purely stationary GP moves + nuisance

    t2 = time()
    for isample = iterlast+1:iterlast+nsamples
        # we do need each remotecall to finish before
        # moving on to the next kind of move
        swap_temps(chains)
        @sync for chain in chains
            # purely nuisance move
            @async chain.misfit = remotecall_fetch(do_mcmc_step, chain.pid,
                                            mn, m, optn, statn,
                                            current_misfit, F,
                                            chain.T, isample, wpn)
        end
        @sync for chain in chains
            # purely stationary GP moves + nuisance
            @async chain.misfit = remotecall_fetch(do_mcmc_step, chain.pid,
                                            m, mn, opt, stat,
                                            current_misfit, F,
                                            chain.T, isample, wp)
        end
        t2 = disptime(isample, t2, iterlast, nsamples)
    end
end

function main(opt_in       ::OptionsStat,
        optn_in      ::OptionsNuisance,
        F_in         ::Operator,
        chainprocs   ::Array{Int, 1};
        nsamples     = 4001,
        nchainsatone = 1,
        Tmax         = 2.5)
    # purely stationary GP moves + nuisance

    chains = Chain(chainprocs, nchainsatone=nchainsatone, Tmax=Tmax)

    m, mn, opt, optn, stat,
    statn, F, current_misfit, wp, wpn,
    iterlast = init_chain_darrays(opt_in, optn_in, F_in, chains)

    domcmciters(iterlast, nsamples, chains, m, mn, opt,
        optn, stat, statn, current_misfit, F, wp, wpn)

    close_history.([wp, wpn])
    nothing
end

function main(opt_in ::OptionsStat,
    optn_in     ::OptionsNuisance,
    F_in         ::Operator;
    nchains      = 4,
    nsamples     = 4001,
    nchainsatone = 1,
    Tmax         = 2.5)
    # unnecessary but for backwards compat
    # purely stationary GP moves + nuisance
    pids = workers()
    @assert length(pids) == nchains
    main(opt_in, optn_in, F_in, workers(), nsamples=nsamples, Tmax=Tmax, nchainsatone=nchainsatone)
end

function init_chain_darrays(opt_in::OptionsStat, F_in::Operator,
    chains::Array{Chain, 1})
    # purely stationary GP moves

    m_, opt_, F_in_, stat_, current_misfit_,
    wp_  = map(x -> Array{Future, 1}(undef, length(chains)), 1:6)

    opt_in.history_mode == "a" && setrestartflag(opt_in)

    iterlast = 0
    @sync for(idx, chain) in enumerate(chains)

        makestatfilenames(opt_in, idx)

        opt_[idx]            = @spawnat chain.pid [opt_in]
        m_[idx]              = @spawnat chain.pid [init(opt_in)]
        @sync wp_[idx]       = @spawnat chain.pid [open_history(opt_in)]
        stat_[idx]           = @spawnat chain.pid [Stats()]
        F_in_[idx]           = @spawnat chain.pid [F_in]
        current_misfit_[idx] = @spawnat chain.pid [[ get_misfit(fetch(m_[idx])[1],
                                fetch(opt_[idx])[1],
                                fetch(F_in_[idx])[1]) ]]
        iterlast = getlastiter(opt_in, chains, idx)
    end

    m, opt, stat, F,
    current_misfit, wp = map(x -> DArray(x), (m_, opt_,
                                    stat_, F_in_, current_misfit_,
                                    wp_))
    @info "initialisation complete"
    return m, opt, stat, F, current_misfit,
        wp, iterlast
end

function domcmciters(iterlast, nsamples, chains, m::DArray{ModelStat},
            opt::DArray{OptionsStat}, stat,
            current_misfit, F, wp)
    # purely stationary GP moves

    t2 = time()
    for isample = iterlast+1:iterlast+nsamples
        swap_temps(chains)
        @sync for chain in chains
            # purely stationary GP moves
            @async chain.misfit = remotecall_fetch(do_mcmc_step, chain.pid,
                                            m, opt, stat,
                                            current_misfit, F,
                                            chain.T, isample, wp)
        end
        t2 = disptime(isample, t2, iterlast, nsamples)
    end
end

function main(opt_in ::OptionsStat,
        F_in         ::Operator,
        chainprocs   ::Array{Int, 1};
        nsamples     = 4001,
        nchainsatone = 1,
        Tmax         = 2.5)
    # purely stationary GP moves

    chains = Chain(chainprocs, nchainsatone=nchainsatone, Tmax=Tmax)

    m, opt, stat,
    F, current_misfit, wp,
    iterlast = init_chain_darrays(opt_in, F_in, chains)

    domcmciters(iterlast, nsamples, chains, m, opt,
        stat, current_misfit, F, wp)

    close_history(wp)
    nothing
end

function main(opt_in ::OptionsStat,
    F_in         ::Operator;
    nchains      = 4,
    nsamples     = 4001,
    nchainsatone = 1,
    Tmax         = 2.5)
    # unnecessary but for backwards compat
    # purely stationary GP moves
    pids = workers()
    @assert length(pids) == nchains
    main(opt_in, F_in, workers(), nsamples=nsamples, Tmax=Tmax, nchainsatone=nchainsatone)
end

function swap_temps(chains::Array{Chain, 1})
    for ichain in length(chains):-1:2
        jchain = rand(1:ichain)
        if ichain != jchain
            logalpha = (chains[ichain].misfit - chains[jchain].misfit) *
                            (1.0/chains[ichain].T - 1.0/chains[jchain].T)
            if log(rand()) < logalpha
                chains[ichain].T, chains[jchain].T = chains[jchain].T, chains[ichain].T
            end
        end
    end
end

function disptime(isample, t2, iterlast, nsamples)
    if mod(isample-1, 1000) == 0
        dt = time() - t2 #seconds
        t2 = time()
        @info("**$dt**sec** $isample out of $(iterlast+nsamples)")
    end
    t2
end

max_stars_repo_name: GeoscienceAustralia/HiQGA

max_issues_repo_name: GeoscienceAustralia/HiQGA


text:
include("../examples/espprc_example.jl")

@info("Testing for OD=($(pg.origin), $(pg.destination))")

print("Pulse     : "); @time pulse = solveESPPRC(pg, method="pulse")
print("Mono      : "); @time mono0= solveESPPRC(pg, method="monodirectional")
print("Mono DSSR : "); @time mono1 = solveESPPRC(pg, method="monodirectional", DSSR=true)
# print("Bi        : "); @time bidi0= solveESPPRC(pg, method="bidirectional")
# print("Bi   DSSR : "); @time bidi1 = solveESPPRC(pg, method="bidirectional", DSSR=true)

@show pulse.cost, pulse.load, pulse.time
@show mono0.cost, mono0.load, mono0.time
@show mono1.cost, mono1.load, mono1.time
# @show bidi0.cost, bidi0.load, bidi0.time
# @show bidi1.cost, bidi1.load, bidi1.time
@show pulse.path
@show mono0.path
@show mono1.path
# @show bidi0.path
# @show bidi1.path

@testset "ESPPRC Example Test" begin
    @test isapprox(pulse.cost, mono0.cost, atol=1e-7)
    @test isapprox(mono0.cost, mono1.cost, atol=1e-7)
    # @test isapprox(pulse.cost, bidi0.cost, atol=1e-7)
    # @test isapprox(bidi0.cost, bidi1.cost, atol=1e-7)
end

# negative cost route return test
max_neg = 10
@info("ESPPRC Example testing with max_neg_routes=$(max_neg)...")

@time pulse, neg_pulse = solveESPPRC_vrp(pg, method="pulse", max_neg_routes=max_neg)
@time mono0, neg_mono0 = solveESPPRC_vrp(pg, method="monodirectional", max_neg_routes=max_neg)
@time mono1, neg_mono1 = solveESPPRC_vrp(pg, method="monodirectional", max_neg_routes=max_neg, DSSR=true)

@show typeof(pulse), typeof(neg_pulse), length(neg_pulse)
@show pulse
@show mono0
@show mono1

@testset "ESPPRC Example Test with max_neg_routes" begin
    @test length(neg_pulse) <= max_neg
    @test length(neg_mono0) <= max_neg
    @test length(neg_mono1) <= max_neg
    for l in neg_pulse
        @test l.cost < 0.0
        @test l.path[1] == pg.origin
        @test l.path[end] == pg.destination
    end
    for l in neg_mono0
        @test l.cost < 0.0
        @test l.path[1] == pg.origin
        @test l.path[end] == pg.destination
    end
    for l in neg_mono1
        @test l.cost < 0.0
        @test l.path[1] == pg.origin
        @test l.path[end] == pg.destination
    end
end

ext: jl


text:
using SparseArrays

"""
    ntrace(X::Space, γ::Mesh)

Compute the normal trace of basis X on mesh γ. γ is assumed to be part of the
boundary of geometry(X).
"""
function ntrace(X::DirectProductSpace{T}, γ) where T
    x = Space{T}[ntrace(s, γ) for s in X.factors]
    return DirectProductSpace(x)
end

function strace(X::DirectProductSpace{T}, γ) where T
    x = Space{T}[strace(s,γ) for s in X.factors]
    return DirectProductSpace(x)
end

function ntrace(X::Space, γ)

    # on_target = overlap_gpredicate(γ)
    # ad = assemblydata(X)
    x = refspace(X)
    E, ad = assemblydata(X)
    igeo = geometry(X)
    @assert dimension(γ) == dimension(igeo)-1
    # Γ = geo
    # Dγ = dimension(γ)
    # Σ = skeleton(Γ,Dγ)

    ogeo = boundary(igeo)
    on_target = overlap_gpredicate(γ)
    ogeo = submesh(ogeo) do face
        on_target(chart(ogeo,face))
    end

    D = copy(transpose(connectivity(ogeo, igeo, abs)))
    rows, vals = rowvals(D), nonzeros(D)

    T = scalartype(X)
    S = Shape{T}
    fns = [Vector{S}() for i in 1:numfunctions(X)]

    for (p,el) in enumerate(E)

        for (q,fc) in enumerate(faces(el))
            on_target(fc) || continue

            # print(Q)
            # @assert norm(Q,Inf) != 0

            r = 0
            for k in nzrange(D,p)
                vals[k] == q && (r = rows[k]; break)
            end
            @assert r != 0

            fc1 = chart(ogeo, cells(ogeo)[r])
            Q = ntrace(x, el, q, fc1)

            for i in 1:size(Q,1)
                for j in 1:size(Q,2)
                    for (m,a) in ad[p,j]
                        # j == q && println("bingo",j,q)
                        v = a*Q[i,j]
                        isapprox(v,0,atol=sqrt(eps(T))) && continue
                        push!(fns[m], Shape(r, i, v))
                    end
                end
            end

        end

    end

    ntrace(X, ogeo, fns)
end

strace(X::Space, γ::Mesh) = strace(X, γ, Val{dimension(γ)+1})

function strace(X::Space, γ, dim1::Type{Val{2}})

    x = refspace(X)
    on_target = overlap_gpredicate(γ)
    E, ad = assemblydata(X)

    geo = geometry(X)
    Γ = geo
    Σ = skeleton(Γ,1)

    # Allows for quick retrieval of the faces in Σ in a given cell of Γ.
    D = copy(transpose(connectivity(Σ, Γ, abs)))
    rows, vals = rowvals(D), nonzeros(D)

    T = scalartype(X)
    @assert numfunctions(X) != 0
    fns = [Shape{T}[] for i in 1:numfunctions(X)]

    for (p,el) in enumerate(E)

        for q in 1:dimension(Γ)+1

            # find the global index in Σ of the q-th face of the p-element
            e = 0
            for k in nzrange(D,p)
                if vals[k] == q
                    e = rows[k]
                    break
                end
            end
            @assert e != 0

            # make sure we use fc as oriented in Σ
            cell = Σ.faces[e]
            fc = chart(Σ, cell)
            on_target(fc) || continue

            Q = strace(x,el,q,fc)
            for i in 1:size(Q,1)
                for j in 1:size(Q,2)
                    for (m,a) in ad[p,j]
                        v = a*Q[i,j]
                        v == 0 && continue
                        push!(fns[m], Shape(e, i, v))
                    end
                end
            end

        end

    end

    strace(X, Σ, fns)
end

"""
currently not working!
"""
# function ttrace(X::Space, γ)
#
#     x = refspace(X)
#     on_target = overlap_gpredicate(γ)
#     # ad = assemblydata(X)
#     E, ad = assemblydata(X)
#
#     geo = geometry(X)
#     Γ = geo
#     Dγ = dimension(γ)
#     Σ = skeleton(Γ,Dγ)
#
#     D = copy(transpose(connectivity(Σ, Γ, abs)))
#     rows, vals = rowvals(D), nonzeros(D)
#
#     T = scalartype(X)
#     fns = [Shape{T}[] for i in 1:numfunctions(X)]
#
#     for (p,el) in enumerate(E)
#
#         for (q,fc) in enumerate(faces(el))
#
#             on_target(fc) || (println("skip"); continue)
#             # print("\n")
#             Q = ttrace(x,el,q,fc)
#             # print(Q,q)
#             # find the global index in Σ of the q-th face of the p-element
#             r = 0
#             for k in nzrange(D,p)
#                 vals[k] == q && (r = rows[k]; break)
#             end
#             @assert r != 0
#             for i in 1:size(Q,1)
#                 for j in 1:size(Q,2)
#                     for (m,a) in ad[p,j]
#                         # if j == q
#                         #     # print("bingo",a,Q[i,j])
#                         # end
#                         v = a*Q[i,j]
#                         v == 0 && continue
#                         push!(fns[m], Shape(r, i, v))
#                     end
#                 end
#             end
#
#         end
#
#     end
#
#     ttrace(X, Σ, fns)
# end

function ttrace(X::Space, γ)

    x = refspace(X)
    E, ad = assemblydata(X)
    igeo = geometry(X)
    @assert dimension(γ) == dimension(igeo)-1

    # ogeo = skeleton(igeo, dimension(γ))
    ogeo = boundary(igeo)
    on_target = overlap_gpredicate(γ)
    ogeo = submesh(ogeo) do face
        on_target(chart(ogeo, face))
    end

    # D = copy(transpose(connectivity(ogeo, igeo, abs)))
    D = connectivity(igeo, ogeo, abs)
    rows, vals = rowvals(D), nonzeros(D)

    T = scalartype(X)
    S = Shape{T}
    fns = [Vector{S}() for i in 1:numfunctions(X)]

    for (p,el) in enumerate(E)

        for (q,fc) in enumerate(faces(el))
            on_target(fc) || continue

            r = 0
            for k in nzrange(D,p)
                vals[k] == q && (r = rows[k]; break)
            end
            @assert r != 0

            fc1 = chart(ogeo, cells(ogeo)[r])
            @assert cartesian(center(fc)) ≈ cartesian(center(fc1))

            Q = ttrace(x, el, q, fc1)
            for i in 1:size(Q,1)
                for j in 1:size(Q,2)
                    for (m,a) in ad[p,j]
                        v = a*Q[i,j]
                        isapprox(v,0,atol=sqrt(eps(T))) && continue
                        push!(fns[m], Shape(r, i, v))
                    end
                end
            end

        end
    end

    ttrace(X, ogeo, fns)
end


text:
"""
Core function for computing splits in MCCs. Include root of MCC if it's not root of tree.
"""
function _splits_in_mcc(
	m::Array{<:AbstractString},
	t::Tree,
	leaves::Array{<:AbstractString},
	leafmap::Dict
)
    # Mask corresponding to leaves in the MCC
    mask = zeros(Bool, length(leaves))
    for l in m
        mask[leafmap[l]] = true
    end
    #
    r = lca(t.lleaves[x] for x in m)
    S = SplitList(r, leaves, mask, leafmap)
    TreeTools.clean!(S, clean_root=false) # if true, do not resolve based on root of MCC

    return S
end
"""
    splits_in_mcc(m::Array{<:AbstractString}, t...)

Compute splits in MCC `m` in tree `t`. Return a `SplitList` object, with `mask` corresponding to leaves in `m`.
"""
function splits_in_mcc(m::Array{<:AbstractString}, t::Tree)
    leaves = sort(collect(keys(t.lleaves)))
    leafmap = Dict(leaf=>i for (i,leaf) in enumerate(leaves))
    return _splits_in_mcc(m, t, leaves, leafmap)
end
splits_in_mcc(m::Array{<:AbstractString}, trees::Vararg{Tree}) = Tuple(splits_in_mcc(m,t) for t in trees)

"""
    splits_in_mccs(MCCs, t::Vararg{Tree})

Compute splits in `MCCs` in tree `t`. Return an array of `SplitList` objects, with `mask` corresponding to leaves in each MCC.
When called with multiple trees, return a tuple of arrays of `SplitList` objects, each of which corresponds to one of the trees (order conserved).
"""
function splits_in_mccs(MCCs, t::Tree)
    leaves = sort(collect(keys(t.lleaves)))
    leafmap = Dict(leaf=>i for (i,leaf) in enumerate(leaves))
    return [_splits_in_mcc(m, t, leaves, leafmap) for m in MCCs]
end
splits_in_mccs(MCCs, trees::Vararg{Tree}) = Tuple(splits_in_mccs(MCCs,t) for t in trees)

"""
    new_splits(MCCs, t1::Tree, t2::Tree)

Given MCCs and two trees, what are the new splits introduced in either of the trees.
Return an array `S` of `SplitList` objects, with `S[i][m]` corresponding to new splits in tree `i` from MCC `m`.

*Note*: Resulting splits are not always unique once constrained to leaves of MCCs!
"""
function new_splits(MCCs, t1::Tree, t2::Tree)
    nS = [
    	new_splits(t1, MCCs, t2), # Using each tree as a reference in turn
    	new_splits(t2, MCCs, t1),
    ]
    return nS
end

"""
    new_splits(tref::Tree, MCCs, t::Tree)

What are the splits introduced in `tref` from `t` using consistent clades `MCCs`?
Return a single `SplitList` object, with unique splits mapped onto the leaves of `tref`.
"""
function new_splits(tref::Tree, MCCs, t::Tree)
    # Splits in `tref`
    S_ref = SplitList(tref)
    # Splits corresponding to each mcc in tree `t`
    MCC_splits = TreeTools.map_splits_to_tree(splits_in_mccs(MCCs, t), tref)
    # Take the new ones only
    _new_splits!(MCC_splits, S_ref)
    # Map them onto leaves of `tref`
    return unique(TreeTools.map_splits_to_tree(MCC_splits, tref), usemask=false)
end

function _new_splits!(MCC_splits::SplitList, tree_splits::SplitList)
    idx = Int64[]
    for (i,s) in enumerate(MCC_splits)
        in(s, tree_splits, MCC_splits.mask) && push!(idx,i)
    end
    deleteat!(MCC_splits.splits, idx)
end

#=
To move to `artificialdata.jl`?
=#
"""
Indices of true splits in `S` w.r. to `Sref`.
"""
function true_splits(S::SplitList, Sref::SplitList, mask=S.mask)
    idx = Int64[]
    for (i,s) in enumerate(S)
        in(s, Sref, mask) && push!(idx, i)
    end
    return idx
end
true_splits(S, tref::Tree, mask=S.mask) = true_splits(S, SplitList(tref), mask)

max_stars_repo_name: PierreBarrat/RecombTools

max_issues_repo_name: PierreBarrat/TreeKnit

max_forks_repo_name: PierreBarrat/TreeKnit


text:
# okaty, this is the code from the sequential decision problem that stumped milton friedman, frmo here https://lectures.quantecon.org/jl/wald_friedman.html

using Distributions
using QuantEcon.compute_fixed_point, QuantEcon.DiscreteRV, QuantEcon.draw, QuantEcon.LinInterp
using Plots
pyplots() #can't remember if thi scauses issues outside of jupyter notebook
using LateXStrings

#okay, our los functions

# for given probability return expected los of choosing model 0

function expect_loss_choose_0(p::Real, L0::Real)
	return (1-p)*L0
end

#for given probability return expected lsos of choosing model 1

function expect_loss_choose_1(p::Real, L1::Real)
	return p*L1
end

#we need to evaluate  the expectation of our bellman equation J. so we need current probability that model 0 is correct, distributions are f0, f1, and a function that can evaluate bellman equatoin

function EJ(p::Real, f0::AbstractVector, f1::AbstractVector, J::LinInterp)
	#get current distribution we believe - i.e p*fo + 1-p*f1

	curr_dist = p*f0 + (1-p)*f1

	#get tomorrows expected distribution through bayes law
	tp1_dist = clamp.((p*f0) ./(p*f0 (1-p)*f1), 0,1)

	#evaluate the expectation
	EJ = dot(curr_dist, J.tp1_dist))
	return EJ
end

# evaluate valiye function for given continuation value

function expect_loss_cont(p::Real, c::Real, f0::AbstractVector, f1::AbstractVector, j::LinInterp)
	return c+EJ(p, f0, f1,J)
end

function bellman_operator(pgrid::AbstractVector, c::Real, f0::AbstractVector, f1::AbstractVector, L0::Real, L1::Real, J::AbstractVector)
	m = length(pgrid)
	@assert m == length(J)

	J_out = zeros(m)
	j_interp = LinInterp(pgrid, J)

	for (p_ind, p) in enumerate(pgrid)
		#payoff of chooseing model 0
		p_c_0 = expect_loss_choose_0(p, L0)
		p_c_1 expect_loss_choose_1(p,L1)
		p_con = expect_loss_cont(p,c,f0,f1,J_interp)

		J_out[p_ind] = min(p_c_0, p_c_1, p_con)
	end
	return J_out
end

# this gets us through the bellman equatino, no wwe just create two istributions over 50 values for k usingdiscrete beta distribution

p_m1 = linspace(0,1,50)
f0 = clamp.(pdf(Beta(1,1),p_m1),1e-8, Inf)
f0 = f0/sum(f0)
f1 = clam.(pdf(Beta(9,9), p_m1), 1e-8, Inf)
f1 = f1/sum(f1)

# to solve
pg = linspace(0,1,251)
J1 = compute_fixed_point(x->bellman_operator(pg, 0.5, f0,f1,5.0,5.0,x), zeros(length(pg)), err_tol=1e-6, print_skip = 5)

# okay, sorted, that's kind of crazy to be honest, but it might or might not work. dagnabbit julia. good ol' solving stuff marting freedman couldn't. there's a better implementation also with more comments and so forth, and we'll do that in a bit. I also need to do my actual phd work. But I kind of understand how this works now, and how powerful bellmans and dynamic programming isfor simple things like this!

max_stars_repo_name: Bmillidgework/Misc-Maths

max_issues_repo_name: Bmillidgework/Misc-Maths

max_forks_repo_name: Bmillidgework/Misc-Maths


text:
function bode_immerse(sys, w, ymin, ymax)
  mag, phase, _ = bode(sys, w)
  I = find(y-> ymin .<= y .<= ymax, mag) # Filter out from range
  mag = mag[I]
  ws = w[I]
  xmajor, xminor, xmajorText, xminorText = get_log_ticks(w)
  ymajor, yminor, ymajorText, yminorText = get_log_ticks([ymin;mag;ymax])
  Immerse.plot(
        major_lines(log10(ws), log10(mag), log10(xmajor), log10(ymajor)),
        layer(x=log10(ws), y=log10(mag), Geom.line),
        Stat.xticks(ticks = log10([xmajor;xminor])),
        Stat.yticks(ticks = log10([ymajor;yminor])),
        Scale.x_continuous(labels= x -> labelf(x), maxvalue=log10(maximum(w))),
        Scale.y_continuous(labels= x -> labelf(x), minvalue=log10(ymin), maxvalue=log10(ymax)),
        Guide.xlabel("Frequency (rad/s)"),
        Guide.ylabel("Magnitude (abs)", orientation=:vertical),
        default_theme)
end
function nyquist_immerse(sys, w, xmin, xmax, ymin, ymax)
  re, im, _ = nyquist(sys, w)
  I = (ymin .<= im .<= ymax) & (xmin .<= re .<= xmax)
  re = re[I]
  im = im[I]
  Immerse.plot(
        major_lines(re, im, [0.0], [0.0]),
        layer(x=[-1], y=[0], Geom.point, Theme(default_color=Colors.RGB(.8,.0,.0), default_point_size=5pt)),
        layer(x=re, y = im, Geom.line(preserve_order=true)),
        Scale.x_continuous(minvalue = xmin, maxvalue=xmax),
        Scale.y_continuous(minvalue = ymin, maxvalue=ymax),
        Guide.xlabel("Real"),
        Guide.ylabel("Imaginary", orientation=:vertical),
        default_theme)
end

function pzmap_immerse(sys, xmin, xmax, ymin, ymax)
  z, p, k = zpkdata(sys)
  z = z[1]; p = p[1];
  I(x) = (ymin .<= imag(x) .<= ymax) & (xmin .<= real(x) .<= xmax)
  zx = real(z[I(z)])
  zy = imag(z[I(z)])
  px = real(p[I(p)])
  py = imag(p[I(p)])
  shapes = xcross(px, py, fill(3mm,length(px)))
  theme = copy(default_theme)
  theme.default_color = Colors.RGB(0,0,1)
  theme.default_point_size = 6pt
  Immerse.plot(
        layer(x=zx, y=zy, Geom.point),
        x = px,
        y = py,
        Scale.x_continuous(minvalue = xmin, maxvalue=xmax),
        Scale.y_continuous(minvalue = ymin, maxvalue=ymax),
        Guide.xlabel("Real"),
        Guide.ylabel("Imaginary", orientation=:vertical),
        theme,
        Guide.annotation(Compose.compose!(Compose.context(), Gadfly.fill(Colors.RGB(1,0,0)), shapes)))
end

function timeplot_immerse(t, y, ymin, ymax)
  Immerse.plot(
        x = t,
        y = y,
        Scale.y_continuous(minvalue = ymin, maxvalue=ymax),
        Guide.xlabel("Time"),
        Guide.ylabel("Output", orientation=:vertical),
        default_theme,
        Geom.line)
end

function xcross(xs::AbstractArray, ys::AbstractArray, rs::AbstractArray)
  n = max(length(xs), length(ys), length(rs))
  polys = Vector{Vector{Tuple{Measures.Measure, Measures.Measure}}}(n)
  s = 1/sqrt(5)
  for i in 1:n
    x = Compose.x_measure(xs[i])
    y = Compose.y_measure(ys[i])
    r = rs[i]
    u = s*r
    polys[i] = Tuple{Measures.Measure, Measures.Measure}[
      (x, y - u), (x + u, y - 2u), (x + 2u, y - u),
      (x + u, y), (x + 2u, y + u), (x + u, y + 2u),
      (x, y + u), (x - u, y + 2u), (x - 2u, y + u),
      (x - u, y), (x - 2u, y - u), (x - u, y - 2u)
    ]
  end
  return Gadfly.polygon(polys)
end


text:
# compute number of samples, average over samples
@everywhere function av_mapper(x)
    n = length(x)
    μ = sum(x) / n
    (n, μ)
end

# combine two (n, μ) tuples
@everywhere function av_reducer(x, y)
    n = x[1] + y[1]
    μ = (x[1] * x[2] + y[1] * y[2]) / n
    (n, μ)
end


text:
abstract type AbstractThermalAnemometer end

"""
    CTASensor(R, Rw)

A structure to manage constant temperature anemometer sensors (CTA)

 * `R` Resistance element
 * `Rw` Operating resistance
 * `gain` Output gain, ``Eo = g⋅E``

The overheat ratio is defined by the ratio between over-resistance and reference
resistance:

    `a = (Rw - Ro) / Ro`
where `a` is the overheat ratio, `Rw` is the operating resistance of the element sensor,
`Ro` is the reference resistance (resistance at reference temperature).

"""
struct CTASensor{ResType <: AbstractResistor} <: AbstractThermalAnemometer
    "Element whose resistance changes with temperature"
    R::ResType
    "Operating temperature of the sensor"
    Rw::Float64
    "Reference temperature"
    T₀::Float64
    "Output gain"
    gain::Float64

end
Base.broadcastable(sensor::CTASensor) = Ref(sensor)

CTASensor(R::ResType, Rw; gain=1.0) where {ResType <: AbstractResistor} = CTASensor(R, Rw, temperature(R), gain)

CTASensor(R::ResType, Rw, T₀; gain=1.0) where {ResType <: AbstractResistor} = CTASensor(R, Rw, T₀, gain)

Wire(R₀, Rw, T₀=293.15, α=0.36e-3, gain=1.0) = CTASensor(Resistor(R₀, α, T₀), Rw; gain=gain)

"Operating temperature of the CTA"
optemperature(w::CTASensor) = temperature(w.R, w.Rw)

"Overheat ratio of the CTA sensor"
overheat_ratio(w::CTASensor) = w.Rw/resistance(w.R) - 1.0

"Over temperature"
overtemp(w::CTASensor) = optemperature(w) - temperature(w.R)

"Resistor element of a thermal anemometer sensor"
resistor(w::AbstractThermalAnemometer) = w.R

"Gain of anemometer output"
gain(w::AbstractThermalAnemometer) = w.gain

"""
    CCASensor(R, a)

A structure to manage constant current anemometer sensors (CCA)

 * `R` Resistance element
 * `I` current in Ampere

"""
struct CCASensor{ResType <: AbstractResistor} <: AbstractThermalAnemometer
    "Element whose resistance changes with temperature"
    R::ResType
    "Operating current in Ampere"
    I::Float64
    "Output gain"
    gain::Float64
end
CCASensor(R::ResType, I; gain=1.0) where {ResType <: AbstractResistor} = CCASensor(R,I,gain)
Base.broadcastable(sensor::CCASensor) = Ref(sensor)

"Current in A flowing through the CCA"
current(w::CCASensor) = w.I


text:
#= The EP implementation is derived from:
Alfredo Braunstein, Anna Paola Muntoni, and Andrea Pagnani, “An Analytic Approximation of the Feasible
Space of Metabolic Networks,” Nature Communications 8 (April 6, 2017): https://doi.org/10.1038/ncomms14915.
Package metabolicEP (https://github.com/anna-pa-m/Metabolic-EP)
It was modified just for implement MaxEnt as described in the Appendix of:
Jorge Fernandez-de-Cossio-Diaz and Roberto Mulet, “Maximum Entropy and Population Heterogeneity in
Continuous Cell Cultures.” PLOS Computational Biology 15, no. 2 (February 27, 2019): e1006823.
https://doi.org/10.1371/journal.pcbi.1006823.
=#

module MetEP

    using SparseArrays
    using ProgressMeter

    import MetNets
    import MetNets: MetNet, AbstractMetState
    import ExtractMacro: @extract
    using LinearAlgebra
    using LinearAlgebra: inv!
    import Printf: @printf
    import SpecialFunctions: erf
    import Distributions: Truncated, Normal, mean, var, MvNormal

    import MetLP
    import MetLP: fva_preprocess

    # TODO: Add Chemostat HR code

    # Types
    include("Types/AbstractEPMat.jl")
    include("Types/EPOut.jl")
    include("Types/EPAlgs.jl")
    include("Types/EPFields.jl")
    include("Types/EPMat.jl")
    include("Types/EPMatT0.jl")
    include("Types/EPModel.jl")

    # EP
    include("EP/compute_mom5d.jl")
    include("EP/converge_ep!.jl")
    include("EP/echelonize.jl")
    include("EP/epconverge.jl")
    include("EP/eponesweep.jl")
    include("EP/eponesweepT0.jl")
    include("EP/fast_maxent_ep.jl")
    include("EP/get_scalefactor.jl")
    include("EP/matchmom.jl")
    include("EP/maxent_ep.jl")
    include("EP/newav.jl")
    include("EP/newμs.jl")
    include("EP/prepare_beta_vec.jl")
    include("EP/prepareinput.jl")
    include("EP/produce_epout.jl")
    include("EP/Q_sigma.jl")
    include("EP/scaleepfield.jl")
    include("EP/status.jl")
    include("EP/toy_maxent_ep.jl")
    include("EP/utils.jl")
    include("EP/get_join.jl")

    # Utils
    include("Utils/summary.jl")
    include("Utils/getters.jl")
    include("Utils/trunc_sample.jl")

end


text:
## Copyright (c) 2013 Miguel Bazdresch
##
## This file is distributed under the 2-clause BSD License.

# Asynchronously reads the specified IO.
# In case of timeout sends :timeout; in case of end of file, sends :eof.
function async_reader(io::IO, timeout_sec)::Channel
    ch = Channel(1)
    task = @async begin
        reader_task = current_task()
        function timeout_cb(timer)
            put!(ch, :timeout)
            Base.throwto(reader_task, InterruptException())
        end
        timeout = Timer(timeout_cb, timeout_sec)
        data = String(readavailable(io))
        if data == ""; put!(ch, :eof); return; end
        timeout_sec > 0 && close(timeout) # Cancel the timeout
        put!(ch, data)
    end
    bind(ch, task)
    return ch
end

# Write plotting data to file.
# `curve` is a Curve.
# `file` is the file to write to.
# If `append` is true, data is appended at the end of the file (for `plot!`)
# TODO: generalize to more formats
function write_data(curve, dims, file; append = false)
    mode = "w"
    append && (mode = "a")
    x = curve.x
    y = curve.y
    z = curve.z
    supp = curve.supp
    open(file, mode) do io
        # 2-d plot
        if dims == 2
            if z !== nothing
                # image; format is "x y z" with reversed "y"
                if ndims(x) == 1 && ndims(y) == 1 && ndims(z) == 2
                    xx = repeat(x,inner=length(y))
                    yy = repeat(reverse(y),length(x))
                    zz = vec(z)
                    writedlm(io,[xx yy zz])
                    # rgbimage; format is "x y r g b" with reversed "y"
                elseif ndims(z) == 3
                    xx = repeat(x,inner=length(y))
                    yy = repeat(reverse(y),length(x))
                    r = vec(z[1,:,:])
                    g = vec(z[2,:,:])
                    b = vec(z[3,:,:])
                    if supp === nothing
                        writedlm(io,[xx yy r g b])
                    else
                        writedlm(io,[xx yy r g b supp])
                    end
                end
            # regular plot
            else
                if supp === nothing
                    data = [x y]
                else
                    data = [x y supp]
                end
                writedlm(io, data)
            end
        # 3-D image
        elseif dims == 3
            # surface plot
            if ndims(x) == 1 && ndims(y) == 1 && ndims(z) == 2
                for (yi,yy) in enumerate(y)
                    for (xi,xx) in enumerate(x)
                        write(io, "$xx $yy $(z[yi,xi])\n")
                    end
                    write(io, "\n")
                end
            # scatter plot
            elseif ndims(x) == 1 && ndims(y) == 1 && ndims(z) == 1
                if supp === nothing
                    for k in 1:length(x)
                        write(io, "$(x[k]) $(y[k]) $(z[k])\n")
                    end
                else
                    for k in 1:length(x)
                        write(io, "$(x[k]) $(y[k]) $(z[k]) ")
                        s = string(supp[k,:])[2:end-1]
                        write(io, s, "\n")
                    end
                end
            # arbitrary plot
            elseif ndims(x) == 2 && ndims(y) == 2 && ndims(z) == 2
                for col = 1:size(x)[2]
                    writedlm(io, [x[:,col] y[:,col] z[:,col]])
                    write(io, "\n")
                end
            end
        end
        write(io,"\n\n")
    end
end

# llplot() is the workhorse plotting function. It handles all communication
# with gnuplot required to actually plot. It can also handle multiplot.
function llplot(F::Figure ; printstring=nothing)
    global gnuplot_state

    isempty(F) && return nothing

    config[:mode] == "null" && return nothing

    gnuplot_state.gnuplot_available || return nothing

    gnuplot_send("reset session")

    # Build terminal setup string
    if printstring === nothing
        handle = F.handle
        gnuplot_send(termstring(handle))
    else
        gnuplot_send(printstring[1])
        gnuplot_send("set output '$(printstring[2])'")
    end

    # If a layout has been defined, enter multiplot mode
    multiplot = false
    if F.layout != (0,0)
        multiplot = true
        gnuplot_send("set multiplot layout $(F.layout[1]),$(F.layout[2]) columnsfirst")
    end

    for sp in F.subplots
        if multiplot
            if (sp === nothing) || (isempty(sp))
                gnuplot_send("set multiplot next")
                continue
            end
        end
        # Send preamble
        prm = config[:preamble]
        if !isempty(prm)
            gnuplot_send(prm)
        end
        # Build figure configuration string
        gnuplot_send(sp.axesconf)
        # send plot command to gnuplot
        gnuplot_send(plotstring(sp))
        # clear session for next plot
        if multiplot
            gnuplot_send("reset session")
        end
    end

    if multiplot
        gnuplot_send("unset multiplot")
    end

    # Close output files, if any
    gnuplot_send("set output")

    # Make sure gnuplot is done.
    err = ""
    gnuplot_state.gp_lasterror = err
    gnuplot_state.gp_error = false

    gnuplot_send("""set print '-'
                    print 'GastonDone'""")

    # Start reading gnuplot's streams in "background"
    ch_out = async_reader(P.gstdout, config[:timeout])
    out = take!(ch_out)
    out === :timeout && @warn("Gnuplot is taking too long to respond.")
    out === :eof     && error("Gnuplot crashed")

    # check for errors while plotting
    yield()
    if bytesavailable(P.gstderr) > 0
        err = String(readavailable(P.gstderr))
        gnuplot_state.gp_lasterror = err
        gnuplot_state.gp_error = true
        @warn("Gnuplot returned an error message:\n  $err")
    end

    return nothing

end


text:
const MAGIC = (UInt8('O'), UInt8('b'), UInt8('j'), 0x01)
const FileHeader = Record{(:magic, :meta, :sync), Tuple{NTuple{4, UInt8}, Dict{String, Vector{UInt8}}, NTuple{16, UInt8}}, 3}
const FileHeaderRecordType = schematype(FileHeader)

struct Block
    count::Int64
    bytes::SubArray{UInt8, 1, Vector{UInt8}, Tuple{UnitRange{Int64}}, true}
    sync::NTuple{16, UInt8}
end

StructTypes.StructType(::Type{Block}) = StructTypes.Struct()
const BlockType = schematype(Block)

"""
    Avro.writetable(io_or_file, tbl; kw...)

Write an input Tables.jl-compatible source table `tbl` out as an avro
object container file. `io_or_file` can be a file name as a `String` or
`IO` argument. If the input table supports `Table.partitions`, each
partition will be written as a separate "block" in the container file.

Because avro data is strictly typed, if the input table doesn't have a
well-defined schema (i.e. `Tables.schema(Tables.rows(tbl)) === nothing`),
then `Tables.dictrowtable(Tables.rows(tbl))` will be called, which scans
the input table, "building up" the schema based on types of values found
in each row.

Compression is supported via the `compress` keyword argument, and can
currently be one of `:zstd`, `:deflate`, `:bzip2`, or `:xz`.
"""
function writetable end

writetable(io_or_file; kw...) = x -> writetable(io_or_file, x; kw...)

function writetable(file::String, tbl; kw...)
    open(file, "w") do io
        writetable(io, tbl; kw...)
    end
    return file
end

function writetable(io::IO, source; compress::Union{Nothing, Symbol}=nothing, kw...)
    comp = get(COMPRESSORS, compress, nothing)
    parts = Tables.partitions(source)
    state = iterate(parts)
    state === nothing && error("no data in input; unable to write avro file")
    part, st = state
    rows = Tables.rows(part)
    sch = Tables.schema(rows)
    dictrow = false
    if sch === nothing || !Base.haslength(rows)
        rows = Tables.dictrowtable(rows)
        sch = Tables.schema(rows)
        dictrow = true
    end
    writewithschema(io, parts, rows, st, sch, dictrow, comp, kw)
    return io
end

function writewithschema(io, parts, rows, st, sch, dictrow, comp, kw)
    schtyp = schematype(sch)
    meta = Dict("avro.schema" => JSON3.write(schtyp))
    sync = _cast(NTuple{16, UInt8}, rand(UInt128))
    buf = write((magic=MAGIC, meta=meta, sync=sync); schema=FileHeaderRecordType)
    Base.write(io, buf)
    @debug 1 "wrote file header from bytes 1:$(pos - 1)"
    i = 1
    while true
        # if rows didn't have schema or length, we materialized w/ Tables.dictrowtable
        nrow = length(rows)
        @debug 1 "writing block count ($nrow) at pos = $pos"
        rowsstate = iterate(rows)
        pos = 1
        if rowsstate === nothing
            bytes = UInt8[]
            pos = 0
        else
            row, rowst = rowsstate
            # calc nbytes on first row, then allocate bytes
            bytesperrow = nbytes(schtyp, row)
            blen = trunc(Int, nrow * bytesperrow * 1.05) # add 5% cushion
            bytes = Vector{UInt8}(undef, blen)
            n = 1
            while true
                pos = writevalue(Binary(), schtyp, row, bytes, pos, blen, kw)
                rowsstate = iterate(rows, rowst)
                rowsstate === nothing && break
                row, rowst = rowsstate
                nb = nbytes(schtyp, row)
                if (pos + nb - 1) > blen
                    # unlikely, but need to resize our buffer
                    rowslefttowrite = nrow - n
                    avgbytesperrowwritten = div(bytesperrow, n)
                    resize!(bytes, pos + (rowslefttowrite * avgbytesperrowwritten))
                end
                bytesperrow += nb
                n += 1
            end
        end
        # compress
        if comp !== nothing
            finalbytes = transcode(comp[Threads.threadid()], unsafe_wrap(Base.Array, pointer(bytes), pos - 1))
        else
            finalbytes = bytes
        end
        block = Block(nrow, view(bytes, 1:(pos - 1)), sync)
        buf = write(block; schema=BlockType)
        Base.write(io, buf)
        state = iterate(parts, st)
        state === nothing && break
        part, st = state
        rows = Tables.rows(part)
        sch = Tables.schema(rows)
        if dictrow
            rows = Tables.dictrowtable(rows)
        end
    end
    return
end

"""
    Avro.Table

A Tables.jl-compatible source returned from `Avro.readtable`.
Conceptually, it can be thought of as an `AbstractVector{Record}`,
where `Record` is the avro version of a "row" or NamedTuple. Thus,
`Avro.Table` supports indexing/iteration like an `AbstractVector`.
"""
struct Table{T} <: AbstractVector{T}
    sch::Tables.Schema
    data::ChainedVector{T}
end

data(x::Table) = getfield(x, :data)
Base.IndexStyle(::Type{Table}) = Base.IndexLinear()
Base.@propagate_inbounds Base.getindex(x::Table, i::Int) = data(x)[i]
Base.size(x::Table) = size(data(x))
Base.IteratorSize(::Type{Table}) = Base.HasLength()
Base.IteratorEltype(::Type{Table}) = Base.HasEltype()
Base.iterate(x::Table) = iterate(data(x))
Base.iterate(x::Table, st) = iterate(data(x), st)

Tables.isrowtable(::Type{Table}) = true
Tables.schema(x::Table) = getfield(x, :sch)

"""
    Avro.readtable(file_or_io) => Avro.Table

Read an avro object container file, returning an [`Avro.Table`](@ref) type,
which is like an array of records, where each record follows the
schema encoded with the file. Any compression will be detected and
decompressed automatically when reading.
"""
function readtable end

readtable(io::IO; kw...) = readtable(Base.read(io), 1, kw)
readtable(io::IOBuffer; kw...) = readtable(take!(io), 1, kw)
readtable(file::String; kw...) = readtable(Mmap.mmap(file), 1, kw)

function readtable(buf, pos, kw)
    header, pos = readvalue(Binary(), FileHeaderRecordType, FileHeader, buf, pos, length(buf), nothing)
    sch = JSON3.read(header.meta["avro.schema"], Schema)
    comp = get(header.meta, "avro.codec", nothing)
    comp = get(DECOMPRESSORS, Symbol(comp), nothing)
    T = juliatype(sch)
    data = readwithschema(T, sch, buf, pos, comp)
    return Table(Tables.Schema(T), data)
end

function readwithschema(::Type{T}, sch, buf, pos, comp) where {T}
    len = length(buf)
    blocks = Array{T}[]
    while pos <= len
        block, pos = readvalue(Binary(), BlockType, Block, buf, pos, len, nothing)
        @debug 1 "reading block: count = $(block.count), size = $(length(block.bytes))"
        values = Vector{Int}(undef, block.count)
        bytes = block.bytes
        # uncompress
        if comp !== nothing
            bytes = transcode(comp, bytes)
        end
        bpos = 1
        blen = length(bytes)
        for i = 1:block.count
            @inbounds values[i] = bpos
            bpos = skipvalue(Binary(), sch, T, bytes, bpos, blen, nothing)
        end
        push!(blocks, Array{T}(bytes, values))
    end
    return ChainedVector(blocks)
end


text:
# MIT license
# Copyright (c) Microsoft Corporation. All rights reserved.
# See LICENSE in the project root for full license information.

# FIXING ERRONEOUS UNBOUND TYPE ERRORS THAT OCCUR WITH VARARG
#=
The set will contain something like this:
    MultiHologramInterface(interfaces::Vararg{HologramInterface{T}, N}) where {T<:Real, N}
To get the signature run:
    methods(OpticSim.MultiHologramInterface).ms[I].sig            where I is typically 1 or 2 depending on the number of methods
This gives:
    Tuple{Type{MultiHologramInterface}, Vararg{HologramInterface{T}, N}} where N where T<:Real
We then have to use which to find the method:
    which(OpticSim.MultiHologramInterface, Tuple{Vararg{HologramInterface{T}, N}} where N where T<:Real)
and pop it from the set
=#

@testset "JuliaLang" begin
    # here we ignore the specific methods which we know are ok but are still failing
    # for some weird reason the unbound args check seems to fail for some (seemingly random) methods with Vararg arguments
    methods_to_ignore = Dict(
        OpticSim => VERSION >= v"1.6.0-DEV" ? [
            (OpticSim.LensAssembly, Tuple{Vararg{Union{CSGTree{T},LensAssembly{T},OpticSim.Surface{T}},N} where N} where {T<:Real}),
            (OpticSim.MultiHologramInterface, Tuple{Vararg{HologramInterface{T},N}} where {N} where {T<:Real}),
        ] : [],
        OpticSim.Zernike => [],
        OpticSim.QType => [],
        OpticSim.Vis => [
            (OpticSim.Vis.drawcurves, Tuple{Vararg{Spline{P,S,N,M},N1} where N1} where {M} where {N} where {S} where {P}),
            (OpticSim.Vis.draw, Tuple{Vararg{S,N} where N} where {S<:Union{OpticSim.Surface{T},OpticSim.TriangleMesh{T}}} where {T<:Real}),
            (OpticSim.Vis.draw!, Tuple{OpticSim.Vis.Makie.LScene,Vararg{S,N} where N} where {S<:Union{OpticSim.Surface{T},OpticSim.TriangleMesh{T}}} where {T<:Real})
        ],
        OpticSim.Examples => [],
        OpticSim.Chebyshev => [],
        OpticSim.GlassCat => [],
        OpticSim.Optimization => [],
    )

    for (mod, ignore_list) in methods_to_ignore
        unbound = setdiff(
            detect_unbound_args(mod),
            [which(method, signature) for (method, signature) in ignore_list]
        )
        # also ignore any generate methods created due to default or keyword arguments
        filter!(m -> !occursin("#", string(m.name)), unbound)
        @test isempty(unbound)

        ambiguous = detect_ambiguities(mod)
        @test isempty(ambiguous)
    end
end # testset JuliaLang


text:
function add(::Type{FunctionUtils}, arg0::Vector{UnivariateDifferentiableFunction})
    return jcall(FunctionUtils, "add", UnivariateDifferentiableFunction, (Vector{UnivariateDifferentiableFunction},), arg0)
end

function add(::Type{FunctionUtils}, arg0::Vector{UnivariateFunction})
    return jcall(FunctionUtils, "add", UnivariateFunction, (Vector{UnivariateFunction},), arg0)
end

function collector(::Type{FunctionUtils}, arg0::BivariateFunction, arg1::UnivariateFunction, arg2::jdouble)
    return jcall(FunctionUtils, "collector", MultivariateFunction, (BivariateFunction, UnivariateFunction, jdouble), arg0, arg1, arg2)
end

function collector(::Type{FunctionUtils}, arg0::BivariateFunction, arg1::jdouble)
    return jcall(FunctionUtils, "collector", MultivariateFunction, (BivariateFunction, jdouble), arg0, arg1)
end

function combine(::Type{FunctionUtils}, arg0::BivariateFunction, arg1::UnivariateFunction, arg2::UnivariateFunction)
    return jcall(FunctionUtils, "combine", UnivariateFunction, (BivariateFunction, UnivariateFunction, UnivariateFunction), arg0, arg1, arg2)
end

function compose(::Type{FunctionUtils}, arg0::Vector{UnivariateDifferentiableFunction})
    return jcall(FunctionUtils, "compose", UnivariateDifferentiableFunction, (Vector{UnivariateDifferentiableFunction},), arg0)
end

function compose(::Type{FunctionUtils}, arg0::Vector{UnivariateFunction})
    return jcall(FunctionUtils, "compose", UnivariateFunction, (Vector{UnivariateFunction},), arg0)
end

function derivative(::Type{FunctionUtils}, arg0::MultivariateDifferentiableFunction, arg1::Vector{jint})
    return jcall(FunctionUtils, "derivative", MultivariateFunction, (MultivariateDifferentiableFunction, Vector{jint}), arg0, arg1)
end

function derivative(::Type{FunctionUtils}, arg0::UnivariateDifferentiableFunction, arg1::jint)
    return jcall(FunctionUtils, "derivative", UnivariateFunction, (UnivariateDifferentiableFunction, jint), arg0, arg1)
end

function equals(obj::Object, arg0::Object)
    return jcall(obj, "equals", jboolean, (Object,), arg0)
end

function fix1st_argument(::Type{FunctionUtils}, arg0::BivariateFunction, arg1::jdouble)
    return jcall(FunctionUtils, "fix1stArgument", UnivariateFunction, (BivariateFunction, jdouble), arg0, arg1)
end

function fix2nd_argument(::Type{FunctionUtils}, arg0::BivariateFunction, arg1::jdouble)
    return jcall(FunctionUtils, "fix2ndArgument", UnivariateFunction, (BivariateFunction, jdouble), arg0, arg1)
end

function get_class(obj::Object)
    return jcall(obj, "getClass", Class, ())
end

function hash_code(obj::Object)
    return jcall(obj, "hashCode", jint, ())
end

function multiply(::Type{FunctionUtils}, arg0::Vector{UnivariateDifferentiableFunction})
    return jcall(FunctionUtils, "multiply", UnivariateDifferentiableFunction, (Vector{UnivariateDifferentiableFunction},), arg0)
end

function multiply(::Type{FunctionUtils}, arg0::Vector{UnivariateFunction})
    return jcall(FunctionUtils, "multiply", UnivariateFunction, (Vector{UnivariateFunction},), arg0)
end

function notify(obj::Object)
    return jcall(obj, "notify", void, ())
end

function notify_all(obj::Object)
    return jcall(obj, "notifyAll", void, ())
end

function sample(::Type{FunctionUtils}, arg0::UnivariateFunction, arg1::jdouble, arg2::jdouble, arg3::jint)
    return jcall(FunctionUtils, "sample", Vector{jdouble}, (UnivariateFunction, jdouble, jdouble, jint), arg0, arg1, arg2, arg3)
end

function to_differentiable(::Type{FunctionUtils}, arg0::MultivariateFunction, arg1::MultivariateVectorFunction)
    return jcall(FunctionUtils, "toDifferentiable", MultivariateDifferentiableFunction, (MultivariateFunction, MultivariateVectorFunction), arg0, arg1)
end

function to_differentiable(::Type{FunctionUtils}, arg0::UnivariateFunction, arg1::Vector{UnivariateFunction})
    return jcall(FunctionUtils, "toDifferentiable", UnivariateDifferentiableFunction, (UnivariateFunction, Vector{UnivariateFunction}), arg0, arg1)
end

function to_string(obj::Object)
    return jcall(obj, "toString", JString, ())
end

function wait(obj::Object)
    return jcall(obj, "wait", void, ())
end

function wait(obj::Object, arg0::jlong)
    return jcall(obj, "wait", void, (jlong,), arg0)
end

function wait(obj::Object, arg0::jlong, arg1::jint)
    return jcall(obj, "wait", void, (jlong, jint), arg0, arg1)
end


text:
# This file is a part of TypeDBClient.  License is MIT: https://github.com/Humans-of-Julia/TypeDBClient.jl/blob/main/LICENSE

mutable struct LogicManager <: AbstractLogicManager
    transaction::AbstractCoreTransaction
end

function get_rule(log_mgr::AbstractLogicManager, label::String)
    res = execute(log_mgr.transaction, LogicManagerRequestBuilder.get_rule_req(label))
    result = if which_oneof(res.logic_manager_res.get_rule_res, :res) == :rule
                Rule(res.logic_manager_res.get_rule_res.rule)
             else
                nothing
             end
    return result
end

function get_rules(log_mgr::AbstractLogicManager)
    dbResult = stream(log_mgr.transaction, LogicManagerRequestBuilder.get_rules_req())
    answers = [res.logic_manager_res_part.get_rules_res_part.rules for res in dbResult]
    rules = map(Rule,answers)
    result = reduce(vcat, rules)
    return result
end

function put_rule(log_mgr::AbstractLogicManager, label::String, when::String, then::String)
    res = execute(log_mgr.transaction, LogicManagerRequestBuilder.put_rule_req(label, when, then))
    result = Rule(res.logic_manager_res.get_rule_res.rule)
    return result
end


text:
# I map some model ids to the enzimatic cost data reported in
# Beg2007, so, I will set the cost of this reactions using this data
# The rest of the model reactions will get the average enzimatic cost also
# reported by Beg.

function load_beg_enz_cost_dict()
    beg_enz_cost_dict = Dict()
    beg_rxns_map = load_beg_rxns_map()
    beg_enz_data = Bd.load_enz_data()
    for (i, beg_rxn) in enumerate(beg_enz_data.Enzyme)
        for model_rnx in beg_rxns_map[beg_rxn]
            if model_rnx != ""
                beg_enz_cost_dict[model_rnx] =
                    beg_enz_data[i, Symbol("ai/xi (h g / mmol)")]
            end
        end
    end
    beg_enz_cost_dict
end

function beg_enz_cost(model_rxn)
    beg_enz_cost_dict = load_beg_enz_cost_dict()
    return get(beg_enz_cost_dict, model_rxn, Bd.AVE_COST)
end


text:
using PastaQ
using ITensors
using Random
using Test

@testset "measure one-body operators" begin
  N = 10
  depth = 10
  ψ = runcircuit(randomcircuit(N, depth))

  @test measure(ψ, ("X", 2)) ≈ inner(ψ, runcircuit(ψ, ("X", 2)))
  @test measure(ψ, ("Y", 1)) ≈ inner(ψ, runcircuit(ψ, ("Y", 1)))
  @test measure(ψ, ("Z", 4)) ≈ inner(ψ, runcircuit(ψ, ("Z", 4)))

  results = measure(ψ, "X")
  for j in 1:length(ψ)
    @test results[j] ≈ inner(ψ, runcircuit(ψ, ("X", j)))
  end

  results = measure(ψ, ("X", 1:2:length(ψ)))
  for j in 1:2:length(ψ)
    @test results[(j + 1) ÷ 2] ≈ inner(ψ, runcircuit(ψ, ("X", j)))
  end

  results = measure(ψ, ("Y", [1, 3, 5]))
  for j in [1, 3, 5]
    @test results[(j + 1) ÷ 2] ≈ inner(ψ, runcircuit(ψ, ("Y", j)))
  end
end


text:
using Plots
using DrWatson
using DataFrames
using CSV

fn = plotsdir("summary_1.csv")
df = CSV.File(fn) |> DataFrame

n=[9437186,37748738]

df_n1=df[df.ngdofs .== n[1], :]
df_n2=df[df.ngdofs .== n[2], :][1:5,:]

x1 = df_n1[:,:np]
y1 = df_n1[:,:Simulation]

x2 = df_n2[:,:np]
y2 = df_n2[:,:Simulation]

function doplot()
  plt = plot(x1,y1,
     thickness_scaling = 1.2,
     xaxis=("Number of cores",:log),
     yaxis=("Wall time [s]",:log),
     xticks=([96,192,384,672,960],
     ["96","192","384","672","960"]),
     yticks=([66,78,87,149,185,200,253,280,336,377],
     ["66","78","87","149","185","200","253","280","336","377"]),
     minorgrid=true,
     title="Galewsky. 10x steps (θ=0.5,dt=480[s]). Newton/GMRES-RAS.",
     titlefont=font(8),
     shape=:square,
     label="256x256x6 [9.4MDoFs]",
     legend=:bottomleft,
     legendfont=font(6),
     markerstrokecolor=:white,
    )

  plot!(x2,y2,
    thickness_scaling = 1.2,
    xaxis=("Number of cores",:log),
    yaxis=("Wall time [s]",:log),
    xticks=([96,192,384,672,960],
    ["96","192","384","672","960"]),
    yticks=([66,78,87,149,185,200,253,280,336,377],
    ["66","78","87","149","185","200","253","280","336","377"]),
    minorgrid=true,
    title="Galewsky. 10x time steps (θ=0.5,dt=480[s]). Newton/GMRES-RAS.",
    titlefont=font(8),
    shape=:square,
    label="512x512x6 [37.7MDoFs]",
    legend=:bottomleft,
    legendfont=font(6),
    color=:green,
    markerstrokecolor=:white,
   )

  function plot_vert(plt,x,y,ngdofs,l,p,color)
     x1 = first(x)
     y1 = first(y)
     s = y1.*(x1./x)
     println(x)
     println(s)
     plot!(x,s,xaxis=:log,yaxis=:log, color=color, label="Ideal"),
     Plots.vline!([l],
           xaxis=:log,
           yaxis=:log,
           linestyle=:dot,seriestype = :vline,
           color=color,
           label="$(p)KDoFs/core")
     plt
  end
  plot_vert(plt,x1,y1,first(df_n1.ngdofs),first(df_n1.ngdofs)/25e3,25,:blue)
  plot_vert(plt,x2,y2,first(df_n2.ngdofs),first(df_n2.ngdofs)/50e3,50,:green)
  plt
end

#unicodeplots()
#display(doplot())

gr()
doplot()
savefig(plotsdir("total_scaling.pdf"))
savefig(plotsdir("total_scaling.png"))

max_stars_repo_name: santiagobadia/GridapGeosciences

max_issues_repo_name: santiagobadia/GridapGeosciences

max_forks_repo_name: santiagobadia/GridapGeosciences


text:
"""
Motif-based community detection using a novel modification of the Louvain
algorithm for modularity optimization.

Based on "Motif-based communities in complex networks" by Alex Arenas,
Alberto Fernandez, Santo Fortunato, and Sergio Gomez. The original paper may
be found at https://arxiv.org/pdf/0710.0059.pdf.
"""

using MatrixNetworks
include("Utilities.jl")
include("Cliques.jl")

"""
Common motif generators
    - `path(k::Int64)`: Returns a path of length `k`
    - `cycle(k::Int64)`: Returns a cycle of length `k`
    - `clique(k::Int64)`: Returns a clique of size `k`
"""
path(k::Int64) = MatrixNetwork([i for i = 1:k], [j+1 for j = 1:k])
cycle(k::Int64) = MatrixNetwork([i for i = 1:k], [[j+1 for j = 1:k-1]; [1]])
clique(k::Int64) = MatrixNetwork([i for i = 1:k for j = 1:k-1], [j for i = 1:k for j = 1:k if j != i])

"""
`motif_modularity`
==================

Calculate the motif modularity of a graph given a motif and community partition.
If no motif is provided, the calculated quantity is the traditional modularity.

Arguments
---------
    - `A::MatrixNetwork`: The graph to calculate motif modularity for
    - `C::Vector{Vector{Int64}}`: The community partition of `A`. The ith subarray of `C` contains the nodes of the ith community.
    - `M::MatrixNetwork (=MatrixNetwork([1],[2]))`: The given motif. A directed edge by default, will measure the standard modularity.

Modularity is a measure of how "good" a community partition is at dividing the
network into a set of tightly intraconnected groups with few interconnections.
Typically, intra- and interconnection are quantified by the number of edges
within and between communities respectively. However, as the author mentions on
pp. 4 of the paper, we can use __any__ connected pattern (weakly connected for a digraph),
called a motif, to measure connection strength by counting the number of motif
occurrences within and between communities.
"""
function motif_modularity(A::MatrixNetwork, C::Vector{Vector{Int64}}; M::MatrixNetwork=MatrixNetwork([1],[2]))
    node_communities = zeros(Int64, A.n) # Easy lookup for which community a node is part of
    for i = 1:size(C,1)
        for n in C[i]
            node_communities[n] = i
        end
    end

    k_delt(n1, n2) = return node_communities[n1] == node_communities[n2] # Kronecker delta function, returns true if nodes are in same community

    a_e = get_edges(A) # Network edges
    m_e = get_edges(M) # Motif edges
    out_w = outweight_distr(A)
    in_w = inweight_distr(A)
    weights = total_edge_weights(A) # Total weights between each pair of adjacent nodes

    shared_strengths = Dict() # Relates sets of communities to the motif strength among them.
                              # E.g., if communities 1,2,3 share a triangle, then
                              # we keep track of that. This way, we can update modularity
                              # without recomputing from scratch.

    function inc_inds!(inds::Vector{Int64}) # Simulate one iteration of an arbitrary number of nested for loops
        an = A.n
        sz = size(inds,1)
        if inds == an * ones(Int64, sz)
            inds .-= inds
            return
        end

        for i = sz:-1:1
            if inds[i] == an
                inds[i] = 1
                continue
            else
                inds[i] += 1
                return
            end
        end
        if length(inds) != length(Set(inds)) # Requires unique indices, avoids "degenerate" motifs
            inc_inds!(inds)
        end
    end

    indices = ones(Int64, M.n)
    z = zeros(Int64, M.n)
    w,W,u,U = 0,0,0,0 # w: strength of motifs fully contained within communities, W: total strength of motifs
                      # u: w for a random, strength-preserving graph, U: u for the same

    while indices != z # While all groups of M.n nodes haven't been checked for containing a motif occurrence
        w_inc = 1 # Controls whether strength gets added to w
        W_inc = 1 # Similar, but for W
        weight = 1

        p = 1 # Amount to add to u and U
        nf = 1 # Similar to w_inc, but for u

        motif_membership = Set() # Communities that the motif is part of

        for e in m_e
            if !(indices[e] in a_e) # Motif edge not in graph
                w_inc = 0
                W_inc = 0
            end
            if !k_delt(indices[e]...) # Edge crosses community lines
                w_inc = 0
                nf = 0
                push!(motif_membership, node_communities[indices[e[1]]]) # Motif shared between these two communities
                push!(motif_membership, node_communities[indices[e[2]]])
            end
            if w_inc == 1
                weight *= weights[[indices[e[1]], indices[e[2]]]]
            end
            p *= out_w[indices[e[1]]] * in_w[indices[e[2]]]
        end
        w += w_inc * weight
        W += W_inc * weight
        u += p * nf
        U += p

        # Check that the group has enough nodes to be a motif
        if !isempty(motif_membership) && reduce(+, length.(C[collect(motif_membership)])) >= M.n
            motif_membership = sort(collect(motif_membership))
            if !haskey(shared_strengths, motif_membership)
                shared_strengths[motif_membership] = [0,0]
            end
            shared_strengths[motif_membership][1] += W_inc * weight
            shared_strengths[motif_membership][2] += p
        end

        inc_inds!(indices)
    end

    return (w/W) - (u/U), shared_strengths, W, U
end

"""
`louvain_motif`
===============

Determine an approximation for the best community division using a modified
Louvain algorithm (https://en.wikipedia.org/wiki/Louvain_modularity) on a given motif.
If no motif provided, performs the traditional Louvain algorithm.

Arguments
---------
    - `A::MatrixNetwork`: The graph in which to find communities
    - `M::MatrixNetwork (=MatrixNetwork([1],[2]))`: The motif that forms the basis of a community

The Louvain algorithm finds a local maximum for modularity by first placing
each node in its own community and then, for each community, performing the merge
that most increases modularity. Now, with a new set of communities, it repeats
this second step until no merge increases modularity.

In this modified version of the Louvain algorithm, instead of merging __pairs__
of communities that share edges, we merge __groups__ of communities that share
motifs. In a very similar process to the above Louvain algorithm, we find the
group of communities that, if combined into one, would most increase
motif modularity as defined above.
"""

function louvain_motif(A::MatrixNetwork; M::MatrixNetwork=MatrixNetwork([1],[2]))
    C = [[i] for i = 1:A.n] # Community partition
    C_last = deepcopy(C) # Copy of above, used to check if any merges occurred in previous round

    mod, strengths, W, U = motif_modularity(A,C; M=M) # Initial values for modularity shortcut

    ### Collapse all communities in the group into one. Updates C and strengths ###

    # E.g. merge_community_group!([1,2,3]) merges 2 and 3 into 1 and accordingly
    # updates the strengths table so that any groups with connections to 2 or 3
    # now have a connection to 1. Additionally, the strengths of 2 and 3 get added
    # to 1, since they are all one community now.

    function merge_community_group!(group)
        new_groups = collect(deepcopy(keys(strengths)))
        cpy = deepcopy(new_groups) # Remember old keys so we can get rid of ones that no longer exist at the end

        new_strengths = Dict()

        i = group[1]

        for k = 1:size(new_groups,1)
            g = new_groups[k]
            for c in group
                g[g.==c] .= i # Replace all other group members with i
            end

            g = sort(collect(Set(g))) # If i already included, remove duplicates and sort for consistency

            if length(g) < 2 continue end # Leave out degenerate groups

            if !haskey(new_strengths, g)
                new_strengths[g] = strengths[cpy[k]]
            else
                new_strengths[g] += strengths[cpy[k]]
            end
        end

        for g in setdiff(cpy, new_groups) # Erase entries that no longer exist
            delete!(new_strengths,g)
        end

        old_keys = collect(deepcopy(keys(new_strengths)))
        new_keys = []
        vals = []

        new_labels = collect(1:length(C))
        for c in group[2:length(group)]
            new_labels[c:length(new_labels)] .-= 1
        end

        for k = 1:length(old_keys) # Relabel community relationships to account for merge
            key = old_keys[k]
            push!(vals, new_strengths[key])
            key = [new_labels[k] for k in key]
            push!(new_keys, key)
        end

        # Update partition
        C[i] = sort(reduce(vcat, C[group]))
        for c in group[2:length(group)] setdiff!(C, [C_last[c]]) end

        strengths = Dict(new_keys[i] => vals[i] for i = 1:length(new_keys)) # Ta da!
    end

    ### Perform one iteration of the modified Louvain algorithm. Finds the group
    ### of communities that, if collapsed into one, would most increase modularity
    ### and performs the merge. ###
    function louvain_iter!()
        max_dQ = 0.0 # The merge must increase modularity, so we set the baseline to 0
        best = -1
        str_k = filter(x -> strengths[x][1] > 0, collect(keys(strengths))) # To increase modularity, strength between communities must be positive
        for k = 1:length(str_k)
            dQ = (strengths[str_k[k]][1]/W) - (strengths[str_k[k]][2]/U)
            if dQ > max_dQ
                max_dQ = dQ
                best = k
            end
        end
        if best != -1 merge_community_group!(str_k[best]) end
        mod += max_dQ
    end
    louvain_iter!() # Get the ball rolling, make C different from C_last

    while C_last != C # Stop merging when no merge increases modularity
        C_last = deepcopy(C)
        louvain_iter!()
    end

    return sort.(C)
end


text:
function em(y;nx::Int=2,u::Array=zeros(2),init=false,ssmi=SSM0(),max_iter::Int=100,tol=1e-6,txo=true,subid=false)
	"""
	estimate A, B, C, D, Q, R, S, m1, P1 using the STATIONARY EM algorithm for model

	x_{t+1} = A*x_{t} + B*u_{t} + w_{t}
	y_{t}   = C*x_{t} + D*u_{t} + v_{t}

	cov(w_{t},v_{t}) = [Q S;S' R]
	x1 -> N(m1,P1)

	Input
	----------------------
		y          : data
		nx         : size of the estimated A matrix
		u != 0			 : input are considered
		init				: if init, ssmi is considered
		ssmi       : initial point for the EM algorithm
		subid=true : SUBspace IDentification algorithm is used, not the EM
		S=true     : matrix S is estimated

	"""

	# output declaration
	ssm = SSM0()

	# initial point
	# ***************************************************************
	if init
		Ai = ssmi.A
		Bi = ssmi.B
		Ci = ssmi.C
		Di = ssmi.D
		Qi = ssmi.Q
		Ri = ssmi.R
		m1i = ssmi.m1
	else
		# initial point with subspace algorithm
		i = nx+1
		Ai,Ci,Qi,Ri,Si = ACQRS_ssi1(y,nx,i)
		m1i = Ci\y[:,1] # x11 = C^{-1}(y_1 - v_1)

		if u != zeros(2)
			# initial values for B and D
			y,ny,nt = byrow(y)
			u,nu,nt = byrow(u)
			Bi = zeros(nx,nu)
			Di = zeros(ny,nu)
		end
	end

	# estimation of parameters
	# *************************************************************
	if subid # subspace algorithm

		ssm.A = Ai
		ssm.C = Ci
		ssm.Q = Qi
		ssm.R = Ri
		ssm.S = Si
		ssm.m1 = m1i

		# loglik
		xtt,Ptt,xtt1,Ptt1,et,St,Kt,loglik = ACQR_kfilter_s(y,Ai,Ci,Qi,Ri,m1i) # change to kfilter with S
		ssm.loglik = loglik
		# aic
		ny = size(Ci,1)
		P = ny*nx + nx*nx + nx*(nx+1)/2 + ny*(ny+1)/2
		ssm.aic = -2*loglik + 2*P

	else # em algorithm

		if u == zeros(2) # em without inputs

			# em
			A,C,Q,R,m1,P1,loglikv,aic = ACQR_em_s(y,nx,Ai,Ci,Qi,Ri,m1i,max_iter,tol,txo)

			# output
			ssm.A = A
			ssm.C = C
			ssm.Q = Q
			ssm.R = R
			ssm.m1 = m1
			ssm.P1 = P1
			ssm.loglik = loglikv[end]
			ssm.aic = aic

		else # em with inputs

			# em
			A,B,C,D,Q,R,m1,P1,loglikv,aic = ABCDQR_em_s(y,u,nx,Ai,Bi,Ci,Di,Qi,Ri,m1i,max_iter,tol,txo)

			# output
			ssm.A = A
			ssm.B = B
			ssm.C = C
			ssm.D = D
			ssm.Q = Q
			ssm.R = R
			ssm.m1 = m1
			ssm.P1 = P1
			ssm.loglik = loglikv[end]
			ssm.aic = aic

		end
	end

	# output
	return ssm

end


text:
# Source comes from git
# Input is URL to git + possible tag, commit, ...
immutable type GitSource <: AbstractSource
    url::String
    commit::String
    GitSource(u::String, c::String="master") = new(u, c)
end
GitSource(u::String, v::VersionNumber) = GitSource(u, "v$v")

function download(package::AbstractPackage, gitter::GitSource, workdir::String; kwargs...)
  const source = sourcedir(package, workdir)
  mkpath(source)
  isdir(joinpath(source, ".git")) || run(`git clone $(gitter.url) $source`)
  cd(source) do
    if length(gitter.commit) > 0
      run(`git checkout $(gitter.commit)`)
    end
  end
end

max_stars_repo_name: UCL/DebbyPacker

max_issues_repo_name: UCL/DebbyPacker

max_forks_repo_name: UCL/DebbyPacker


text:
export play, scene_and_node

using Makie
using Colors

const SNAKE_GAME_SCEENS = IdDict()

function scene_and_node(game::SnakeGame)
    if haskey(SNAKE_GAME_SCEENS, game)
        SNAKE_GAME_SCEENS[game]
    else
        node = Node(game)
        scene = init_screen(node)
        get!(SNAKE_GAME_SCEENS, game, (scene, node))
    end
end

function Base.display(game::SnakeGame)
    scene, node = scene_and_node(game)
    node[] = game
    display(scene)
end

function init_screen(game::Observable{<:SnakeGame{2}}; resolution=(1000,1000))
    SNAKE_COLORS = range(HSV(60,1,1), stop=HSV(300,1,1), length=length(game[].snakes)+1)
    scene = Scene(resolution = resolution, raw = true, camera = campixel!)

    area = scene.px_area
    poly!(scene, area)

    grid_size = @lift((widths($area)[1] / size($game)[1], widths($area)[2] / size($game)[2]))

    walls = get_walls(game[])
    if length(walls) > 0
        wall_boxes = @lift([FRect2D((w.I .- (1,1)) .* $grid_size , $grid_size) for w in walls])
        poly!(scene, wall_boxes, color=:gray)
    end

    for i in 1:length(game[].snakes)
        snake_boxes = @lift([FRect2D((p.I .- (1,1)) .* $grid_size , $grid_size) for p in $game.snakes[i]])
        poly!(scene, snake_boxes, color=SNAKE_COLORS[i], strokewidth = 5, strokecolor = :black)

        snake_head_box = @lift(FRect2D(($game.snakes[i][1].I .- (1,1)) .* $grid_size , $grid_size))
        poly!(scene, snake_head_box, color=:black)
        snake_head = @lift((($game.snakes[i][1].I .- 0.5) .* $grid_size))
        scatter!(scene, snake_head, marker='◉', color=SNAKE_COLORS[i], markersize=@lift(minimum($grid_size)))
    end

    food_position = @lift([(f.I .- (0.5, 0.5)) .* $grid_size for f in $game.foods])
    scatter!(scene, food_position, color=:red, marker='♥', markersize=@lift(minimum($grid_size)))

    display(scene)
    scene
end

play() = play(SnakeGame())

function play(game::SnakeGame{2};f_name="test.gif",framerate = 2)
    @assert length(game.snakes) <= 3 "At most three players are supported in interactive mode"
    scene, game_node = scene_and_node(game)

    LEFT = CartesianIndex(-1, 0)
    RIGHT = CartesianIndex(1, 0)
    UP = CartesianIndex(0, 1)
    DOWN = CartesianIndex(0, -1)

    actions = [rand([LEFT,RIGHT,UP,DOWN]) for _ in game.snakes]
    is_exit = Ref{Bool}(false)

    on(scene.events.keyboardbuttons) do but
        if ispressed(but, Keyboard.left)
            actions[1] != -LEFT && (actions[1] = LEFT)
        elseif ispressed(but, Keyboard.up)
            actions[1] != -UP && (actions[1] = UP)
        elseif ispressed(but, Keyboard.down)
            actions[1] != -DOWN && (actions[1] = DOWN)
        elseif ispressed(but, Keyboard.right)
            actions[1] != -RIGHT && (actions[1] = RIGHT)
        elseif ispressed(but, Keyboard.a)
            actions[2] != -LEFT && (actions[2] = LEFT)
        elseif ispressed(but, Keyboard.w)
            actions[2] != -UP && (actions[2] = UP)
        elseif ispressed(but, Keyboard.s)
            actions[2] != -DOWN && (actions[2] = DOWN)
        elseif ispressed(but, Keyboard.d)
            actions[2] != -RIGHT && (actions[2] = RIGHT)
        elseif ispressed(but, Keyboard.j)
            actions[3] != -LEFT && (actions[3] = LEFT)
        elseif ispressed(but, Keyboard.i)
            actions[3] != -UP && (actions[3] = UP)
        elseif ispressed(but, Keyboard.k)
            actions[3] != -DOWN && (actions[3] = DOWN)
        elseif ispressed(but, Keyboard.l)
            actions[3] != -RIGHT && (actions[3] = RIGHT)
        elseif ispressed(but, Keyboard.q)
            is_exit[] = true
        end
    end

    record(scene, f_name; framerate=framerate) do io
        while true
            sleep(1)
            is_success = game(actions)
            game_node[] = game
            recordframe!(io)
            is_success || break
            is_exit[] && break
        end
    end
    println("game over")
end


text:
function take_action(::Discord.Client, ::Discord.Message, action::BotAction)
    @error "Unsupported action $action"
end

function take_action(
    client::Discord.Client,
    message::Discord.Message,
    create_trial_post::CreateTrialPost
)
    trial_post = """
    Trial - $(create_trial_post.trialname)
    $(create_trial_post.trialdate) :: $(create_trial_post.trialtime)

    Tanks:

    Healers:

    DPS:

    Or something, idfk
    """

    create_message(client, message.channel_id, content=trial_post)
end

"""
Create a message as a reply to the provided discord message.
"""
function take_action(client::Discord.Client, message::Discord.Message, reply::Reply)
    create_message(client, message.channel_id, content=reply.content)
end


text:
# Testing Glov pretrained embedding with LSTM model
include("../src/ModelPreTrained.jl")
model, test_data = train()
Flux.reset!(model)
test(model,test_data)

max_stars_repo_name: akaysh/Rango

max_issues_repo_name: akaysh/Rango

max_forks_repo_name: akaysh/Rango


text:
module TestTextAnalysis
using SparseArrays
using Test
using Languages
using TextAnalysis
using WordTokenizers

# @testset "TextAnalysis" begin

println("Running tests:")

include("tokenizer.jl")
include("ngramizer.jl")
include("document.jl")
include("metadata.jl")
include("corpus.jl")
include("preprocessing.jl")
include("dtm.jl")
include("stemmer.jl")
include("tf_idf.jl")
include("lda.jl")
include("summarizer.jl")
include("sentiment.jl")
include("bayes.jl")
include("rouge.jl")
include("averagePerceptronTagger.jl")
include("utils.jl")

# end
end


text:
module SnoopCompileBot

@warn """This package is deprecated and renamed to CompileBot. See https://github.com/aminya/CompileBot.jl for more information."""

using Reexport

@reexport using CompileBot

end # module


text:
obj_path = joinpath(@__DIR__, "035_power_drill/textured_simple.obj")
texture_path = joinpath(@__DIR__, "035_power_drill/texture_map.png")

include("test_depth.jl")
include("test_rgb_basic.jl")
include("test_rgb.jl")
include("test_texture.jl")
include("test_texture_mixed.jl")


text:
# Figures 6.3a and 6.3b in "Infinite GMRES for parameterized linear systems"

using SparseArrays
using LinearAlgebra
using PyPlot
using Statistics
pygui(true)
include("infgmres.jl");
include("create_helmholtz_2d.jl");

data=matread("remove_mu.mat")
A0=data["A0"]; A1=data["A1"];
A2=data["A2"]; A3=data["A3"];
A4=data["A3"]; b=vec(data["bvec_2d"]);

# Number of derivatives
N=51;
mu=.1;
kk=10;

# Discretize the problem
myproblem=create_helmholtz2D(N,A0,A1,A2,A3,A4,b);

how_many = 5;
timestamps_mat=zeros(kk,how_many); timestamps=zeros(kk);
timestamps_mat2=zeros(kk,how_many); timestamps2=zeros(kk);
for i=1:how_many
    global ff, hist1, timestamps
    global ff2, hist2, timestamps2
    (ff,hist1)=infgmres(myproblem,kk,method=:tiar);
    (ff2,hist2)=infgmres(myproblem,kk,method=:iar);
    timestamps=hist1[:time_count][1:kk] .+ hist1[:precomp_time]
    timestamps_mat[:,i]=timestamps;
    timestamps2=hist2[:time_count][1:kk] .+ hist2[:precomp_time]
    timestamps_mat2[:,i]=timestamps2;
end
for i=1:kk
    timestamps[i]=float(median(timestamps_mat[i,:]))
    timestamps2[i]=float(median(timestamps_mat2[i,:]))
end

errv=zeros(kk);
errv2=zeros(kk);
for k=1:kk-1
    xx=ff(mu,k)
    xx2=ff2(mu,k)
    errv[k]=norm(myproblem.Afun(mu)*xx-b)/norm(b)
    errv2[k]=norm(myproblem.Afun(mu)*xx2-b)/norm(b)
end
its=range(1,stop=kk,length=kk);

figure(1)
ax = PyPlot.gca()
ax.set_yscale("log")
xlabel("its")
ylabel("res")
semilogy(its,errv,"o")
semilogy(its,errv2,"+")
legend(["Tensor","Alg 1"])

figure(2)
ax = PyPlot.gca()
ax.set_yscale("log")
xlabel("CPU-time")
ylabel("res")
semilogy(timestamps,errv,"o")
semilogy(timestamps2,errv2,"+")
legend(["Tensor","Alg 1"])

max_stars_repo_name: siobhanie/InfGMRES

max_issues_repo_name: siobhanie/InfGMRES

max_forks_repo_name: siobhanie/InfGMRES