text:
export alias_elimination

function flatten(sys::ODESystem)
    if isempty(sys.systems)
        return sys
    else
        return ODESystem(equations(sys),
                         independent_variable(sys),
                         states(sys),
                         parameters(sys),
                         observed=observed(sys))
    end
end

using SymbolicUtils: Rewriters

function fixpoint_sub(x, dict)
    y = substitute(x, dict)
    while !isequal(x, y)
        y = x
        x = substitute(y, dict)
    end

    return x
end

function substitute_aliases(diffeqs, dict)
    lhss(diffeqs) .~ fixpoint_sub.(rhss(diffeqs), (dict,))
end

# Note that we reduce parameters, too
# i.e. `2param = 3` will be reduced away
isvar(s::Sym) = true
isvar(s::Term) = isvar(operation(s))
isvar(s::Any) = false

function get_α_x(αx)
    if isvar(αx)
        return 1, αx
    elseif αx isa Term && operation(αx) === (*)
        args = arguments(αx)
        nums = []
        syms = []
        for arg in args
            isvar(arg) ? push!(syms, arg) : push!(nums, arg)
        end

        if length(syms) == 1
            return prod(nums), syms[1]
        end
    else
        return nothing
    end
end

function alias_elimination(sys::ODESystem)
    eqs = vcat(equations(sys), observed(sys))
    subs = Pair[]
    diff_vars = filter(!isnothing, map(eqs) do eq
            if isdiffeq(eq)
                eq.lhs.args[1]
            else
                nothing
            end
        end) |> Set

    # only substitute when the variable is algebraic
    del = Int[]
    for (i, eq) in enumerate(eqs)
        isdiffeq(eq) && continue
        res_left = get_α_x(eq.lhs)
        if !isnothing(res_left) && !(res_left[2] in diff_vars)
            # `α x = rhs` => `x = rhs / α`
            α, x = res_left
            push!(subs, x => _isone(α) ? eq.rhs : eq.rhs / α)
            push!(del, i)
        else
            res_right = get_α_x(eq.rhs)
            if !isnothing(res_right) && !(res_right[2] in diff_vars)
                # `lhs = β y` => `y = lhs / β`
                β, y = res_right
                push!(subs, y => _isone(β) ? eq.lhs : β * eq.lhs)
                push!(del, i)
            end
        end
    end
    deleteat!(eqs, del)

    eqs′ = substitute_aliases(eqs, Dict(subs))
    alias_vars = first.(subs)

    newstates = setdiff(states(sys), alias_vars)
    ODESystem(eqs′, sys.iv, newstates, parameters(sys), observed=alias_vars .~ last.(subs))
end


text:
# Solving ESPPRC
# Modified from the example given by Google OR-Tools https://developers.google.com/optimization/routing/vrp

using Routing
using Random

# This example is modifed from the VRPTW example from OR-Tools
# But ESPPRC can be used outside the VRP context
# This example describes non-VRP context

# Number of nodes
n_nodes = 17

# cost matrix:
# The objective of ESPPRC is to minimize the total cost
# Randomly generated in this example
# cost allows negative values
cost = - rand(n_nodes, n_nodes)

# travel_time matrix:
# Some elements are Inf, meaning disconnected
# travel_time is used to check the time windows constraints
travel_time = Float64[
    0 8 3 2 6 8 4 8 8 13 7 5 8 12 10 Inf 6;
    8 0 11 10 6 3 9 5 8 4 15 14 13 9 18 9 9;
    3 11 0 1 7 10 6 10 10 14 6 7 9 14 6 16 8;
    2 10 1 0 6 9 4 8 9 13 4 6 8 12 8 14 7;
    6 6 7 6 0 2 3 2 2 7 9 7 7 6 12 8 3;
    8 3 10 9 2 0 6 2 5 4 Inf 10 10 6 Inf 5 6;
    4 9 6 4 3 6 0 4 4 8 5 4 3 7 8 10 2;
    8 5 10 8 2 2 4 0 3 4 9 8 7 3 13 6 3;
    8 8 10 9 2 5 4 3 0 4 6 5 4 3 9 5 2;
    13 4 14 13 7 4 8 4 4 0 10 9 8 4 13 4 6;
    7 Inf 6 4 9 Inf 5 9 6 10 0 1 3 7 3 10 6;
    5 14 7 6 7 10 4 8 5 9 1 0 2 6 4 8 4;
    8 13 9 8 7 10 3 7 4 8 3 2 0 4 5 6 4;
    12 9 14 12 6 6 7 3 3 4 7 6 4 0 9 2 5;
    10 18 6 8 12 15 8 13 9 13 3 4 5 9 0 9 9;
    Inf 9 16 14 8 5 10 6 5 4 10 8 6 2 9 0 7;
    6 9 8 7 3 6 2 3 2 6 6 4 4 5 9 7 0
]

# load matrix
# The load information is for each arc, hence given as a matrix
# If you have load on each node, then it should be converted to a matrix
# For example, load[i, j] = load_on_node[j]
load_on_node = [1, 1, 2, 4, 2, 4, 8, 8, 1, 2, 1, 2, 4, 4, 8, 8, 0]
load = zeros(n_nodes, n_nodes)
for i in 1:n_nodes, j in 1:n_nodes
    load[i, j] = load_on_node[j]
end

# Time windows, for all nodes, including origin and destination
time_windows = [
    (0, 12),    #  1
    (0, 15),    #  2
    (16, 28),   #  3
    (10, 13),   #  4
    (0, 5),     #  5
    (5, 10),    #  6
    (0, 4),     #  7
    (5, 10),    #  8
    (0, 3),     #  9
    (10, 16),   # 10
    (10, 15),   # 11
    (0, 5),     # 12
    (5, 10),    # 13
    (7, 8),     # 14
    (10, 15),   # 15
    (11, 15),   # 16
    (0, 35),    # 17
]
# If early_time is set for the origin, it means the earliest possible departure time.
early_time = [time_windows[i][1] for i in 1:n_nodes]
late_time = [time_windows[i][2] for i in 1:n_nodes]

# vehicle capacity
capacity = 15

# service_time is just set to all zeros in this example
service_time = zeros(n_nodes)

# origin and destination
origin, destination = 6, 3

pg = ESPPRC_Instance(
    origin,
    destination,
    capacity,
    cost,
    travel_time,
    load,
    early_time,
    late_time,
    service_time
)

# Solve the ESPPRC
solution_label = solveESPPRC(pg)

@show solution_label.path
@show solution_label.cost
@show solution_label.time
@show solution_label.load

ext: jl


text:
# ---
# title: 593. Valid Square
# id: problem593
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Math
# link: <https://leetcode.com/problems/valid-square/description/>
# hidden: true
# ---
#
# Given the coordinates of four points in 2D space `p1`, `p2`, `p3` and `p4`,
# return `true` _if  the four points construct a square_.
#
# The coordinate of a point `pi` is represented as `[xi, yi]`. The input is
# **not** given in any order.
#
# A **valid square** has four equal sides with positive length and four equal
# angles (90-degree angles).
#
#
#
# **Example 1:**
#
#
#
#     Input: p1 = [0,0], p2 = [1,1], p3 = [1,0], p4 = [0,1]
#     Output: true
#
#
# **Example 2:**
#
#
#
#     Input: p1 = [0,0], p2 = [1,1], p3 = [1,0], p4 = [0,12]
#     Output: false
#
#
# **Example 3:**
#
#
#
#     Input: p1 = [1,0], p2 = [-1,0], p3 = [0,1], p4 = [0,-1]
#     Output: true
#
#
#
#
# **Constraints:**
#
#   * `p1.length == p2.length == p3.length == p4.length == 2`
#   * `-104 <= xi, yi <= 104`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
# X "mime" types

# These are possible return types from the X11 clipboard which are not
# actual MIME types (for the most part), but which are text and can be
# treated as such.

for mime in ["STRING", "UTF8_STRING", "text/plain;charset=utf-8", "text/plain;charset=UTF-8", "TEXT", "COMPOUND_TEXT"]
    @eval @textmime $mime
    @eval decodemime(m::MIME{symbol($mime)}, x) = bytestring(x)
end

for mime in ["text/vnd.graphviz", "text/latex", "text/calendar", "text/n3", "text/richtext", "text/x-setext", "text/sgml", "text/tab-separated-values", "text/x-vcalendar", "text/x-vcard", "text/cmd", "text/css", "text/csv", "text/html", "text/javascript", "text/plain", "text/vcard", "text/xml", "application/atom+xml", "application/ecmascript", "application/json", "application/rdf+xml", "application/rss+xml", "application/xml-dtd", "application/postscript", "image/svg+xml", "application/x-latex", "application/xhtml+xml", "application/javascript", "application/xml", "model/x3d+xml", "model/x3d+vrml", "model/vrml"]
    @eval decodemime(m::MIME{symbol($mime)}, x) = bytestring(x)
end

decodemime(m::MIME, x) = x


text:
function calc_nu_J(reactor, R_0, B_0)
  tmp_I_M = calc_I_M(reactor, R_0, B_0)

  q_star = (
    5 * reactor.epsilon ^ 2 * reactor.kappa ^ 1.27 *
    B_0 * R_0 / tmp_I_M
  )

  1 - sqrt( q_star / 4 )
end


text:
module Fields

export Face, Center
export AbstractField, AbstractDataField, Field
export CenterField, XFaceField, YFaceField, ZFaceField
export ReducedField, AveragedField, ComputedField, KernelComputedField, BackgroundField
export interior, data
export xnode, ynode, znode, location
export set!, compute!, @compute
export VelocityFields, TracerFields, tracernames, PressureFields, TendencyFields
export interpolate, FieldSlicer

using Oceananigans.Architectures
using Oceananigans.Grids
using Oceananigans.BoundaryConditions

include("abstract_field.jl")
include("reduced_getindex_setindex.jl")
include("field.jl")
include("zero_field.jl")
include("reduced_field.jl")
include("averaged_field.jl")
include("computed_field.jl")
include("kernel_computed_field.jl")
include("pressure_field.jl")
include("function_field.jl")
include("set!.jl")
include("tracer_names.jl")
include("validate_field_tuple_grid.jl")
include("field_tuples.jl")
include("background_fields.jl")
include("interpolate.jl")
include("field_slicer.jl")
include("show_fields.jl")
include("broadcasting_abstract_fields.jl")
include("mapreduce_abstract_fields.jl")

# Fallback: cannot infer boundary conditions.
boundary_conditions(field) = nothing
boundary_conditions(f::Union{Field, ReducedField, ComputedField, KernelComputedField}) = f.boundary_conditions

end


text:
# load a dataset
using RDatasets
iris = dataset("datasets", "iris");
println(iris)
# load the StatPlots recipes (for DataFrames) available via:
# Pkg.add("StatPlots")
using StatPlots

# Scatter plot with some custom settings
@df iris scatter(:SepalLength, :SepalWidth, group=:Species,
        title = "Iris awesome plot from dataset",
        xlabel = "Length", ylabel = "Width",
        m=(0.5, [:cross :hex :star7], 12),
        bg=RGB(.2,.2,.2))

# save a png
png("iris")

max_stars_repo_name: MarsBighead/mustang

max_issues_repo_name: MarsBighead/mustang

max_forks_repo_name: MarsBighead/mustang


text:
@inline function two_sum_tail(a::T, b::T, x::T) where {T}
    v = x - a
    return (a - (x - v)) + (b - v)
end

@inline function fast_two_sum_tail(a::T, b::T, x::T) where {T}
    v = x - a
    return b - v
end

@inline function two_sum(a::T, b::T) where {T}
    x = a + b
    return two_sum_tail(a, b, x), x
end

@inline function fast_two_sum(a::T, b::T) where {T}
    x = a + b
    return fast_two_sum_tail(a, b, x), x
end

@inline function two_diff_tail(a::T, b::T, x::T) where {T}
    b1 = a - x
    a1 = x + b1
    return (a - a1) + (b1 - b)
end

@inline function two_diff(a::T, b::T) where {T}
    x = a - b
    return two_diff_tail(a, b, x), x
end

@inline two_prod_tail(a::T, b::T, x::T) where {T} = fma(a, b, -x)

@inline function two_prod(a::T, b::T) where {T}
    x = a * b
    return two_prod_tail(a, b, x), x
end

@inline grow_expansion(a::T, b::T) where {T} = two_sum(a, b)
@inline grow_expansion(a) = a

@generated function grow_expansion(e::NTuple{N,T}, b::T) where {N,T}
    h = [Symbol(:h, i) for i in 1:N]
    body = [:(($(h[i]), b) = two_sum(b, e[$i])) for i in 1:N]
    return quote
        $(body...)
        return ($(h...), b)
    end
end

@generated function expansion_sum(e::NTuple{N,T}, f::NTuple{M,T}) where {N,M,T}
    h = [Symbol(:h, i) for i in 1:(N+M)]
    body = [
        :(($(h[i:i+N]...),) = grow_expansion(($(h[i:i+N-1]...),),f[$i])) for i in 1:M
    ]
    return quote
        ($(h[1:N]...),) = e
        $(body...)
        return ($(h...),)
    end
end

@inline function expansion_diff(e::NTuple{2,T}, b::T) where {T}
    x1, x2 = two_diff(e[1], b)
    x2, x3 = two_sum(e[2], x2)
    return x1, x2, x3
end

@inline function expansion_diff(e::NTuple{2,T}, f::NTuple{2,T}) where {T}
    x1, x2, x3 = expansion_diff((e[2], e[1]), f[1])
    x2, x3, x4 = expansion_diff((x3, x2), f[2])
    return x1, x2, x3, x4
end

@generated function expansion_scale(e::NTuple{N,T}, b::T) where {N,T}
    h = [Symbol(:h, i) for i in 1:(2N-1)]
    Q = [Symbol(:Q, i) for i in 2:2N]
    t = [Symbol(:t, i) for i in 2:N]
    s = [Symbol(:s, i) for i in 2:N]
    body = [
        quote
            ($(t[i-1]), $(s[i-1])) = two_prod(e[$i], b)
            ($(h[2i-2]), $(Q[2i-2])) = two_sum($(t[i-1]), $(Q[2i-3]))
            ($(h[2i-1]), $(Q[2i-1])) = fast_two_sum($(s[i-1]), $(Q[2i-2]))
        end for i in 2:N
    ]
    return quote
        (h1, Q2) = two_prod(e[1], b)
        $(body...)
        return ($(h...), $(Q[end]))
    end
end

approximate(e) = sum(e)

max_stars_repo_name: PaulMatlashewski/DelaunayTriangulations

max_issues_repo_name: PaulMatlashewski/DelaunayTriangulations

max_forks_repo_name: PaulMatlashewski/DelaunayTriangulations


text:
using SMLMData

"""
    clusterdata = organizeclusters(smld::SMLMData.SMLD2D)

Organize pre-clusters into a vector indexing distinct clusters.

# Description
Pre-clusters in the input `smld`, as related by their shared integer value k_off
`smld.connectID`, are organized into a vector of matrices `clusterdata`.  Each
index of `clusterdata` corresponds to a distinct cluster.
"""
function organizeclusters(smld::SMLMData.SMLD2D)
    # Isolate and sort some arrays from the smd structure.
    connectID = smld.connectID
    sortindices = sortperm(connectID)
    combineddata = [smld.x smld.y smld.σ_x smld.σ_y smld.framenum smld.datasetnum smld.connectID]
    combineddata = combineddata[sortindices, :]

    # Loop over preclusters and organize member localizations into a more
    # accessible format.
    npercluster = counts(connectID)
    npercluster = npercluster[npercluster.!=0]
    ncumulative = [0; cumsum(npercluster)]
    uniqueIDs = unique(connectID)
    nclusters = length(uniqueIDs)
    clusterdata = Vector{Matrix{Float32}}(undef, nclusters)
    for nn = 1:nclusters
        clusterind = (1:npercluster[nn]) .+ ncumulative[nn]
        clusterdata[nn] = [combineddata[clusterind, :] sortindices[clusterind]]
    end

    return clusterdata
end


text:
using Documenter, SIIPExamples

makedocs(
    modules = [SIIPExamples],
    format = Documenter.HTML(),
    sitename = "SIIPExamples.jl",
    authors = "Clayton Barrows",
    pages = [
        "Home" => "index.md",
        "Adding a New Tutorial" => "new.md",
        "Notes" => "notes.md",
        "Function Index" => "api.md",
    ],
)

Documenter.deploydocs(repo = "github.com/NREL-SIIP/SIIPExamples.jl.git", target = "build")


text:
# first run Bolivia_plotinterpolatedPTR3data.jl

################ now do same comparison as above but for single Cnr, now divide in barchart by H-nr! ##############################################
singleCnrs = collect(1:20)
for singleCnr in singleCnrs
	global singleCnrMask = findall(x -> x == singleCnr, composition[:,1])
	global singleCnrtraces = ptr3Traces[:,singleCnrMask]
	Hnrs = MasslistFunctions.masslistPos(round(singleCnr/2)*2-4):1:2*singleCnr+3
	global Hnrtraces = DatasetFunctions.nansum(singleCnrtraces[:,findall(x -> x == Hnrs[1], composition[singleCnrMask,3])],2)
	for n in Hnrs[2:end]
		global Hnrtraces = hcat(Hnrtraces, DatasetFunctions.nansum(singleCnrtraces[:,findall(x -> x == n, composition[singleCnrMask,3])],2))
	end

	n = 1
	global meanHnrContrRel = DatasetFunctions.nanmean((Hnrtraces./DatasetFunctions.nansum(singleCnrtraces,2))[findall(x -> avTimes[n,1] <= x <= avTimes[n,2], ptr3unixTimes),:],1)
	for n in 2:size(avTimes)[1]
		global meanHnrContrRel = vcat(meanHnrContrRel, DatasetFunctions.nanmean((Hnrtraces./DatasetFunctions.nansum(singleCnrtraces,2))[findall(x -> avTimes[n,1] <= x <= avTimes[n,2], ptr3unixTimes),:],1))
	end
	n = 1
	global meanHnrContrAbs = DatasetFunctions.nanmean((Hnrtraces)[findall(x -> avTimes[n,1] <= x <= avTimes[n,2], ptr3unixTimes),:],1)
	for n in 2:size(avTimes)[1]
		global meanHnrContrAbs = vcat(meanHnrContrAbs, DatasetFunctions.nanmean((Hnrtraces)[findall(x -> avTimes[n,1] <= x <= avTimes[n,2], ptr3unixTimes),:],1))
	end

	# absolute
	figure()
	global x = collect(1:size(Hnrs)[1])
	width = 0.35
	bar(x .- width*0.5, (meanHnrContrAbs[1,:] .+ meanHnrContrAbs[2,:])./2, width, label = "mean of $(citydays[1]). and $(citydays[2]).")
	bar(x .+ width*0.5, (meanHnrContrAbs[3,:] .+ meanHnrContrAbs[4,:])./2, width, label = "mean of $(ruraldays[1]). and $(ruraldays[2]).")
	yscale("log")
	ylabel("conc. [ppt]")
	xlabel("#H")
	xticks(x, Int32.(Hnrs))
	legend()
	#text(0.2, 0.2, L"$[\overline{C_xH_y}]_{city}-[\overline{C_xH_y}]_{rural}$")	#, xy=[0.2;0.2], xycoords="axes fraction", fontsize=14.0, ha="right", va="bottom")
	if excludeCHON
		title("summed C$(singleCnr) peaks by Hnr (only CHO)")
		savefig("$(fp)_HnrShifts_rural_city_Cnr$(singleCnr)_abs_$(firstday)-$(lastday)May_CHOonly.png")
	elseif excludeCHO
		title("summed C$(singleCnr) peaks by Hnr (only CHON)")
		savefig("$(fp)_HnrShifts_rural_city_Cnr$(singleCnr)_abs_$(firstday)-$(lastday)May_CHONonly.png")
	elseif showboth
		title("summed C$(singleCnr) peaks by Hnr (all)")
		savefig("$(fp)_HnrShifts_rural_city_Cnr$(singleCnr)_abs_$(firstday)-$(lastday)May.png")
	end

	#relative
	figure()
	global x = collect(1:size(Hnrs)[1])
	width = 0.7
	bar(x, 100 .* ((meanHnrContrRel[1,:] .+ meanHnrContrRel[2,:])./2 .- (meanHnrContrRel[3,:] .+ meanHnrContrRel[4,:])./2), width, label = "\"city\" - \"rural\", $(firstday).-$(lastday)")
	yscale("linear")
	ylabel("change of rel. contribution [%]")
	title("relative contribution of summed peaks by Hnr of C$(singleCnr) compounds")
	xticks(x, Int32.(Hnrs))
	xlabel("#H")
	axhline(0)
	annotate(L"$\frac{[\overline{C_xH_y}]_{city}-[\overline{C_xH_y}]_{rural}}{[\overline{\sum_{y} C_xH_y}]_{tot}}$", xy=[0.4;0.2], xycoords="axes fraction", fontsize=14.0, ha="right", va="bottom")
	legend()
	if excludeCHON
		title("relative contribution of summed peaks by Hnr of C$(singleCnr) compounds (only CHO)")
		savefig("$(fp)_HnrShifts_rural_city_Cnr$(singleCnr)_rel_$(firstday)-$(lastday)May_CHOonly.png")
	elseif excludeCHO
		title("relative contribution of summed peaks by Hnr of C$(singleCnr) compounds (only CHON)")
		savefig("$(fp)_HnrShifts_rural_city_Cnr$(singleCnr)_rel_$(firstday)-$(lastday)May_CHONonly.png")
	elseif showboth
		title("relative contribution of summed peaks by Hnr of C$(singleCnr) compounds (all)")
		savefig("$(fp)_HnrShifts_rural_city_Cnr$(singleCnr)_rel_$(firstday)-$(lastday)May.png")
	end
end


text:
include("conway_test.jl")
using .ConwayTest #Local modules have to be preceded with a dot
turbine_test()
blinker_test()


text:
using ABC
using Distributions
using Base.Test

tic()
srand(1234)
include("test1.jl")
@time @test test1()
#=
err_code = test1()
if err_code !=0
  println("# Test 1 FAILED with error code ", err_code)
else
  println("# Test 1 passed.")
end
=#
toc()


text:
module DocMetaTests
using Test
using Documenter

module TestMod
module Submodule end
end

@testset "DocMeta" begin
    @test DocMeta.getdocmeta(TestMod) == Dict()
    @test DocMeta.getdocmeta(TestMod, :DocTestSetup) === nothing
    @test DocMeta.getdocmeta(TestMod, :DocTestSetup, 42) === 42
    @test DocMeta.setdocmeta!(TestMod, :DocTestSetup, :foo) === nothing
    @test DocMeta.getdocmeta(TestMod) == Dict(:DocTestSetup => :foo)
    @test DocMeta.getdocmeta(TestMod, :DocTestSetup) == :foo
    @test DocMeta.getdocmeta(TestMod, :DocTestSetup, 42) == :foo
    # bad key
    @test_throws ArgumentError DocMeta.setdocmeta!(TestMod, :FooBar, 0)
    # bad argument type
    @test_throws ArgumentError DocMeta.setdocmeta!(TestMod, :DocTestSetup, 42)
    # setting again works
    @test DocMeta.setdocmeta!(TestMod, :DocTestSetup, :foo; warn=false) === nothing
    # recursive setting
    @test DocMeta.getdocmeta(TestMod, :DocTestSetup) == :foo
    @test DocMeta.getdocmeta(TestMod.Submodule, :DocTestSetup) === nothing
    @test DocMeta.setdocmeta!(TestMod, :DocTestSetup, :foo; recursive=true, warn=false) === nothing
    @test DocMeta.getdocmeta(TestMod, :DocTestSetup) == :foo
    @test DocMeta.getdocmeta(TestMod.Submodule, :DocTestSetup) == :foo
end

end


text:
"""
    uniform(dims...)

Return an `Array` of size `dims` containing random variables taken from a uniform
distribution in the interval ``[-1, 1]``.

# Examples
```jldoctest; setup = :(using Random; Random.seed!(0))
julia> uniform(2, 3)
2×3 Array{Float32,2}:
 0.548721  -0.524116  -0.743921
 0.822376   0.735768   0.0525029
```
"""
uniform(dims...) = (rand(dims...) .- 0.5f0) * 2

"""
    SIREN_init(SIREN_init(omega_0::Real, is_first::Bool, dims::Integer...)

Return an `Array` of size `dims` containing random variables.
The distribution is supposed to be an implementation of the SIREN initialization
from the original paper (Sec 3.2, Appdx 1, Appdx 1.5). My explanation will probably
be wrong so check out their explanation.

# Examples
```jldoctest; setup = :(using Random; Random.seed!(0))
julia> SIREN_init(30,false,3,2)
3×2 Array{Float64,2}:
  0.0305138  -0.0304217
  0.0386888  -0.0208474
 -0.031625   -0.0279565
```
"""
function SIREN_init(omega_0::Real, is_first::Bool, dims::Integer...)
    w = uniform(dims...)

    fan_in = 1
    s = size(w)
    for i in 1:(length(s)-1)
      fan_in *= s[i]
    end

    w .*= sqrt(6/fan_in)

    if !is_first
      w ./= omega_0
    end

    return w
end

max_issues_repo_name: MacKenzieHnC/SIREN

max_forks_repo_name: MacKenzieHnC/SIREN


text:
module train
using ..gpflow
import ..gpflow: instantiate!, minimize!
export
    ScipyOptimizer,
    instantiate!,
    minimize!

mutable struct ScipyOptimizer<:Optimizer
    o::Union{PyObject, Nothing}
end

function ScipyOptimizer()
    out = ScipyOptimizer(nothing)
    instantiate!(out)
    out
end

function minimize!(opt::Union{ScipyOptimizer, Nothing}, m::Model)
opt.o.minimize(m.o)
end

function instantiate!(o::Union{ScipyOptimizer, Nothing})
    if o === nothing return nothing end
    if typeof(o.o)<:PyObject return o.o end
    o.o = py_gpflow.train.ScipyOptimizer()
end

end # module


text:
# in this overload we will work on some meta/debug information management
function CC.typeinf(interp::JETInterpreter, frame::InferenceState)
    # # print debug info before typeinf
    # depth = interp.depth[]
    # io = stdout::IO
    # color = RAIL_COLORS[(depth+1)%N_RAILS+1]
    # print_rails(io, depth)
    # printstyled(io, "┌ @ "; color)
    # print(io, frame.linfo)
    # if is_constant_propagated(frame)
    #     printstyled(io, " (constant prop': ", frame.result.argtypes, ')'; color = :cyan)
    # end
    # file, line = get_file_line(frame.linfo)
    # print(io, ' ', file, ':', line)
    # println(io)

    interp.depth[] += 1 # for debug

    ret = @invoke typeinf(interp::AbstractInterpreter, frame::InferenceState)

    push!(ANALYZED_LINFOS, frame.linfo) # analyzed !

    interp.depth[] -= 1 # for debug

    # # print debug info after typeinf
    # print_rails(io, depth)
    # printstyled(io, "└─→ "; color)
    # result = get_result(frame)
    # isa(result, InferenceState) || printstyled(io, get_result(frame); color = TYPE_ANNOTATION_COLOR)
    # printlnstyled(io, " ($(length(interp.reports)) reports)"; color = ERROR_COLOR)

    return ret
end

# in this overload we can work on `frame.src::CodeInfo` (and also `frame::InferenceState`)
# where type inference (and also optimization if applied) already ran on
function CC._typeinf(interp::JETInterpreter, frame::InferenceState)
    linfo = frame.linfo

    # some methods like `getproperty` can't propagate accurate types without actual values,
    # and constant prop' plays a somewhat critical role in those cases by overwriteing the
    # previously-inferred lousy result; JET.jl also needs that to reduce false positive reports,
    # and so here we will throw-away previously-collected error reports that are "lineage" of
    # this frame, when it is being re-inferred with constant-prop'ed inputs
    #
    # NOTE:
    # - constant prop' only happens after inference with non-constant abstract values (i.e. types)
    # - xref to track the change in the native constant propagation logic:
    #   https://github.com/JuliaLang/julia/blob/a108d6cb8fdc7924fe2b8d831251142386cb6525/base/compiler/abstractinterpretation.jl#L153
    # IDEA:
    # we may want to keep some "serious" error reports like `GlobalUndefVarErrorReport`
    # even when constant prop' reveals it never happens given the current constant arguments
    if is_constant_propagated(frame)
        # check from `frame.parent` is because there might be other analysis for this
        # `frame` with different constants in this abstract interpretation pass
        filter!(r->!is_lineage(r.lineage, frame.parent, linfo), interp.reports)
    end

    before = Set(interp.reports)

    ret = @invoke _typeinf(interp::AbstractInterpreter, frame::InferenceState)

    stmts = frame.src.code

    # report (local) undef var error
    # this only works when optimization is enabled, just because `:throw_undef_if_not` and
    # `:(unreachable)` are introduced by `optimize`
    if may_optimize(interp)
        for (idx, stmt) in enumerate(stmts)
            if isa(stmt, Expr) && stmt.head === :throw_undef_if_not
                sym, _ = stmt.args
                next_idx = idx + 1
                if checkbounds(Bool, stmts, next_idx) && is_unreachable(@inbounds stmts[next_idx])
                    # the optimization so far has found this statement is never "reachable";
                    # JET reports it since it will invoke undef var error at runtime, or will just
                    # be dead code otherwise

                    report!(interp, LocalUndefVarErrorReport(interp, frame, sym))
                # else
                    # by excluding this pass, JET accepts some false negatives (i.e. don't report
                    # those that may actually happen on actual execution)
                end
            end
        end
    end

    after = Set(interp.reports)
    reports_for_this_linfo = setdiff(after, before)

    # report `throw` calls "appropriately"
    if get_result(frame) === Bottom
        # if the return type here is `Bottom` annotated, this _may_ mean there're uncaught
        # `throw` calls
        # XXX: well, it's possible that the `throw` calls within them are all caught but the
        # other critical errors make the return type `Bottom`
        # NOTE: to reduce the false positive `UncaughtExceptionReport`s described above, we count
        # `throw` calls here after optimization, since it may have eliminated "unreachable"
        # `throw` calls
        codelocs    = frame.src.codelocs
        linetable   = frame.src.linetable
        throw_locs  = LineInfoNode[]
        throw_calls = Any[]
        for r in reports_for_this_linfo
            if isa(r, ExceptionReport)
                push!(throw_locs, r.lin)
            end
        end
        for (i, stmt) in enumerate(stmts)
            is_throw_call′(stmt) || continue
            # if this `throw` is already reported, don't duplciate
            linetable[codelocs[i]] in throw_locs && continue
            push!(throw_calls, stmt)
        end
        isempty(throw_calls) || push!(interp.uncaught_exceptions, UncaughtExceptionReport(interp, frame, throw_calls))
    else
        # the non-`Bottom` result here means `throw` calls within the children frames
        # reported so far (if exist) are caught and not propagated to the result;
        # we don't want to cache for this frame and its parents, so just filter them away
        # NOTE: this is critical for performance issue https://github.com/aviatesk/JET.jl/issues/71
        isempty(interp.uncaught_exceptions) || empty!(interp.uncaught_exceptions)
    end

    isempty(interp.uncaught_exceptions) || push!(reports_for_this_linfo, interp.uncaught_exceptions...)

    if !isempty(reports_for_this_linfo)
        if is_constant_propagated(frame)
            argtypes = frame.result.argtypes
            cache = interp.cache

            # if haskey(cache, argtypes)
            #     @warn "this control flow shouldn't happen ..." linfo argtypes
            # end
            local_cache = cache[argtypes] = InferenceErrorReportCache[]

            for report in reports_for_this_linfo
                cache_report!(report, linfo, local_cache)
            end
        else
            # if haskey(JET_GLOBAL_CACHE, linfo)
            #     # can happen when the inference for this frame was manually invoked,
            #     # e.g. calling `@report_call sum("julia")`  two times from REPL
            #     # but we will just update the cache anyway
            #     @warn "deplicated analysis happened" linfo
            # end
            global_cache = JET_GLOBAL_CACHE[linfo] = InferenceErrorReportCache[]

            for report in reports_for_this_linfo
                cache_report!(report, linfo, global_cache)
            end
        end
    end

    return ret
end

is_unreachable(@nospecialize(_)) = false
is_unreachable(rn::ReturnNode)   = !isdefined(rn, :val)

is_throw_call′(@nospecialize(_)) = false
is_throw_call′(e::Expr)          = is_throw_call(e)

"""
    function overload_typeinf_edge!()
        ...
    end
    push_inithook!(overload_typeinf_edge!)

the aims of this overload are:
1. invalidate code cache for a `MethodInstance` that is not yet analyzed by JET;
   it happens when the `MethodInstance` is cached within the system image itself
2. when cache is hit (i.e. any further inference won't occcur for the cached frame),
   append cached reports associated with the cached `MethodInstance`
"""
function overload_typeinf_edge!()

# %% for easier interactive update of typeinf_edge
Core.eval(CC, quote

# compute (and cache) an inferred AST and return the current best estimate of the result type
function typeinf_edge(interp::$(JETInterpreter), method::Method, @nospecialize(atypes), sparams::SimpleVector, caller::InferenceState)
    mi = specialize_method(method, atypes, sparams)::MethodInstance
    #=== typeinf_edge patch-point 1 start ===#
    code = $(∉)(mi, $(ANALYZED_LINFOS)) ? nothing : get(code_cache(interp), mi, nothing)
    #=== typeinf_edge patch-point 1 end ===#
    if code isa CodeInstance # return existing rettype if the code is already inferred
        #=== typeinf_edge patch-point 2 start ===#
        # cache hit, now we need to append cached reports associated with this `MethodInstance`
        global_cache = $(get)($(JET_GLOBAL_CACHE), mi, nothing)
        if isa(global_cache, $(Vector{InferenceErrorReportCache}))
            $(foreach)(global_cache) do cached
                $(restore_cached_report!)(cached, interp, caller)
            end
        end
        #=== typeinf_edge patch-point 2 end ===#
        update_valid_age!(caller, WorldRange(min_world(code), max_world(code)))
        if isdefined(code, :rettype_const)
            if isa(code.rettype_const, Vector{Any}) && !(Vector{Any} <: code.rettype)
                return PartialStruct(code.rettype, code.rettype_const), mi
            else
                return Const(code.rettype_const), mi
            end
        else
            return code.rettype, mi
        end
    end
    if ccall(:jl_get_module_infer, Cint, (Any,), method.module) == 0
        return Any, nothing
    end
    if !caller.cached && caller.parent === nothing
        # this caller exists to return to the user
        # (if we asked resolve_call_cyle, it might instead detect that there is a cycle that it can't merge)
        frame = false
    else
        frame = resolve_call_cycle!(interp, mi, caller)
    end
    if frame === false
        # completely new
        lock_mi_inference(interp, mi)
        result = InferenceResult(mi)
        frame = InferenceState(result, #=cached=#true, interp) # always use the cache for edge targets
        if frame === nothing
            # can't get the source for this, so we know nothing
            unlock_mi_inference(interp, mi)
            return Any, nothing
        end
        if caller.cached || caller.limited # don't involve uncached functions in cycle resolution
            frame.parent = caller
        end
        typeinf(interp, frame)
        update_valid_age!(frame, caller)
        return widenconst_bestguess(frame.bestguess), frame.inferred ? mi : nothing
    elseif frame === true
        # unresolvable cycle
        return Any, nothing
    end
    # return the current knowledge about this cycle
    frame = frame::InferenceState
    update_valid_age!(frame, caller)
    return widenconst_bestguess(frame.bestguess), nothing
end

end) # Core.eval(CC, quote
# %% for easier interactive update of typeinf_edge

end # function overload_typeinf_edge!()
push_inithook!(overload_typeinf_edge!)

const ANALYZED_LINFOS  = IdSet{MethodInstance}() # keeps `MethodInstance`s analyzed by JET

const JET_GLOBAL_CACHE = IdDict{MethodInstance,Vector{InferenceErrorReportCache}}()


text:
using FizzBuzzParametricDispatch
using Test

@testset "FizzBuzzParametricDispatch.jl" begin
    @test fizzbuzz.(1:15) == [
        "1",
        "2",
        "Fizz",
        "4",
        "Buzz",
        "Fizz",
        "7",
        "8",
        "Fizz",
        "Buzz",
        "11",
        "Fizz",
        "13",
        "14",
        "FizzBuzz"
    ]
end


text:
# This file is a part of SimilaritySearch.jl
export KnnResult
export maxlength, maxlength, getdist, getid, idview, distview, reuse!

"""
    KnnResult(ksearch::Integer)

Creates a priority queue with fixed capacity (`ksearch`) representing a knn result set.
It starts with zero items and grows with [`push!(res, id, dist)`](@ref) calls until `ksearch`
size is reached. After this only the smallest items based on distance are preserved.
"""
struct KnnResult # <: AbstractVector{Tuple{IdType,DistType}}
    id::Vector{Int32}
    dist::Vector{Float32}
    k::Int  # number of neighbors
end

function KnnResult(k::Integer)
    @assert k > 0
    res = KnnResult(Vector{Int32}(undef, 0), Vector{Float32}(undef, 0), k)
    sizehint!(res.id, k)
    sizehint!(res.dist, k)
    res
end

"""
    _shifted_fixorder!(res::KnnResult, sp, ep)

Sorts the result in place; the possible element out of order is on the last entry always.
It implements a kind of insertion sort that it is efficient due to the expected
distribution of the items being inserted (it is expected just a few elements smaller than the current ones)
"""
function _shifted_fixorder!(res::KnnResult, sp::Int, ep::Int)
    id, dist = res.id, res.dist
    @inbounds i, d = id[ep], dist[ep]
    pos = _find_inspos(dist, sp, ep, d)
    _shift_vector(id, pos, ep, i)
    _shift_vector(dist, pos, ep, d)

    nothing
end

@inline function _find_inspos(dist::Vector{Float32}, sp::Int, ep::Int, d::Float32)
    @inbounds while ep > sp
        ep -= 1
        d < dist[ep] || return ep + 1
    end

    ep
end

@inline function _shift_vector(arr::Vector, sp::Int, ep::Int, val)
    @inbounds while ep > sp
        arr[ep] = arr[ep-1]
        ep -= 1
    end

    @inbounds arr[ep] = val
end

"""
    push!(res::KnnResult, item::Pair)
    push!(res::KnnResult, id::Integer, dist::Real)

Appends an item into the result set
"""
@inline function Base.push!(res::KnnResult, id::Integer, dist::Real; sp=1, k=maxlength(res))
    len = length(res)

    if len < k
        push!(res.id, id)
        push!(res.dist, dist)

        _shifted_fixorder!(res, sp, len+1)
        return true
    end

    dist >= last(res.dist) && return false

    @inbounds res.id[end], res.dist[end] = id, dist
    _shifted_fixorder!(res, sp, len)
    true
end

#@inline Base.push!(res::KnnResult, id::Integer, dist::Real) = push!(res, convert(Int32, id), convert(Float32, dist))
@inline Base.push!(res::KnnResult, p::Pair) = push!(res, p.first, p.second)

"""
    popfirst!(p::KnnResult)

Removes and returns the nearest neeighboor pair from the pool, an O(length(p.pool)) operation
"""
@inline function Base.popfirst!(res::KnnResult)
    popfirst!(res.id) => popfirst!(res.dist)
end

"""
    pop!(res::KnnResult)

Removes and returns the last item in the pool, it is an O(1) operation
"""
@inline function Base.pop!(res::KnnResult)
    pop!(res.id) => pop!(res.dist)
end

"""
    maxlength(res::KnnResult)

The maximum allowed cardinality (the k of knn)
"""
@inline maxlength(res::KnnResult) = res.k
@inline Base.length(res::KnnResult) = length(res.id)

"""
    reuse!(res::KnnResult)
    reuse!(res::KnnResult, k::Integer)

Returns a result set and a new initial state; reuse the memory buffers
"""
@inline function reuse!(res::KnnResult, k::Integer=res.k)
    @assert k > 0
    empty!(res.id); empty!(res.dist)
    # sizehint!(res.id, k), sizehint!(res.dist, k)
    KnnResult(res.id, res.dist, k)
end

"""
    getindex(res::KnnResult, i)

Access the i-th item in `res`
"""
@inline function Base.getindex(res::KnnResult, i)
    @inbounds res.id[i] => res.dist[i]
end

@inline getid(res::KnnResult, i) = @inbounds res.id[i]
@inline getdist(res::KnnResult, i) = @inbounds res.dist[i]

@inline Base.last(res::KnnResult) = last(res.id) => last(res.dist)
@inline Base.first(res::KnnResult) = @inbounds first(res.id) => first(res.dist)
@inline Base.maximum(res::KnnResult) = last(res.dist)
@inline Base.minimum(res::KnnResult) = @inbounds res.dist[1]
@inline Base.argmax(res::KnnResult) = last(res.id)
@inline Base.argmin(res::KnnResult) = @inbounds res.id[1]
@inline Base.firstindex(res::KnnResult) = 1
@inline Base.lastindex(res::KnnResult) = length(res)

@inline idview(res::KnnResult) = res.id
@inline distview(res::KnnResult) = res.dist

@inline Base.eachindex(res::KnnResult) = firstindex(res):lastindex(res)
Base.eltype(res::KnnResult) = Pair{Int32,Float32}

##### iterator interface
### KnnResult
"""
    Base.iterate(res::KnnResult, state::Int=1)

Support for iteration
"""
function Base.iterate(res::KnnResult, i::Int=1)
    n = length(res)
    (n == 0 || i > n) && return nothing
    @inbounds res[i], i+1
end

max_stars_repo_name: sadit/SimilaritySearch

max_issues_repo_name: sadit/SimilaritySearch

max_forks_repo_name: sadit/SimilaritySearch


text:
#Written By Johnathan Bizzano
export @nullc, eqerror, getbit, setbit, bytesof,
matrixbytesof, objectof, matrixobjectof, issupertype, matrixbitsof, bitsof,
arraybytesof, arrayobjectof, arraybitsof, pass_func, proto, findOrIsMethod, FakeMethod, abstract

"Puts abstract implementation error"
macro abstract(expr)
      return :($expr = error("Abstract Implementation!"))
end

struct FakeMethod{T}
      obj::T
      FakeMethod(obj::T) where T = new{T}(obj)
      (x::FakeMethod)(args...) = x.obj
end

function findOrIsMethod(CheckMod::Module, MethodOrName)
      MethodOrName isa Method && return FakeMethod(MethodOrName)
      return isdefined(CheckMod, MethodOrName) ? getproperty(CheckMod, MethodOrName) : nothing
end

function findOrIsMethod(CheckMods::Tuple{Module}, MethodOrName)
      MethodOrName isa Method && return MethodOrName
      for mod in CheckMods
            isdefined(mod, MethodOrName) && return getproperty(mod, MethodOrName)
      end
      return nothing
end

macro proto(expr)
      Expr(:(=), Expr(:call, expr), Expr(:block))
end

macro nullc(value, if_null)
      return esc(:($value == nothing ? $if_null : $value))
end

function eqerror(error_message, nums...)
      if length(nums) == 0; return; end
      num = nums[1]
      for number in 2:nums
            if num != number; error(error_message); end
      end
end

@inline function getbit(byte::UInt8, index)::Bool
      return (byte & (UInt8(1) << index)) != 0
end

@inline function setbit(byte::UInt8, index, value::Bool)::UInt8
      return value ? (byte | (UInt8(1) << index)) : (byte & ~(UInt8(1) << index))
end

@inline function bytesof(v::T)::Array{UInt8, 1} where T
      return reinterpret(UInt8, isa(v, AbstractArray) ? v : T[v])
end

function bitsof(v::T)::BitArray where T
      bytes = bytesof(v)
      byte_copy = zeros(UInt8, length(bytes) + (8 - length(bytes) % 8))
      copyto!(byte_copy, bytes)
      bits = BitArray(undef, length(bytes) * 8)
      bits.chunks = reinterpret(UInt64, byte_copy)
      bits
end

function matrixbitsof(v::AbstractArray{T})::BitMatrix where T
      array = BitMatrix(undef, length(v), sizeof(T) * 8)
      for i in 1:length(v)
           array[i, 1:end] = bitsof(v[i])[1:(sizeof(T) * 8)]
      end
      array
end

function matrixbytesof(v::AbstractArray{T})::Matrix{UInt8} where T
      array = zeros(UInt8, length(v), sizeof(T))
      for i in 1:length(v)
           array[i, 1:end] = bytesof(v[i])
      end
      array
end

@inline function objectof(T::Type, bytes::AbstractArray{UInt8})
      obj_arr = reinterpret(T, bytes)
      return T == typeof(obj_arr) ? obj_arr : obj_arr[1]
end

function objectof(T::Type, bits::BitArray)
      bytes = reinterpret(UInt8, bits.chunks)
      object_bytes = zeros(UInt8, sizeof(T))
      copyto!(object_bytes, 1, bytes, 1, length(object_bytes))
      bytes = object_bytes
      objectof(T, bytes)
end

function matrixobjectof(T::Type, v::BitMatrix)
      array = Array{T}(undef, size(v, 1))
      for i in 1:size(v, 1)
            array[i] = objectof(T, v[i, 1:end])
      end
      array
end

function matrixobjectof(T::Type, v::AbstractMatrix{UInt8})
      array = Array{T}(undef, size(v, 1))
      for i in 1:size(v, 1)
            array[i] = objectof(T, v[i, 1:end])
      end
      array
end

function arrayobjectof(T::Type, v::AbstractArray{Array{UInt8, 1}})
      [objectof(T, data) for data in v]
end

function arrayobjectof(T::Type, v::AbstractArray{BitArray})
      [objectof(T, data) for data in v]
end

function arraybytesof(v::AbstractArray{T})::Array{Array{UInt8, 1}, 1} where T
      [bytesof(data) for data in v]
end

function arraybitsof(v::AbstractArray{T})::Array{BitArray, 1} where T
      [bitsof(data) for data in v]
end

function issupertype(type, check)::Bool
 while type != Any
      if type == check; return true; end
      type = supertype(type)
 end
 return false
end


text:
export
    particlesmoother_ffbs,
    particlesmoother_lffbs,
    particlesmoother_bbis,
    particlesmoother_lbbis,
    particlesmoother_llbbis,
    particlesmoother_fearnhead_lg

"""
    particlesmoother_ffbs

Particle smoother based on the Forward Filtering Backward Smoothing algorithm.
The complexity of this algorithm is O(KN^2)
"""
function particlesmoother_ffbs(hmm::HMM, psf::ParticleSet)
    K = length(psf)
    N = length(psf.p[1])
    # particle set smoother (storage)
    psw = deepcopy(psf)
    #
    pk = psw.p[K]

    for k=(K-1):-1:1
        pkp1 = pk
        pk   = psw.p[k]
        # denominator precomputation
        # -- complexity O(N^2) (linear storage)
        ds = [ sum( pk.w[l] *
                      exp(hmm.transloglik(k,pk.x[l],pkp1.x[j]))
                        for l in 1:N ) for j in 1:N ]
        # FFBS formula, in place computations
        # -- complexity O(N^2)
        for i in 1:N
            pk.w[i] *= sum( pkp1.w[j] *
                              exp(hmm.transloglik(k,pk.x[i],pkp1.x[j])) / ds[j]
                                for j in 1:N )
        end
        # normalisation
        pk.w /= sum(pk.w)
        # store the updated weights
        psw.p[k].w = copy(pk.w)
    end
    psw
end

"""
    particlesmoother_lffbs

Particle smoother based on the Forward Filtering Backward Smoothing algorithm.
The complexity of this algorithm is O(KNM)
"""
function particlesmoother_lffbs(hmm::HMM, psf::ParticleSet, M::Int)
    K = length(psf)
    N = length(psf.p[1])
    # particle set smoother (storage)
    psw = deepcopy(psf)
    #
    pk = psw.p[K]

    for k=(K-1):-1:1
        pkp1 = pk
        pk   = psw.p[k]

        r = rand(Multinomial(M, pkp1.w))
        r = [i for k in 1:N for i in ones(Int,r[k])*k]

        # denominator precomputation
        # -- complexity O(NM) (linear storage)
        ds = [ sum( pk.w[l] *
                      exp(hmm.transloglik(k,pk.x[l],pkp1.x[r[j]]))
                        for l in 1:N ) for j in 1:M ]
        # FFBS formula, in place computations
        # -- complexity O(NM)
        for i in 1:N
            pk.w[i] *= sum( pkp1.w[r[j]] *
                              exp(hmm.transloglik(k,pk.x[i],pkp1.x[r[j]])) / ds[j]
                                for j in 1:M )
        end
        # normalisation
        pk.w /= sum(pk.w)
        # store the updated weights
        psw.p[k].w = copy(pk.w)
    end
    psw
end

"""
  particlesmoother_bbis

Bootstrap backward information smoother. (CONSISTENT)
"""
function particlesmoother_bbis(hmm::HMM, observations::Matrix{Float},
                               psf::ParticleSet, bootstrap::Proposal;
                               resampling::Function=multinomialresampling,
                               essthresh::Float=0.5
                               )::Tuple{ParticleSet,Vector{Float}}
    K = length(psf)
    N = length(psf.p[1])
    # particle set smoother (storage)
    pss = ParticleSet(N, hmm.dimx, K)
    ess = zeros(K)
    # Particles at last step need to be resampled
    (pK,eK)  = resample(psf.p[K], essthresh, resampling, N)
    pss.p[K] = pK
    ess[K]   = eK
    #
    for k=(K-1):-1:2
        obsk   = observations[:,k]
        psskp1 = pss.p[k+1] # smoothing particles from previous step (k+1)
        psfk   = psf.p[k]   # filtering particles for PD_{k+1}(x_k+1)
        psfkm1 = psf.p[k-1] # filtering particles for PD_{k}(x_k)
        # preparing indices (for each j sample from 1 of the mixture component)
        randmult = rand(Multinomial(N, psfkm1.w))
        indices  = [i for k in 1:N for i in ones(Int,randmult[k])*k] # unroll
        # precompute denominator of the update factor PD_{k+1}(x_k+1)
        # -- complexity O(N^2) (linear storage)
        denj = [ sum( psfk.w[l] *
                        exp(hmm.transloglik(l, psfk.x[l], psskp1.x[j]))
                          for l in 1:N ) for j in 1:N ]
        xk    = similar(psskp1.x)
        logak = zeros(N)
        # sample and update weight for each smoothing particle
        # -- complexity O(N)
        for j = 1:N
            # sampling from corresponding element (see multinomial step)
            xk[j] = bootstrap.mean(k, psfkm1.x[indices[j]]) + bootstrap.noise()
            # weight update factor
            logak[j] = hmm.transloglik(k, xk[j], psskp1.x[j]) +
                        hmm.obsloglik(k, obsk, xk[j]) -
                          log(denj[j])
        end
        # normalise weights
        # -- complexity O(N)
        Wk  = log.(psskp1.w) + logak
        Wk -= minimum(Wk)
        wk  = exp.(Wk)
        wk /= sum(wk)

        (pk, ek) = resample(Particles(xk, wk), essthresh, resampling)

        pss.p[k] = pk
        ess[k]   = ek
    end

    # ----------------------------------------
    # LAST STEP (sampling from gamma1 = prior)

    obsk   = observations[:,1]
    psskp1 = pss.p[2] # smoothing particles from previous step (k+1)
    psfk   = psf.p[1]   # filtering particles for PD_{k+1}(x_k+1)
    denj   = [ sum( psfk.w[l] *
                      exp(hmm.transloglik(l, psfk.x[l], psskp1.x[j]))
                        for l in 1:N ) for j in 1:N ]
    xk     = similar(psskp1.x)
    logak  = zeros(N)

    for j = 1:N
        xk[j]    = bootstrap.mu0 + bootstrap.noise() # sampling from prior
        logak[j] = hmm.transloglik(1, xk[j], psskp1.x[j]) +
                    hmm.obsloglik(1, obsk, xk[j]) -
                      log(denj[j])
    end
    # normalise weights
    Wk  = log.(psskp1.w) + logak
    Wk -= minimum(Wk)
    wk  = exp.(Wk)
    wk /= sum(wk)

    (pss.p[1], ess[1]) = resample(Particles(xk,wk), essthresh, resampling)

    (pss, ess)
end

"""
  particlesmoother_lbbis

Bootstrap backward information smoother with linear complexity (APPROX BBIS).
"""
function particlesmoother_lbbis(hmm::HMM, observations::Matrix{Float},
                                psf::ParticleSet, bootstrap::Proposal;
                                resampling::Function=multinomialresampling,
                                essthresh::Float=0.5
                                )::Tuple{ParticleSet,Vector{Float}}
    K = length(psf)
    N = length(psf.p[1])
    # particle set smoother (storage)
    pss = ParticleSet(N, hmm.dimx, K)
    ess = zeros(K)
    # Particles at last step need to be resampled
    (pK,eK)  = resample(psf.p[K], essthresh, resampling, N)
    pss.p[K] = pK
    ess[K]   = eK
    #
    for k=(K-1):-1:2
        obsk   = observations[:,k]
        psskp1 = pss.p[k+1] # smoothing particles from previous step (k+1)
        psfk   = psf.p[k]   # filtering particles for PD_{k+1}(x_k+1)
        psfkm1 = psf.p[k-1] # filtering particles for PD_{k}(x_k)

        # preparing indices for the forward part
        randmult_fwd = rand(Multinomial(N, psfkm1.w))
        randmult_bwd = rand(Multinomial(N, psskp1.w))
        # unroll indices into masks
        indices_fwd = [i for k in 1:N for i in ones(Int,randmult_fwd[k])*k]
        indices_bwd = [i for k in 1:N for i in ones(Int,randmult_bwd[k])*k]

        xk    = similar(psskp1.x)
        logak = zeros(N)
        # sample and update weight for each smoothing particle
        for j = 1:N
            # sampling from corresponding element (see multinomial step)
            xk[j] = bootstrap.mean(k, psfkm1.x[indices_fwd[j]]) +
                      bootstrap.noise()
            # weight update factor p(x_t+1|x_t)p(y_t|x_t)
            logak[j] = hmm.transloglik(k, xk[j], psskp1.x[indices_bwd[j]]) +
                         hmm.obsloglik(k, obsk, xk[j])
        end
        # normalise weights
        Wk  = log.(psskp1.w) + logak
        Wk -= minimum(Wk)
        wk  = exp.(Wk)
        wk /= sum(wk)

        (pk, ek) = resample(Particles(xk, wk), essthresh, resampling)

        pss.p[k] = pk
        ess[k]   = ek
    end

    # ----------------------------------------
    # LAST STEP (sampling from gamma1 = prior)

    obsk   = observations[:,1]
    psskp1 = pss.p[2] # smoothing particles from previous step (k+1)

    randmult_bwd = rand(Multinomial(N, psskp1.w ))
    indices_bwd  = [i for k in 1:N for i in ones(Int,randmult_bwd[k])*k]

    xk     = similar(psskp1.x)
    logak  = zeros(N)

    for j = 1:N
        xk[j]    = bootstrap.mu0 + bootstrap.noise() # sampling from prior
        logak[j] = hmm.transloglik(1, xk[j], psskp1.x[indices_bwd[j]]) +
                    hmm.obsloglik(1, obsk, xk[j])
    end
    # normalise weights
    Wk  = log.(psskp1.w) + logak
    Wk -= minimum(Wk)
    wk  = exp.(Wk)
    wk /= sum(wk)

    (pss.p[1], ess[1]) = resample(Particles(xk,wk), essthresh, resampling)

    (pss, ess)
end

"""
  particlesmoother_llbbis

Bootstrap backward information smoother with log-log resampling (CONSISTENT)
"""
function particlesmoother_llbbis(hmm::HMM, observations::Matrix{Float},
                                 psf::ParticleSet, M::Int, bootstrap::Proposal;
                                 resampling::Function=multinomialresampling,
                                 essthresh::Float=0.5
                                 )::Tuple{ParticleSet,Vector{Float}}
    K = length(psf)
    N = length(psf.p[1])
    # particle set smoother (storage)
    pss = ParticleSet(N, hmm.dimx, K)
    ess = zeros(K)
    # Particles at last step need to be resampled
    (pK,eK)  = resample(psf.p[K], essthresh, resampling, N)
    pss.p[K] = pK
    ess[K]   = eK
    # force sub-resampling of forward components
    (psfkm1, e) = resample(psf.p[K-1], Inf, resampling, M)
    #
    for k=(K-1):-1:2
        obsk   = observations[:,k]
        # smoothing particles from previous step (k+1) + force sub-resampling
        psskp1 = pss.p[k+1]
        # filtering particles for PD_{k+1}(x_k+1) + recycled from previous step
        psfk   = psfkm1
        # filtering particles for PD_{k}(x_k)
        (psfkm1, e) = resample(psf.p[k-1], Inf, resampling, M)
        # preparing indices (for each j sample from 1 of the mixture component)
        # since we've done the resampling all mixture components have equal w
        randmult = rand(Multinomial(N, M))
        indices  = [i for k in 1:M for i in ones(Int,randmult[k])*k] # unroll
        # precompute denominator of the update factor PD_{k+1}(x_k+1)
        # -- complexity O(N*M) (N storage)
        denj = [ sum( psfk.w[l] *
                        exp(hmm.transloglik(l, psfk.x[l], psskp1.x[j]))
                          for l in 1:M ) for j in 1:N ]
        xk    = similar(psskp1.x)
        logak = zeros(N)
        # sample and update weight for each smoothing particle
        # -- complexity O(N)
        for j = 1:N
            # sampling from corresponding element (see multinomial step)
            xk[j] = bootstrap.mean(k, psfkm1.x[indices[j]]) + bootstrap.noise()
            # weight update factor
            logak[j] = hmm.transloglik(k, xk[j], psskp1.x[j]) +
                        hmm.obsloglik(k, obsk, xk[j])
        end
        # normalise weights
        # -- complexity O(N)
        Wk  = log.(psskp1.w) + logak - log.(denj)
        Wk -= minimum(Wk)
        wk  = exp.(Wk)
        wk /= sum(wk)

        (pk, ek) = resample(Particles(xk, wk), essthresh, resampling)

        pss.p[k] = pk
        ess[k]   = ek
    end

    # ----------------------------------------
    # LAST STEP (sampling from gamma1 = prior)

    obsk   = observations[:,1]
    psskp1 = pss.p[2] # smoothing particles from previous step (k+1)
    psfk   = psfkm1   # filtering particles for PD_{k+1}(x_k+1)
    denj   = [ sum( psfk.w[l] *
                      exp(hmm.transloglik(l, psfk.x[l], psskp1.x[j]))
                        for l in 1:M ) for j in 1:N ]
    xk     = similar(psskp1.x)
    logak  = zeros(N)

    for j = 1:N
        xk[j]    = bootstrap.noise() # sampling from prior
        logak[j] = hmm.transloglik(1, xk[j], psskp1.x[j]) +
                    hmm.obsloglik(1, obsk, xk[j]) -
                      log(denj[j])
    end
    # normalise weights
    Wk  = log.(psskp1.w) + logak
    Wk -= minimum(Wk)
    wk  = exp.(Wk)
    wk /= sum(wk)

    (pss.p[1], ess[1]) = resample(Particles(xk,wk), essthresh, resampling)

    (pss, ess)
end

"""
    particlesmoother_fearnhead_lg

Linear complexity smoother a la Fearnhead with normalising densities obtained
via the simple dynamic. This is just for reproduction purposes and should
not really be used.
"""
function particlesmoother_fearnhead_lg(lg::LinearGaussian,
                                       observations::Matrix{Float},
                                       psf::ParticleSet;
                                       resampling=multinomialresampling,
                                       essthresh::Float=0.5)::ParticleSet
    K = length(psf)
    N = length(psf.p[1])

    Qi   = inv(lg.Q)
    QiA  = Qi * lg.A
    AQiA = lg.A' * QiA

    Ri   = inv(lg.R)
    RiB  = Ri * lg.B
    BRiB = lg.B' * RiB

    ############################################
    # Computation of the means and covariances
    # of the normalising densities γ_k
    ############################################
    # initialise
    gamma_mu = zeros(lg.dimx, K)
    gamma_S  = zeros(lg.dimx, lg.dimx, K)

    gamma_mu[:, 1] = zeros(lg.dimx)    # this could be passed to function
    gamma_S[:,:,1] = 0.2*eye(lg.dimx)  # this could be passed to function

    for k = 2:K
        # Cf thesis appendix for full development, just compl of quadratics
        temp_S         = inv(inv(gamma_S[:,:,k-1]) + AQiA)
        gamma_S[:,:,k] = inv(Qi - QiA * temp_S * QiA')
        temp_mu        = gamma_S[:,:,k-1] \ gamma_mu[:,k-1]
        gamma_mu[:,k]  = gamma_S[:,:,k] * QiA * temp_S * temp_mu
    end
    # ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

    ############################################
    # Targeting the Backward Information Filter
    # with the normalising densities computed
    # above. Since everything is tractable,
    # proposal is the optimal proposal.
    ############################################

    bif = ParticleSet(N, lg.dimx, K)

    # initial: sample from γ_K p(y_K | x_K), completion of squares again

    yK      = observations[:,K]
    bif_SK  = inv( inv(gamma_S[:,:,K]) + BRiB )
    bif_muK = bif_SK * ( (gamma_S[:,:,K] \ gamma_mu[:, K]) + RiB' * yK )

    noise    = chol(Symmetric(bif_SK))' * randn(lg.dimx, N)
    bif.p[K] = Particles( [ bif_muK + noise[:,j] for j in 1:N] ,
                          ones(N)/N)

    mvnpdf(x,mu,S) = exp(-0.5dot(S\(x-mu),x-mu))
    #
    for k = K-1:-1:1
        # observation
        yk = observations[:,k]

        # construct the gaussian γ_t * p(x_k+1 | x_k) p(y_k | x_k)
        # > covariance
        optpropS   = Symmetric(inv( inv(gamma_S[:,:,k]) + BRiB + AQiA ))
        # > mean (fixed part not depending on x_k+1)
        optpropmu_ = optpropS * ( gamma_S[:,:,k] \ gamma_mu[:,k] + RiB' * yk )
        # > mean (variable part, depending on x_k+1)
        optpropmu_j = [optpropS * QiA' * bif.p[k+1].x[j] for j in 1:N]

        # sample from optimal prop
        noise   = chol(optpropS)' * randn(lg.dimx, N)
        samples = [optpropmu_ + optpropmu_j[j] + noise[:,j] for j in 1:N]

        # compute γ_k+1 for each x_k+1
        gammakp1 = [ mvnpdf( bif.p[k+1].x[j], gamma_mu[:,k+1],
                       gamma_S[:,:,k+1]) for j in 1:N ]

        # compute the weights (previous weights div by gammakp1)
        weights  = bif.p[k+1].w ./ gammakp1
        # > normalise
        weights /= sum(weights)
        # > resample
        (bif.p[k], ek) = resample(Particles(samples, weights),
                                    essthresh, resampling)
    end

    # combination
    pss = ParticleSet(N, lg.dimx, K)
    pss.p[1] = deepcopy(bif.p[1])

    propS  = Symmetric(inv(Qi + AQiA + BRiB))
    cpropS = chol(propS)
    for k = 2:K-1
        # observation
        yk = observations[:, k]
        # preparing indices for the forward part
        randmult_fwd = rand(Multinomial(N, psf.p[k-1].w))
        randmult_bwd = rand(Multinomial(N, bif.p[k].w))
        # ^ note bif.p[k] and not k+1 since kth weight is scaled by γ_t+1
        # unroll indices into masks
        indices_fwd = [i for k in 1:N for i in ones(Int,randmult_fwd[k])*k]
        indices_bwd = [i for k in 1:N for i in ones(Int,randmult_bwd[k])*k]
        # # get samples in one shot
        noise = cpropS' * randn(lg.dimx, N)
        mus = zeros(lg.dimx, N)
        for j = 1:N
            mus[:, j] = QiA' * bif.p[k+1].x[indices_bwd[j]] +
                          QiA * psf.p[k-1].x[indices_fwd[j]] + RiB' * yk
        end
        mus = propS * mus + noise
        xk  = [mus[:, j] for j in 1:N]
        # no resampling since perfect sampling
        pss.p[k] = Particles(xk, ones(N)/N)
    end
    pss.p[K] = deepcopy(psf.p[K])
    pss
end

# """
#     particle_fbpf
#
# Forward Backward Particle Filter
# """
# function particlesmoother_fbpf(hmm::HMM, observations::Matrix{Float},
#                                N::Int, M::Int,
#                                propfwd::Proposal, propbwd::Proposal,
#                                npasses::Int=3 ;
#                                resampling::Function=multinomialresampling,
#                                essthresh::Float=0.5
#                                )::Tuple{ParticleSet,Vector{Float}}
#     # PASS 1 = particle filter
#     (pf, ess) = particlefilter(hmm, observations, N, propfwd)
#     # PASS 2 = filter using PF as guide
#     (ps, ess) = particlesmoother_llbbis(hmm, observations, pf, M, propfwd)
#     # PASS 3 = filter using PS as guide etc...
#     K = length(pf)
#
#     pss = ParticleSet(N, hmm.dimx, K)
#     ess = zeros(K)
#
#     (p1,e1) = resample(
#                 Particles( [proposal.noise() for i in 1:N], ones(N)/N),
#                 essthresh )
#
#     # store
#     pss.p[1] = p1
#     ess[1]   = e1
#
#     for k=2:K-1
#         pkp1  = ps.p[k+1]
#         obsk  = observations[:,k]
#         logak = zeros(N)
#         xk    = similar(pkp1.x)
#         # sample from multinomial
#         r     = rand(Multinomial())
#
#     end
#
# end


text:
"""
    SynchronousExecution

Functor object for using synchronous execution in a progressive-hedging algorithm (assuming multiple Julia cores are available). Create by supplying a [`Synchronous`](@ref) object through `execution` in the `ProgressiveHedgingSolver` factory function and then pass to a `StochasticPrograms.jl` model.

"""
struct SynchronousExecution{T <: AbstractFloat,
                            A <: AbstractVector,
                            PT <: AbstractPenaltyterm} <: AbstractProgressiveHedgingExecution
    subworkers::Vector{SubWorker{T,A,PT}}

    function SynchronousExecution(::Type{T}, ::Type{A}, ::Type{PT}) where {T <: AbstractFloat,
                                                                           A <: AbstractVector,
                                                                           PT <: AbstractPenaltyterm}
        return new{T,A,PT}(Vector{SubWorker{T,A,PT}}(undef, nworkers()))
    end
end

function initialize_subproblems!(ph::AbstractProgressiveHedging,
                                 execution::SynchronousExecution,
                                 scenarioproblems::DistributedScenarioProblems,
                                 penaltyterm::AbstractPenaltyterm)
    # Create subproblems on worker processes
    initialize_subproblems!(ph,
                            execution.subworkers,
                            scenarioproblems,
                            penaltyterm)
    # Initial reductions
    update_iterate!(ph)
    update_dual_gap!(ph)
    return nothing
end

function restore_subproblems!(::AbstractProgressiveHedging, execution::SynchronousExecution)
    restore_subproblems!(execution.subworkers)
    return nothing
end

function resolve_subproblems!(ph::AbstractProgressiveHedging, execution::SynchronousExecution{T}) where T <: AbstractFloat
    partial_solutions = Vector{SubproblemSolution{T}}(undef, nworkers())
    @sync begin
        for (i,w) in enumerate(workers())
            @async partial_solutions[i] = remotecall_fetch(resolve_subproblems!,
                                                           w,
                                                           execution.subworkers[w-1],
                                                           ph.ξ,
                                                           penalty(ph))
        end
    end
    return sum(partial_solutions)
end

function update_iterate!(ph::AbstractProgressiveHedging, execution::SynchronousExecution{T,A}) where {T <: AbstractFloat, A <: AbstractVector}
    partial_primals = Vector{A}(undef, nworkers())
    @sync begin
        for (i,w) in enumerate(workers())
            @async partial_primals[i] = remotecall_fetch(collect_primals,
                                                         w,
                                                         execution.subworkers[w-1],
                                                         length(ph.ξ))
        end
    end
    ξ_prev = copy(ph.ξ)
    ph.ξ .= sum(partial_primals)
    # Update δ₁
    ph.data.δ₁ = norm(ph.ξ - ξ_prev, 2) ^ 2
    return nothing
end

function update_subproblems!(ph::AbstractProgressiveHedging, execution::SynchronousExecution)
    # Update dual prices
    @sync begin
        for w in workers()
            @async remotecall_fetch(
                w,
                execution.subworkers[w-1],
                ph.ξ,
                penalty(ph)) do sw, ξ, r
                    subproblems = fetch(sw)
                    if length(subproblems) > 0
                        update_subproblems!(subproblems, ξ, r)
                    end
                end
        end
    end
    return nothing
end

function update_dual_gap!(ph::AbstractProgressiveHedging, execution::SynchronousExecution)
    return update_dual_gap!(ph, execution.subworkers)
end

function calculate_objective_value(ph::AbstractProgressiveHedging, execution::SynchronousExecution)
    return calculate_objective_value(execution.subworkers)
end

function scalar_subproblem_reduction(value::Function, execution::SynchronousExecution{T}) where T <: AbstractFloat
    partial_results = Vector{T}(undef, nworkers())
    @sync begin
        for (i,w) in enumerate(workers())
            @async partial_results[i] = remotecall_fetch(w, execution.subworkers[w-1], value) do sw, value
                subproblems = fetch(sw)
                return mapreduce(+, subproblems, init = zero(T)) do subproblem
                    π = subproblem.probability
                    return π * value(subproblem)
                end
            end
        end
    end
    return sum(partial_results)
end

function vector_subproblem_reduction(value::Function, execution::SynchronousExecution{T,A}, n::Integer) where {T <: AbstractFloat, A <: AbstractVector}
    partial_results = Vector{A}(undef, nworkers())
    @sync begin
        for (i,w) in enumerate(workers())
            @async partial_results[i] = remotecall_fetch(w, execution.subworkers[w-1], value, n) do sw, value, n
                subproblems = fetch(sw)
                return mapreduce(+, subproblems, init = zero(T, n)) do subproblem
                    π = subproblem.probability
                    return π * value(subproblem)
                end
            end
        end
    end
    return sum(partial_results)
end

# API
# ------------------------------------------------------------
function (execution::Synchronous)(::Type{T}, ::Type{A}, ::Type{PT}) where {T <: AbstractFloat,
                                                                           A <: AbstractVector,
                                                                           PT <: AbstractPenaltyterm}
    return SynchronousExecution(T,A,PT)
end

function str(::Synchronous)
    return "Synchronous "
end


text:
# TODO: turn this into an ITensors.jl function `originalsiteinds`
# that generically returns the site indices that would be used to
# make an object of the same type with the same indices.
hilbertspace(M::MPS)  = siteinds(first, M; plev=0)
hilbertspace(M::MPO)  = dag.(siteinds(first, M; plev=0))
hilbertspace(L::LPDO) = dag.(siteinds(first, L.X; plev=0, tags=!purifier_tags(L)))

"""
    convertdatapoint(datapoint::Array,basis::Array;state::Bool=false)

0 1 0 0 1 -> Z+ Z- Z+ Z+ Z-
"""
function convertdatapoint(datapoint::Array{Int64}; state::Bool=false)
  newdata = []
  basis = ["Z" for _ in 1:length(datapoint)]
  for j in 1:length(datapoint)
    if datapoint[j] == 0
      push!(newdata, "Z+")
    elseif datapoint[j] == 1
      push!(newdata, "Z-")
    else
      error("non-binary data")
    end
  end
  return newdata
end

"""

(Z+, X-) -> (Z => 0), (X => 1)

"""
function convertdatapoint(datapoint::Array{String})
  basis = []
  outcome = []
  for j in 1:length(datapoint)
    push!(basis, string(datapoint[j][1]))

    if datapoint[j][2] == Char('+')
      push!(outcome, 0)
    elseif datapoint[j][2] == Char('-')
      push!(outcome, 1)
    else
      error("non-binary data")
    end
  end
  return basis .=> outcome
end

function convertdatapoints(datapoints::Array{String})
  nshots = size(datapoints)[1]
  newdata = Matrix{Pair{String,Int64}}(undef, nshots, size(datapoints)[2])

  for n in 1:nshots
    newdata[n, :] = convertdatapoint(datapoints[n, :])
  end
  return newdata
end

"""

(Z, 0), (X, 1)  / (Z+, X-)

"""
function convertdatapoint(outcome::Array{Int64}, basis::Array{String}; state::Bool=false)
  @assert length(outcome) == length(basis)
  newdata = []
  if state
    basis = basis
  end
  for j in 1:length(outcome)
    if outcome[j] == 0
      push!(newdata, basis[j] * "+")
    elseif outcome[j] == 1
      push!(newdata, basis[j] * "-")
    else
      error("non-binary data")
    end
  end
  return newdata
end

function convertdatapoint(datapoint::Array{Pair{String,Int64}}; state::Bool=false)
  return convertdatapoint(last.(datapoint), first.(datapoint); state=state)
end

function convertdatapoints(datapoints::Array{Pair{String,Int64}}; state::Bool=false)
  nshots = size(datapoints)[1]
  newdata = Matrix{String}(undef, nshots, size(datapoints)[2])

  for n in 1:nshots
    newdata[n, :] = convertdatapoint(datapoints[n, :]; state=state)
  end
  return newdata
end

function convertdatapoints(outcome::Matrix{Int64}, basis::Matrix{String}; state::Bool=false)
  nshots = size(datapoints)[1]
  newdata = Matrix{String}(undef, nshots, size(datapoints)[2])

  for n in 1:nshots
    newdata[n, :] = convertdatapoint(datapoints[n, :]; state=state)
  end
  return newdata
end

"""

(Z => 0), (X => 1)  / (Z+, X-)

"""

"""
    split_dataset(data::Matrix; train_ratio::Float64 = 0.9, randomize::Bool = true)

Split a data set into a `train` and `test` sets, given a `train_ratio` (i.e. the
percentage of data in `train_data`. If `randomize=true` (default), the `data` is
randomly shuffled before splitting.
"""
function split_dataset(data::Matrix; train_ratio::Float64=0.9, randomize::Bool=true)
  ndata = size(data, 1)
  ntrain = Int(ndata * train_ratio)
  ntest = ndata - ntrain
  if randomize
    data = data[shuffle(1:end), :]
  end
  train_data = data[1:ntrain, :]
  test_data = data[(ntrain + 1):end, :]
  return train_data, test_data
end

"""
    ischoi(M::LPDO)

Check whether a given LPDO{MPO}
"""
function ischoi(M::LPDO{MPO})
  !haschoitags(M) && return false
  for j in 1:length(M)
    length(inds(M.X[j], "Site")) != 2 && return false
  end
  return true
end

function ischoi(M::MPO)
  !haschoitags(M) && return false
  for j in 1:length(M)
    length(inds(M[j], "Site")) != 4 && return false
  end
  return true
end

ischoi(M::ITensor) = haschoitags(M)

"""
    haschoitags(L::LPDO)
    haschoitags(M::Union{MPS,MPO})

Check whether the TN has input/output Choi tags
"""
function haschoitags(L::LPDO)
  for j in 1:length(L)
    !(hastags(L.X[j], "Input") && hastags(L.X[j], "Output")) && return false
  end
  return true
end
function haschoitags(M::Union{MPS,MPO})
  for j in 1:length(M)
    !(hastags(M[j], "Input") && hastags(M[j], "Output")) && return false
  end
  return true
end

function haschoitags(M::ITensor)
  N = nsites(M)
  for j in 1:N
    !(hastags(M,"Input,n=$j") && hastags(M,"Output,n=$j")) && return false
  end
  return true
end

"""
    choitags(U::MPO)

Assign the input/output tags defined for a Choi matrix to an MPO.

  σ₁ -o- σ₁′       σ₁ⁱ -o- σ₁ᴼ
      |                 |
  σ₂ -o- σ₂′  ⟶    σ₂ⁱ -o- σ₂ᴼ
      |                 |
  σ₃ -o- σ₃′       σ₃ⁱ -o- σ₃ᴼ

"""
function choitags(U::MPO)
  haschoitags(U) && return U
  U = addtags(siteinds, U, "Input"; plev=0)
  U = addtags(siteinds, U, "Output"; plev=1)
  return noprime(siteinds, U)
end

function choitags(T::ITensor)
  haschoitags(T) && return T
  T = addtags(T, "Input"; plev=0)
  T = addtags(T, "Output"; plev=1)
  return noprime(T)
end

choitags(O::ITensorOperator) =
  ITensorState(choitags(O.T))

"""
    mpotags(U::MPO)

Inverse of choitags.
"""
function mpotags(U::MPO)
  U = prime(U; tags="Output")
  U = removetags(U, "Input")
  U = removetags(U, "Output")
  return U
end

mpotags(M::Union{MPS,MPO}) = mpotags(LPDO(M)).X

"""
    unitary_mpo_to_choi_mps(U::MPO)

     MPO          MPS (vectorized MPO)
  σ₁ -o- σ₁′       o= (σ₁ⁱ,σ₁ᴼ)
      |            |
  σ₂ -o- σ₂′  ⟶    o= (σ₂ⁱ,σ₂′ᴼ)
      |            |
  σ₃ -o- σ₃′       o= (σ₃ⁱ,σ₃′ᴼ)

Transforms a unitary MPO into a Choi MPS with appropriate tags.
"""
unitary_mpo_to_choi_mps(U::MPO) = convert(MPS, choitags(U))
unitary_mpo_to_choi_mps(T::ITensor) = choitags(T)
unitary_mpo_to_choi_mps(L::LPDO{MPO}) = unitary_mpo_to_choi_mps(L.X)

"""
    unitary_mpo_to_choi_mpo(U::MPO)

     MPO                   MPO
  σ₁ -o- σ₁′     (σ₁ⁱ,σ₁ᴼ) =o= (σ₁′ⁱ,σ₁′ᴼ)
      |                     |
  σ₂ -o- σ₂′  ⟶  (σ₂ⁱ,σ₂ᴼ) =o= (σ₂′ⁱ,σ₂′ᴼ)
      |                     |
  σ₃ -o- σ₃′     (σ₃ⁱ,σ₃ᴼ) =o= (σ₃′ⁱ,σ₃′ᴼ)

Convert a unitary MPO to a Choi matrix represented as an MPO with 4 site indices.
"""
unitary_mpo_to_choi_mpo(U::MPO) = MPO(LPDO(convert(MPS, choitags(U))))

"""
    choi_mps_to_unitary_mpo(Ψ::MPS)

Transforms a Choi MPS into an MPO with appropriate tags.
Inverse of `unitary_mpo_to_choi_mps`.
"""
choi_mps_to_unitary_mpo(Ψ::MPS) = mpotags(convert(MPO, Ψ))
choi_mps_to_unitary_mpo(L::LPDO{MPS}) = choi_mps_to_unitary_mpo(L.X)

function nsites(T::ITensor)
  s1 = inds(T,tags="Site,n=1")
  # Wavefunction
  if length(s1) == 1 || length(s1) == 2
    return length(inds(T,plev=0))
  # Choi matrix
  elseif length(s1) == 4
    return length(inds(T,plev=0)) ÷ 2
  else
    error("Indices not recognized")
  end
end

nsites(M::Union{MPS,MPO,LPDO}) =
  length(M)


text:
#=
 Example : reconstruct glacial oxygen
1. diagnose sat o2 percentage modern ocean
2. using sat o2 percentage and glacial T for glacial o2 boundary condition
3. propagate into interior with stoichiometric ratio
4. difference with modern
5. plot it
=#
using Revise
using TMI, BenchmarkTools, PyPlot, PyCall, GibbsSeaWater

modernversion = "modern_90x45x33_GH10_GH12"
A, Alu, γ, modernfile, L, B = config_from_nc(modernversion)

O₂solmodern,O₂fractionmodern = surface_oxygensaturation(modernfile)

# modern section
O₂modern = readtracer(modernfile,"O₂")
Osection = section(O₂modern,lon_section,γ)
titlelabel = "O₂ Modern "*string(lon_section)*"°E"
# make a plot of dye in the ocean
dyeplot(γ.lat,γ.depth,Osection,lims,titlelabel)
figname = plotsdir("oxygen_30W_modern.eps")
savefig(figname)

# next get LGM temperature
lgmnames = ("G14","G14A","GPLS1","GPLS2","OG18")
for lgmname in lgmnames

    LGMversion = "LGM_90x45x33_"*lgmname
    O₂lgm = oxygen(LGMversion,O₂fractionmodern)

    # plot a section at 330 east longitude (i.e., 30 west)
    lon_section = 330;
    lims = 50:10:350

    # LGM section
    Osection = section(O₂lgm,lon_section,γ)
    titlelabel = "O₂ "*lgmname*string(lon_section)*"°E"
    # make a plot of dye in the ocean
    dyeplot(γ.lat,γ.depth,Osection,lims,title)
    figname = plotsdir("oxygen_30W_"*lgmname*".eps")
    savefig(figname)

    # LGM - modern
    lims = -100:10:100
    Osection = section(O₂lgm - O₂modern,lon_section,γ)
    titlelabel = "O₂ "*lgmname*"-modern "*string(lon_section)*"°E"
    # make a plot of dye in the ocean
    dyeplot(γ.lat,γ.depth,Osection,lims,titlelabel)
    figname = plotsdir("oxygen_30W_"*lgmname*"-Modern.png")
    savefig(figname)

end


text:
@inline function setblockdiag!(l::Lattice, A::AbstractSparseMatrix, row::Int, col::Int, B::AbstractVector)
  rstart = (row-1)*l.sites+1
  colstart = (col-1)*l.sites+1

  @inbounds for shift in 0:length(B)-1
    A[rstart+shift,colstart+shift] = B[shift+1]
  end
  nothing
end

@inline function setblockdiag_minus!(l::Lattice, A::AbstractSparseMatrix, row::Int, col::Int, B::AbstractVector)
  rstart = (row-1)*l.sites+1
  colstart = (col-1)*l.sites+1

  @inbounds for shift in 0:length(B)-1
    A[rstart+shift,colstart+shift] = -B[shift+1] # minus here!
  end
  nothing
end

# interaction_matrix_exp = exp(- power delta_tau V(slice)), with power = +- 1.
function interaction_matrix_exp(mc::AbstractDQMC, slice::Int, power::Float64=1.)
  G = geltype(mc)

  eV = spzeros(G, mc.p.flv * mc.l.sites, mc.p.flv * mc.l.sites)
  interaction_matrix_exp!(mc, slice, power, eV)
  return eV
end

# interaction_matrix_exp = exp(- power delta_tau V(slice)), with power = +- 1.
function interaction_matrix_exp!(mc::AbstractDQMC, slice::Int, power::Float64=1., eV::AbstractSparseMatrix=mc.s.eV)
  @mytimeit mc.a.to "interactions (combined)" begin
  @mytimeit mc.a.to "interaction_matrix_exp!" begin
  p = mc.p
  l = mc.l
  N = l.sites
  G = geltype(mc)

  C = mc.s.C
  S = mc.s.S
  R = mc.s.R

  # C = zeros(G, N) #1, 1)
  # S = zeros(G, N) #1, 2)
  # if p.opdim == 3
    # R = zeros(G, N) #1, 4)
  # end
  @inbounds @simd for i in 1:N
    @views n = norm(p.hsfield[:,i,slice])
    sh = sinh(p.lambda * p.delta_tau * n)/n
    C[i] = cosh(p.lambda * p.delta_tau * n)
    if p.opdim == 3
      S[i] = (im * p.hsfield[2,i,slice] - p.hsfield[1,i,slice]) * power * sh
      R[i] = (-p.hsfield[3,i,slice]) * power * sh
    elseif p.opdim == 2
      S[i] = (im * p.hsfield[2,i,slice] - p.hsfield[1,i,slice]) * power * sh
    else # O(1)
      S[i] = - p.hsfield[1,i,slice] * power * sh
    end
  end

  setblockdiag!(l,eV,1,1,C)
  setblockdiag!(l,eV,1,2,S)

  cS = conj(S) # Once julia 1.2 is out, use lazy ConjVector here.
  setblockdiag!(l,eV,2,1,cS)
  setblockdiag!(l,eV,2,2,C)

  if p.opdim == 3
    setblockdiag!(l,eV,1,4,R)

    setblockdiag_minus!(l,eV,2,3,R)

    setblockdiag_minus!(l,eV,3,2,R)
    setblockdiag!(l,eV,3,3,C)
    setblockdiag!(l,eV,3,4,cS)

    setblockdiag!(l,eV,4,1,R)
    setblockdiag!(l,eV,4,3,S)
    setblockdiag!(l,eV,4,4,C)
  end
  end #timeit
  end #timeit
end

# calculate p.flv x p.flv (4x4 for O(3) model) interaction matrix exponential for given op
function interaction_matrix_exp_op(mc::AbstractDQMC, op::Vector{Float64}, power::Float64=1.)
  eVop = Matrix{geltype(mc)}(undef,mc.p.flv,mc.p.flv)
  interaction_matrix_exp_op!(mc,op,power,eVop)
  return eVop
end

# calculate p.flv x p.flv (4x4 for O(3) model) interaction matrix exponential for given op
function interaction_matrix_exp_op!(mc::AbstractDQMC{C,G}, op::Vector{Float64}, power::Float64=1., eVop::Matrix{G}=mc.s.eVop1) where {C,G}
  @mytimeit mc.a.to "interactions (combined)" begin
  @mytimeit mc.a.to "interaction_matrix_exp_op!" begin
  n = norm(op)
  sh = power * sinh(mc.p.lambda * mc.p.delta_tau*n)/n
  Cii = cosh(mc.p.lambda * mc.p.delta_tau*n)
  if mc.p.opdim == 3
    Sii = (im * op[2] - op[1]) * sh
    Rii = (-op[3]) * sh
  elseif mc.p.opdim == 2
    Sii = (im * op[2] - op[1]) * sh
  else # O(1)
    Sii = - op[1] * sh
  end

  eVop[1,1] = Cii
  eVop[1,2] = Sii
  eVop[2,1] = conj(Sii)
  eVop[2,2] = Cii

  if mc.p.opdim == 3
    eVop[1,3] = zero(G)
    eVop[1,4] = Rii

    eVop[2,3] = -Rii
    eVop[2,4] = zero(G)

    eVop[3,1] = zero(G)
    eVop[3,2] = -Rii
    eVop[3,3] = Cii
    eVop[3,4] = conj(Sii)

    eVop[4,1] = Rii
    eVop[4,2] = zero(G)
    eVop[4,3] = Sii
    eVop[4,4] = Cii
  end
  end #timeit
  end #timeit
end

max_stars_repo_name: crstnbr/dqmc

max_issues_repo_name: crstnbr/dqmc

max_forks_repo_name: crstnbr/dqmc


text:
################################################################################
#                            Planar and Radial Flows                           #
#             Ref: Variational Inference with Normalizing Flows,               #
#               D. Rezende, S. Mohamed(2015) arXiv:1505.05770                  #
################################################################################

###############
# RadialLayer #
###############

mutable struct RadialLayer{T1<:Union{Real, AbstractVector{<:Real}}, T2<:AbstractVector{<:Real}} <: Bijector{1}
    α_::T1
    β::T1
    z_0::T2
end
function Base.:(==)(b1::RadialLayer, b2::RadialLayer)
    return b1.α_ == b2.α_ && b1.β == b2.β && b1.z_0 == b2.z_0
end

function RadialLayer(dims::Int, wrapper=identity)
    α_ = wrapper(randn(1))
    β = wrapper(randn(1))
    z_0 = wrapper(randn(dims))
    return RadialLayer(α_, β, z_0)
end

# all fields are numerical parameters
Functors.@functor RadialLayer

h(α, r) = 1 ./ (α .+ r)     # for radial flow from eq(14)
#dh(α, r) = .- (1 ./ (α .+ r)) .^ 2   # for radial flow; derivative of h()

# An internal version of transform that returns intermediate variables
function _transform(flow::RadialLayer, z::AbstractVecOrMat)
    return _radial_transform(first(flow.α_), first(flow.β), flow.z_0, z)
end
function _radial_transform(α_, β, z_0, z)
    α = LogExpFunctions.log1pexp(α_)            # from A.2
    β_hat = -α + LogExpFunctions.log1pexp(β)    # from A.2
    if z isa AbstractVector
        r = norm(z .- z_0)
    else
        r = vec(sqrt.(sum(abs2, z .- z_0; dims = 1)))
    end
    transformed = z .+ β_hat ./ (α .+ r') .* (z .- z_0)   # from eq(14)
    return (transformed = transformed, α = α, β_hat = β_hat, r = r)
end

(b::RadialLayer)(z::AbstractMatrix{<:Real}) = _transform(b, z).transformed
(b::RadialLayer)(z::AbstractVector{<:Real}) = vec(_transform(b, z).transformed)

function forward(flow::RadialLayer, z::AbstractVecOrMat)
    transformed, α, β_hat, r = _transform(flow, z)
    # Compute log_det_jacobian
    d = size(flow.z_0, 1)
    h_ = h(α, r)
    if transformed isa AbstractVector
        T = eltype(transformed)
    else
        T = typeof(vec(transformed))
    end
    log_det_jacobian::T = @. (
        (d - 1) * log(1 + β_hat * h_)
        + log(1 +  β_hat * h_ + β_hat * (- h_ ^ 2) * r)
    )   # from eq(14)
    return (rv = transformed, logabsdetjac = log_det_jacobian)
end

function (ib::Inverse{<:RadialLayer})(y::AbstractVector{<:Real})
    flow = ib.orig
    z0 = flow.z_0
    α = LogExpFunctions.log1pexp(first(flow.α_))            # from A.2
    α_plus_β_hat = LogExpFunctions.log1pexp(first(flow.β))  # from A.2

    # Compute the norm ``r`` from A.2.
    y_minus_z0 = y .- z0
    r = compute_r(y_minus_z0, α, α_plus_β_hat)
    γ = (α + r) / (α_plus_β_hat + r)

    return z0 .+ γ .* y_minus_z0
end

function (ib::Inverse{<:RadialLayer})(y::AbstractMatrix{<:Real})
    flow = ib.orig
    z0 = flow.z_0
    α = LogExpFunctions.log1pexp(first(flow.α_))            # from A.2
    α_plus_β_hat = LogExpFunctions.log1pexp(first(flow.β))  # from A.2

    # Compute the norm ``r`` from A.2 for each column.
    y_minus_z0 = y .- z0
    rs = mapvcat(eachcol(y_minus_z0)) do c
        return compute_r(c, α, α_plus_β_hat)
    end
    γ = reshape((α .+ rs) ./ (α_plus_β_hat .+ rs), 1, :)

    return z0 .+ γ .* y_minus_z0
end

"""
    compute_r(y_minus_z0::AbstractVector{<:Real}, α, α_plus_β_hat)

Compute the unique solution ``r`` to the equation
```math
\\|y_minus_z0\\|_2 = r \\left(1 + \\frac{α_plus_β_hat - α}{α + r}\\right)
```
subject to ``r ≥ 0`` and ``r ≠ α``.

Since ``α > 0`` and ``α_plus_β_hat > 0``, the solution is unique and given by
```math
r = (\\sqrt{(α_plus_β_hat - γ)^2 + 4 α γ} - (α_plus_β_hat - γ)) / 2,
```
where ``γ = \\|y_minus_z0\\|_2``. For details see appendix A.2 of the reference.

# References

D. Rezende, S. Mohamed (2015): Variational Inference with Normalizing Flows.
arXiv:1505.05770
"""
function compute_r(y_minus_z0::AbstractVector{<:Real}, α, α_plus_β_hat)
    γ = norm(y_minus_z0)
    a = α_plus_β_hat - γ
    r = (sqrt(a^2 + 4 * α * γ) - a) / 2
    return r
end

logabsdetjac(flow::RadialLayer, x::AbstractVecOrMat) = forward(flow, x).logabsdetjac


text:
"""
    function gather_post_mean_infection_rates_by_type(model::KenyaCoVaccines.CoVAreaModel)

Gather the posterior mean predictions for infections by: age, dose, primary/reinfection episode, variant.
"""
function gather_post_mean_infection_rates_by_type(model::KenyaCoVaccines.CoVAreaModel, projection_date::Date)
    infections_by_age_dose = model |> flatten_chain .|> θ -> project_daily_infections_by_age_and_vaccine(θ, model, projection_date)

    uncertaintyinfections_by_age_unvac_first = [[infs[1][:, a, 1] for infs in infections_by_age_dose] for a = 1:size(infections_by_age_dose[1][1], 2)]
    uncertaintyinfections_by_age_unvac_reinf = [[infs[2][:, a, 1] .- infs[1][:, a, 1] for infs in infections_by_age_dose] for a = 1:size(infections_by_age_dose[1][1], 2)]
    uncertaintyinfections_by_age_1dose_first = [[infs[1][:, a, 2] for infs in infections_by_age_dose] for a = 1:size(infections_by_age_dose[1][1], 2)]
    uncertaintyinfections_by_age_1dose_reinf = [[infs[2][:, a, 2] .- infs[1][:, a, 2] for infs in infections_by_age_dose] for a = 1:size(infections_by_age_dose[1][1], 2)]
    uncertaintyinfections_by_age_2dose_first = [[infs[1][:, a, 3] for infs in infections_by_age_dose] for a = 1:size(infections_by_age_dose[1][1], 2)]
    uncertaintyinfections_by_age_2dose_reinf = [[infs[2][:, a, 3] .- infs[1][:, a, 3] for infs in infections_by_age_dose] for a = 1:size(infections_by_age_dose[1][1], 2)]

    n = size(uncertaintyinfections_by_age_unvac_first[1][1], 1)
    θs = model |> flatten_chain
    var_frqs = θs .|> θ -> get_variant_prediction(θ, model, n)
    wt_frqs = [data.prop_inc_wt for data in var_frqs]
    alpha_beta_frqs = [data.prop_inc_alpha_beta for data in var_frqs]
    delta_frqs = [data.prop_inc_delta for data in var_frqs]

    return (infections_by_age_unvac_first_pred = [get_credible_intervals(uncertaintyinfections_by_age_unvac_first[age]).pred for age = 1:6],
        infections_by_age_1dose_first_pred = [get_credible_intervals(uncertaintyinfections_by_age_1dose_first[age]).pred for age = 1:6],
        infections_by_age_2dose_first_pred = [get_credible_intervals(uncertaintyinfections_by_age_2dose_first[age]).pred for age = 1:6],
        infections_by_age_unvac_reinf_pred = [get_credible_intervals(uncertaintyinfections_by_age_unvac_reinf[age]).pred for age = 1:6],
        infections_by_age_1dose_reinf_pred = [get_credible_intervals(uncertaintyinfections_by_age_1dose_reinf[age]).pred for age = 1:6],
        infections_by_age_2dose_reinf_pred = [get_credible_intervals(uncertaintyinfections_by_age_2dose_reinf[age]).pred for age = 1:6],
        wt_frq_pred = get_credible_intervals(wt_frqs).pred,
        alpha_beta_frq_pred = get_credible_intervals(alpha_beta_frqs).pred,
        delta_frq_pred = get_credible_intervals(delta_frqs).pred)

end

"""
    function gather_post_mean_infection_rates_by_type_novac(model::KenyaCoVaccines.CoVAreaModel)

Gather the posterior mean predictions for infections by: age, primary/reinfection episode, variant. Only to be used with zero vaccination rates.
"""
function gather_post_mean_infection_rates_by_type_novac(model::KenyaCoVaccines.CoVAreaModel, projection_date::Date)
    firstinfections_by_age_dose = model |> flatten_chain .|> θ -> project_daily_infections_by_age_and_vaccine(θ, model, projection_date)

    uncertaintyfirstinfections_by_age_first = [[infs[1][:, a, 1] for infs in firstinfections_by_age_dose] for a = 1:size(firstinfections_by_age_dose[1][1], 2)]
    uncertaintyfirstinfections_by_age_reinf = [[infs[2][:, a, 1] .- infs[1][:, a, 1] for infs in firstinfections_by_age_dose] for a = 1:size(firstinfections_by_age_dose[1][1], 2)]

    n = size(uncertaintyfirstinfections_by_age_first[1][1], 1)
    θs = model |> flatten_chain
    var_frqs = θs .|> θ -> get_variant_prediction(θ, model, n)
    wt_frqs = [data.prop_inc_wt for data in var_frqs]
    alpha_beta_frqs = [data.prop_inc_alpha_beta for data in var_frqs]
    delta_frqs = [data.prop_inc_delta for data in var_frqs]

    return (infections_by_age_first_pred = [get_credible_intervals(uncertaintyfirstinfections_by_age_first[age]).pred for age = 1:6],
        infections_by_age_reinf_pred = [get_credible_intervals(uncertaintyfirstinfections_by_age_reinf[age]).pred for age = 1:6],
        wt_frq_pred = get_credible_intervals(wt_frqs).pred,
        alpha_beta_frq_pred = get_credible_intervals(alpha_beta_frqs).pred,
        delta_frq_pred = get_credible_intervals(delta_frqs).pred)

end

"""
    function pred_death_rate_county(mort_scale, mort_age, mort_var, mort_reinf, predictions, VE_severe_disease_risk, η, p_ID)

Given various parameters, convert the infection predictions into a death rate predictions.
"""
function pred_death_rate_county(mort_scale, mort_age, mort_var, mort_reinf, predictions, VE_severe_disease_risk, η, p_ID)
    pred_mean = vecvec_to_mat((1 - VE_severe_disease_risk[1]) .* mort_age .* predictions.infections_by_age_unvac_first_pred .+ (1 - VE_severe_disease_risk[2]) .* mort_age .* predictions.infections_by_age_1dose_first_pred .+ (1 - VE_severe_disease_risk[3]) .* mort_age .* predictions.infections_by_age_2dose_first_pred)
    pred_mean .+= mort_reinf * vecvec_to_mat((1 - VE_severe_disease_risk[1]) .* mort_age .* predictions.infections_by_age_unvac_reinf_pred .+ (1 - VE_severe_disease_risk[2]) .* mort_age .* predictions.infections_by_age_1dose_reinf_pred .+ (1 - VE_severe_disease_risk[3]) .* mort_age .* predictions.infections_by_age_2dose_reinf_pred)
    pred_mean .*= mort_scale .* (predictions.wt_frq_pred .+ mort_var[1] .* predictions.alpha_beta_frq_pred .+ mort_var[2] .* predictions.delta_frq_pred)'
    _p_ID = normalize(p_ID .* [exp(η * t) for t = 1:length(p_ID)], 1)
    for a = 1:6
        pred_mean[a, :] = [sum(pred_mean[a, max(1, t - 99):t] .* reverse(_p_ID[1:min(t, 100)])) for t = 1:size(pred_mean, 2)]
    end
    return pred_mean
end

"""
    function pred_death_rate_county_novac(mort_scale, mort_age, mort_var, mort_reinf, predictions, VE_severe_disease_risk, η, p_ID)

Given various parameters, convert the infection predictions into a death rate predictions. Only to be used with zero vaccination rates.
"""
function pred_death_rate_county_novac(mort_scale, mort_age, mort_var, mort_reinf, predictions, η, p_ID)
    pred_mean = vecvec_to_mat(mort_age .* predictions.infections_by_age_first_pred)
    pred_mean .+= mort_reinf * vecvec_to_mat(mort_age .* predictions.infections_by_age_reinf_pred)
    pred_mean .*= mort_scale .* (predictions.wt_frq_pred .+ mort_var[1] .* predictions.alpha_beta_frq_pred .+ mort_var[2] .* predictions.delta_frq_pred)'
    _p_ID = normalize(p_ID .* [exp(η * t) for t = 1:length(p_ID)], 1)
    for a = 1:6
        pred_mean[a, :] = [sum(pred_mean[a, max(1, t - 99):t] .* reverse(_p_ID[1:min(t, 100)])) for t = 1:size(pred_mean, 2)]
    end
    return pred_mean
end

"""
    function poisson_loss_death_pred_county(county_deaths, mort_scale, mort_age, mort_var, mort_reinf, predictions, η, p_ID)

Return the negative log-likelihood of the death data given the daily predictions.
"""
function poisson_loss_death_pred_county(county_deaths, mort_scale, mort_age, mort_var, mort_reinf, predictions, η, p_ID)
    preds = pred_death_rate_county_novac(mort_scale, mort_age, mort_var, mort_reinf, inf_predictions[1], η, p_ID)
    T = eltype(mort_scale)
    loss = T(0.0)
    for a = 1:size(county_deaths, 1), t = 31:size(county_deaths, 2)
        loss -= logpdf(Poisson(preds[a, t]), county_deaths[a, t])
    end
    return loss, sum(preds, dims = 1)[:]
end


text:
#!/usr/bin/env julia

const ITERS = 3e8
const Tlist = [0.9, 0.8, 0.7]#, 0.6, 0.5]
const Llist = [20]#[4, 6, 8, 12, 16, 20]

const param = "-m 100 -t 1e4 -T 1e8 -n 10 -L"

if isdir("jobs")
    error("Cowardly refusing to use existing directory 'jobs'")
else
    mkdir("jobs")
end

function write_job_tuples(;T=0.9, iters=3e8, L=8)
    globalid    = join(rand('a':'z',5))
    open("jobs/$(globalid).job", "w") do f
        println(f, "#!/bin/sh")
        println(f, "#\$ -o \$HOME/$(globalid).out -j y")
        println(f, "#\$ -N $(globalid)")
        println(f, "#\$ -pe mp16 16")
        println(f)
        println(f, "cd tfg/")
        println(f)
        for i in 1:16/2
            id    = join(rand('a':'z',5))
            meta  = "quench_T_$(T)_L_$(L)_$(id)"
            jseed = rand(UInt32)
            println(f, "./exe/annealer -j $jseed $param -p -l $L -i $iters $T \"$(meta)_RA.net\"")
            println(f, "./exe/annealer -j $jseed $param -p -l $L -i $iters $T \"$(meta)_RB.net\"")
        end
    end
end

f = 0
for T in Tlist, L in Llist
    write_job_tuples(iters=ITERS, T=T, L=L)
    f+=1
    print("$f files generated  \r")
end
println()

max_stars_repo_name: RedPointyJackson/tfg

max_issues_repo_name: RedPointyJackson/tfg

max_forks_repo_name: RedPointyJackson/tfg


text:
module LinQu

using ITensors, LinearAlgebra, Random

import Base: length,
			 copy,
			 push!,
			 isless,
			 delete!,
			 getindex,
			 setindex,
			 findfirst,
			 deleteat!,
			 print,
			 show,
			 replace!,
			 reverse,
			 size,
			 +,
			 -,
			 *,
			 /,
			 ==,
			 isapprox,
			 pop!,
			 range,
			 clamp
import LinearAlgebra: norm,
					  normalize!
import ITensors:getindex,
				 noprime,
				 noprime!,
				 prime,
				 svd,
				 qr,
				 position!,
				 setindex!,
				 commonindex,
				 iterate,
				 ITensor,
				 Index,
				 norm,
				 order,
				 *,
				 dag

# functions from Base, ITensors
export linkindex,
	   getindex,
	   noprime,
	   noprime!,
	   prime,
	   svd,
	   position!,
	   setindex!,
	   svd,
	   commonindex,
	   leftLim,
	   rightLim,
	   iterate,
	   ITensor,
	   Index,
	   randomITensor,
	   qr,
	   IndexSet,
	   order,
	   norm
# Types
export  QState,
		QGate,
		QCircuit,
		MPSState,
		ExactState,
		QGateBlock

include("extension.jl")
export getindex,
		contract,
		clamp

include("qgates/gatekernel.jl")
export  GateKernel,
		func,
		name
include("qgates/commongates/constants.jl")
export  IntFloat,
		Location,
		I_DATA,
		X_DATA,
		Y_DATA,
		Z_DATA,
		H_DATA,
		T_DATA,
		TDAG_DATA,
		S_DATA,
		SDAG_DATA,
		CNOT_DATA,
		SWAP_DATA,
		TOFFOLI_DATA,
		FREDKIN_DATA,
		CSWAP_DATA,
		CZ_DATA
include("qgates/commongates/commonkernel.jl")

include("qgates/qgate.jl")
export  QGate

include("qgates/constgate.jl")
export 	ConstGate,
		qubits,
		data,
		copy,
		control,
		randomConstGate

include("qgates/vargate.jl")
export  VarGate,
		qubits,
		param,
		func,
		copy,
		data,
		control,
		==,
		randomVarGate

include("qgates/measuregate.jl")
export MeasureGate,
	   qubits,
	   copy,
	   MEASURE_RESET,
	   MEASURE

include("qgates/commongates/commongates.jl")
export 	I,
		X,
		Y,
		Z,
		H,
		T,
		TDAG,
		S,
		SDAG,
		CNOT,
		SWAP,
		TOFFOLI,
		FREDKIN,
		CZ,
		Rx,
		Ry,
		Rz,
		Rϕ,
		ROTATE

include("qstates/qstate.jl")
export 	QState,
		probability

include("qstates/mpsstate/qubitsitemap.jl")
export 	QubitSiteMap,
		updateMap!

include("qstates/mpsstate/mpsstate.jl")
export MPSState,
	   swapSites!,
	   applyLocalGate!,
	   localizeQubits!,
	   centerAtSite!,
	   centerAtQubit!,
	   moveQubit!,
	   moveSite!,
	   orderQubits!,
	   apply!,
	   toExactState,
	   oneShot,
	   measure!,
	   collapse!,
	   showData,
	   normalize!,
	   dag,
	   ρ,
	   logprobability,
	   project!

include("qstates/exactstate.jl")
export 	ExactState,
		findindex,
		ITensor,
		toMPSState,
		isapprox

include("qgates/qgateblock.jl")
export  QGateBlock,
	    gates,
	    length,
	    push!,
	    add!,
	    addCopy!,
	    flatten,
	    extractParams,
	    insertParams!,
	    inverse,
	    randomQGateBlock

include("qgates/qcircuit.jl")
export	add!,
		QCircuit,
		flatten,
		randomQCircuit,
		showDetail

include("inter_qstate_qgate/exactstate_qgate.jl")
include("inter_qstate_qgate/mpsstate_qgate.jl")
include("inter_qstate_qgate/qstate_block_circuit.jl")
export  apply!

include("interface.jl")
export 	getindex,
		apply!

include("stat/sampling.jl")
export MCprobability!

# Following files are documentations
include("../doc/YQ_doc.jl")
end #module


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct VirtualMachineScaleSetVMInstanceView <: SwaggerModel
    platformUpdateDomain::Any # spec type: Union{ Nothing, Int32 } # spec name: platformUpdateDomain
    platformFaultDomain::Any # spec type: Union{ Nothing, Int32 } # spec name: platformFaultDomain
    rdpThumbPrint::Any # spec type: Union{ Nothing, String } # spec name: rdpThumbPrint
    vmAgent::Any # spec type: Union{ Nothing, VirtualMachineAgentInstanceView } # spec name: vmAgent
    maintenanceRedeployStatus::Any # spec type: Union{ Nothing, MaintenanceRedeployStatus } # spec name: maintenanceRedeployStatus
    disks::Any # spec type: Union{ Nothing, Vector{DiskInstanceView} } # spec name: disks
    extensions::Any # spec type: Union{ Nothing, Vector{VirtualMachineExtensionInstanceView} } # spec name: extensions
    vmHealth::Any # spec type: Union{ Nothing, VirtualMachineHealthStatus } # spec name: vmHealth
    bootDiagnostics::Any # spec type: Union{ Nothing, BootDiagnosticsInstanceView } # spec name: bootDiagnostics
    statuses::Any # spec type: Union{ Nothing, Vector{InstanceViewStatus} } # spec name: statuses
    assignedHost::Any # spec type: Union{ Nothing, String } # spec name: assignedHost
    placementGroupId::Any # spec type: Union{ Nothing, String } # spec name: placementGroupId

    function VirtualMachineScaleSetVMInstanceView(;platformUpdateDomain=nothing, platformFaultDomain=nothing, rdpThumbPrint=nothing, vmAgent=nothing, maintenanceRedeployStatus=nothing, disks=nothing, extensions=nothing, vmHealth=nothing, bootDiagnostics=nothing, statuses=nothing, assignedHost=nothing, placementGroupId=nothing)
        o = new()
        validate_property(VirtualMachineScaleSetVMInstanceView, Symbol("platformUpdateDomain"), platformUpdateDomain)
        setfield!(o, Symbol("platformUpdateDomain"), platformUpdateDomain)
        validate_property(VirtualMachineScaleSetVMInstanceView, Symbol("platformFaultDomain"), platformFaultDomain)
        setfield!(o, Symbol("platformFaultDomain"), platformFaultDomain)
        validate_property(VirtualMachineScaleSetVMInstanceView, Symbol("rdpThumbPrint"), rdpThumbPrint)
        setfield!(o, Symbol("rdpThumbPrint"), rdpThumbPrint)
        validate_property(VirtualMachineScaleSetVMInstanceView, Symbol("vmAgent"), vmAgent)
        setfield!(o, Symbol("vmAgent"), vmAgent)
        validate_property(VirtualMachineScaleSetVMInstanceView, Symbol("maintenanceRedeployStatus"), maintenanceRedeployStatus)
        setfield!(o, Symbol("maintenanceRedeployStatus"), maintenanceRedeployStatus)
        validate_property(VirtualMachineScaleSetVMInstanceView, Symbol("disks"), disks)
        setfield!(o, Symbol("disks"), disks)
        validate_property(VirtualMachineScaleSetVMInstanceView, Symbol("extensions"), extensions)
        setfield!(o, Symbol("extensions"), extensions)
        validate_property(VirtualMachineScaleSetVMInstanceView, Symbol("vmHealth"), vmHealth)
        setfield!(o, Symbol("vmHealth"), vmHealth)
        validate_property(VirtualMachineScaleSetVMInstanceView, Symbol("bootDiagnostics"), bootDiagnostics)
        setfield!(o, Symbol("bootDiagnostics"), bootDiagnostics)
        validate_property(VirtualMachineScaleSetVMInstanceView, Symbol("statuses"), statuses)
        setfield!(o, Symbol("statuses"), statuses)
        validate_property(VirtualMachineScaleSetVMInstanceView, Symbol("assignedHost"), assignedHost)
        setfield!(o, Symbol("assignedHost"), assignedHost)
        validate_property(VirtualMachineScaleSetVMInstanceView, Symbol("placementGroupId"), placementGroupId)
        setfield!(o, Symbol("placementGroupId"), placementGroupId)
        o
    end
end # type VirtualMachineScaleSetVMInstanceView

const _property_map_VirtualMachineScaleSetVMInstanceView = Dict{Symbol,Symbol}(Symbol("platformUpdateDomain")=>Symbol("platformUpdateDomain"), Symbol("platformFaultDomain")=>Symbol("platformFaultDomain"), Symbol("rdpThumbPrint")=>Symbol("rdpThumbPrint"), Symbol("vmAgent")=>Symbol("vmAgent"), Symbol("maintenanceRedeployStatus")=>Symbol("maintenanceRedeployStatus"), Symbol("disks")=>Symbol("disks"), Symbol("extensions")=>Symbol("extensions"), Symbol("vmHealth")=>Symbol("vmHealth"), Symbol("bootDiagnostics")=>Symbol("bootDiagnostics"), Symbol("statuses")=>Symbol("statuses"), Symbol("assignedHost")=>Symbol("assignedHost"), Symbol("placementGroupId")=>Symbol("placementGroupId"))
const _property_types_VirtualMachineScaleSetVMInstanceView = Dict{Symbol,String}(Symbol("platformUpdateDomain")=>"Int32", Symbol("platformFaultDomain")=>"Int32", Symbol("rdpThumbPrint")=>"String", Symbol("vmAgent")=>"VirtualMachineAgentInstanceView", Symbol("maintenanceRedeployStatus")=>"MaintenanceRedeployStatus", Symbol("disks")=>"Vector{DiskInstanceView}", Symbol("extensions")=>"Vector{VirtualMachineExtensionInstanceView}", Symbol("vmHealth")=>"VirtualMachineHealthStatus", Symbol("bootDiagnostics")=>"BootDiagnosticsInstanceView", Symbol("statuses")=>"Vector{InstanceViewStatus}", Symbol("assignedHost")=>"String", Symbol("placementGroupId")=>"String")
Base.propertynames(::Type{ VirtualMachineScaleSetVMInstanceView }) = collect(keys(_property_map_VirtualMachineScaleSetVMInstanceView))
Swagger.property_type(::Type{ VirtualMachineScaleSetVMInstanceView }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_VirtualMachineScaleSetVMInstanceView[name]))}
Swagger.field_name(::Type{ VirtualMachineScaleSetVMInstanceView }, property_name::Symbol) =  _property_map_VirtualMachineScaleSetVMInstanceView[property_name]

function check_required(o::VirtualMachineScaleSetVMInstanceView)
    true
end

function validate_property(::Type{ VirtualMachineScaleSetVMInstanceView }, name::Symbol, val)
end


text:
function truncate!(λ, D)
    # truncate variational standard deviation
    for d in 1:D
        λ[d, 2] = λ[d, 2] < invsoftplus(1e-5) ?
                  invsoftplus(1e-5) : λ[d, 2]
        λ[d, 2] = λ[d, 2] > softplus(log(prevfloat(Inf))) ?
                  invsoftplus(log(prevfloat(Inf))) : λ[d, 2]
    end
end

function sigmoid(z)
    return 1.0 / (1.0 + exp(-z))
end