text:
#!/usr/bin/env julia

const OUTPUT = haskey(ENV, "OUTPUT")

immutable Node
    starting::Int32
    no_of_edges::Int32
end

function Usage()
    prog = basename(Base.source_path())
    println(STDERR,"Usage: ", prog, " <input_file>")
    exit(1)
end

function parseline{T<:Number}(::Type{T}, f)
    while true
        line = chomp(readline(f))
        if length(line) > 0
            return map(s -> parse(T, s), split(line))::Vector{T}
        end
    end
    return T[]
end

function main(args)
    if length(args) != 1
        Usage()
    end

    info("Reading file")
    f = open(args[1],"r")

    # Read graph data
    no_of_nodes, = parseline(Int,f)
    h_graph_nodes = Array{Node, 1}(no_of_nodes)
    for i in 1:no_of_nodes
        start, edgeno = parseline(Int, f)
        # add one to every node index to accomodate for 1-based arrays
        h_graph_nodes[i] = Node(start+1, edgeno)
    end

    # initialize node structures
    # the type specification should not be necessary, but is at least on 0.4.5
    h_graph_mask::SharedArray{Bool,1} = SharedArray(Bool,(no_of_nodes,),init = h_graph_mask -> h_graph_mask[Base.localindexes(h_graph_mask)] = false)
    h_updating_graph_mask::SharedArray{Bool,1} = SharedArray(Bool,(no_of_nodes,),init = h_updating_graph_mask -> h_updating_graph_mask[Base.localindexes(h_updating_graph_mask)] = false)
    h_graph_visited = zeros(Bool,no_of_nodes)

    # read the source node from the file (+1 to accomodate 1-based arrays)
    source, = parseline(Int,f)
    source += 1

    # set the source node as true in the mask
    h_graph_mask[source]=true
    h_graph_visited[source]=true

    # skip empty line
    edge_list_size, = parseline(Int,f)
    h_graph_edges_with_cost::Array{Int32,2} = readdlm(f,' ',UInt32; use_mmap=true, dims=(edge_list_size,2), quotes=false, comments=false)
    # add one to every node index to accomodate for 1-based arrays
    h_graph_edges = [h_graph_edges_with_cost[row,1]+1 for row in 1:size(h_graph_edges_with_cost,1)]

    close(f)

    h_cost::SharedArray{Int32,1} = SharedArray(Int32,(no_of_nodes,), init = h_cost -> h_cost[Base.localindexes(h_cost)] = -1)
    h_cost[source] = 0

    info("Start traversing the tree")

    tic()
    again = true
    while again
        # if no one changes this value then the loop stops
        again = false

        @sync @parallel for maskiter in eachindex(h_graph_mask)
            if h_graph_mask[maskiter]
                h_graph_mask[maskiter] = false
                for nodeiter in h_graph_nodes[maskiter].starting:h_graph_nodes[maskiter].starting+h_graph_nodes[maskiter].no_of_edges-1
                    id = h_graph_edges[nodeiter]
                    if !h_graph_visited[id]
                        h_cost[id] = h_cost[maskiter]+1
                        h_updating_graph_mask[id]=true
                    end
                end
            end
        end

        for updatingmaskiter in eachindex(h_updating_graph_mask)
            if h_updating_graph_mask[updatingmaskiter]
                h_graph_mask[updatingmaskiter]=true
                h_graph_visited[updatingmaskiter]=true
                again=true
                h_updating_graph_mask[updatingmaskiter]=false
            end
        end
    end
    toc()

    if OUTPUT
        f = open("output.txt","w")
        for i in eachindex(h_cost)
            println(f,i-1,") cost:",h_cost[i])
        end
        close(f)
        info("Result stored in output.txt")
    end
end


text:
export OutcomeSamplingMCCFRPolicy

using Random
using StatsBase: sample, Weights

"""
    OutcomeSamplingMCCFRPolicy

This implementation uses stochasticaly-weighted averaging.

Ref:

- [MONTE CARLO SAMPLING AND REGRET MINIMIZATION FOR EQUILIBRIUM COMPUTATION AND DECISION-MAKING IN LARGE EXTENSIVE FORM GAMES](http://mlanctot.info/files/papers/PhD_Thesis_MarcLanctot.pdf)
- [Monte Carlo Sampling for Regret Minimization in Extensive Games](https://papers.nips.cc/paper/3713-monte-carlo-sampling-for-regret-minimization-in-extensive-games.pdf)
"""
struct OutcomeSamplingMCCFRPolicy{S,T,R<:AbstractRNG} <: AbstractPolicy
    nodes::Dict{S,InfoStateNode}
    behavior_policy::QBasedPolicy{TabularLearner{S,T},WeightedExplorer{true,R}}
end

(p::OutcomeSamplingMCCFRPolicy)(env::AbstractEnv) = p.behavior_policy(env)

RLBase.get_prob(p::OutcomeSamplingMCCFRPolicy, env::AbstractEnv) =
    get_prob(p.behavior_policy, env)

function OutcomeSamplingMCCFRPolicy(;
    env::AbstractEnv,
    n_iter::Int,
    rng = Random.GLOBAL_RNG,
    ϵ = 0.6,
)
    @assert NumAgentStyle(env) isa MultiAgent
    @assert DynamicStyle(env) === SEQUENTIAL
    @assert RewardStyle(env) === TERMINAL_REWARD
    @assert ChanceStyle(env) === EXPLICIT_STOCHASTIC
    @assert DefaultStateStyle(env) === Information{String}()

    nodes = init_info_state_nodes(env)

    for i in 1:n_iter
        for p in get_players(env)
            if p != get_chance_player(env)
                outcome_sampling(copy(env), p, nodes, ϵ, 1.0, 1.0, 1.0, rng)
            end
        end
    end

    behavior_policy = QBasedPolicy(;
        learner = TabularLearner{String}(),
        explorer = WeightedExplorer(; is_normalized = true, rng = rng),
    )

    for (k, v) in nodes
        s = sum(v.cumulative_strategy)
        if s != 0
            update!(behavior_policy, k => v.cumulative_strategy ./ s)
        end
    end

    OutcomeSamplingMCCFRPolicy(nodes, behavior_policy)
end

function outcome_sampling(env, i, nodes, ϵ, πᵢ, π₋ᵢ, s, rng)
    current_player = get_current_player(env)

    if get_terminal(env)
        get_reward(env, i) / s, 1.0
    elseif current_player == get_chance_player(env)
        env(rand(rng, get_actions(env)))
        outcome_sampling(env, i, nodes, ϵ, πᵢ, π₋ᵢ, s, rng)
    else
        I = get_state(env)
        node = nodes[I]
        regret_matching!(node)
        σ, rI, sI = node.strategy, node.cumulative_regret, node.cumulative_strategy
        n = length(node.strategy)

        if i == current_player
            aᵢ = rand(rng) >= ϵ ? sample(rng, Weights(σ, 1.0)) : rand(rng, 1:n)
            pᵢ = σ[aᵢ] * (1 - ϵ) + ϵ / n
            πᵢ′, π₋ᵢ′, s′ = πᵢ * pᵢ, π₋ᵢ, s * pᵢ
        else
            aᵢ = sample(rng, Weights(σ, 1.0))
            pᵢ = σ[aᵢ]
            πᵢ′, π₋ᵢ′, s′ = πᵢ, π₋ᵢ * pᵢ, s * pᵢ
        end

        env(get_legal_actions(env)[aᵢ])
        u, πₜₐᵢₗ = outcome_sampling(env, i, nodes, ϵ, πᵢ′, π₋ᵢ′, s′, rng)

        if i == current_player
            w = u * π₋ᵢ
            rI .+= w * πₜₐᵢₗ .* ((1:n .== aᵢ) .- σ[aᵢ])
        else
            sI .+= π₋ᵢ / s .* σ
        end

        u, πₜₐᵢₗ * σ[aᵢ]
    end
end

ext: jl


text:
module HW2

using Obfuscatee
using Nettle
using JSON
using ProgressMeter
using Random

using POMDPs
using StaticArrays
using Distributions
using Discretizers
using POMDPModelTools
using POMDPPolicies
using POMDPModels: SimpleGridWorld

import UnicodePlots
using Compose: compose, context, text, polygon, hcenter, hleft, hright, vtop

import DMUStudent

export
    UnresponsiveACASMDP,
    ACASState,
    transition_matrices,
    reward_vectors,
    grid_world,
    render

# Problem 4: Grid World

const grid_world = SimpleGridWorld()

# Problem 5: ACAS

const v = 44000 # 500 mph in ft per minute
const hlim = (27500, 32500) # ft
const hdotlim = (-3750, 3750)
const h_i_stdev = 500 # ft per sqrt(minute)
const base_dbins = 10
const base_hbins = 5
const nhdotbins = 5
const delta_hdot = 1500 # ft/min
const timespan = 150/60 # minutes
const maxd = timespan*2*v

struct ACASState <: FieldVector{4, Float64}
    h_o::Float64
    hdot_o::Float64
    h_i::Float64
    d::Float64
end

struct UnresponsiveACASMDP <: MDP{ACASState, Int}
    dt::Float64 # minutes
    hbins::LinearDiscretizer{Float64, Int}
    hdotbins::LinearDiscretizer{Float64, Int}
    dbins::LinearDiscretizer{Float64, Int}
    _cached_cdf::Vector{Float64}
end

function UnresponsiveACASMDP(n)
    nhbins = base_hbins*n
    hbins = LinearDiscretizer(range(hlim[1], hlim[2], length=nhbins+1))
    @assert nlabels(hbins) == nhbins

    ndbins = base_dbins*n
    dbins = LinearDiscretizer(range(0.0, maxd, length=ndbins+1))
    @assert nlabels(dbins) == ndbins

    hdotbins = LinearDiscretizer(range(hdotlim[1], hdotlim[2], length=nhdotbins+1))

    dt = timespan/ndbins
    ccdf = cached_cdf(hbins, sqrt(dt*h_i_stdev^2))
    return UnresponsiveACASMDP(dt, hbins, hdotbins, dbins, ccdf)
end

POMDPs.actions(m::UnresponsiveACASMDP) = (-delta_hdot, 0, delta_hdot)
POMDPs.states(m::UnresponsiveACASMDP) = [convert_s(ACASState, si, m) for si in 1:prod(map(nlabels, (m.hbins, m.hdotbins, m.hbins, m.dbins)))]

POMDPs.discount(::UnresponsiveACASMDP) = 0.99

function POMDPs.reward(m::UnresponsiveACASMDP, s::ACASState, a, sp::ACASState)
    r = 0.0
    if is_nmac(m, sp)
        r -= 100.0
    end
    r -= abs(sign(a))
    return r
end
POMDPs.reward(m::UnresponsiveACASMDP, v::AbstractVector, a, vp::AbstractVector) = reward(m, convert(ACASState, v), a, convert(ACASState, vp))

function POMDPs.transition(m::UnresponsiveACASMDP, s::ACASState, a)
    d = s.d - 2*v*m.dt
    hdot_o = clamp(s.hdot_o + a, -3000, 3000)
    h_o = bincenter(m.hbins, encode(m.hbins, s.h_o + hdot_o*m.dt))

    nhbins = nlabels(m.hbins)
    h_i_bin = encode(m.hbins, s.h_i)

    states = Vector{ACASState}(undef, nhbins)
    for (i, h) in enumerate(bincenters(m.hbins))
        states[i] = ACASState(h_o, hdot_o, h, d)
    end

    start = nhbins-h_i_bin + 1
    @assert 1 <= start <= length(m._cached_cdf)
    finish = start + nhbins - 2
    @assert finish <= length(m._cached_cdf)
    probs = Vector{Float64}(undef, nhbins)
    probs[1] = m._cached_cdf[start]
    probs[2:nhbins-1] = view(m._cached_cdf, start+1:finish) - view(m._cached_cdf, start:finish-1)
    probs[end] = 1.0 - m._cached_cdf[finish]
    @assert length(probs) == nhbins
    @assert length(states) == nhbins
    @assert sum(probs) ≈ 1

    return SparseCat(states, probs)
end
POMDPs.transition(m::UnresponsiveACASMDP, s::AbstractVector, a) = transition(m::UnresponsiveACASMDP, convert(ACASState, s), a)

function POMDPs.isterminal(m::UnresponsiveACASMDP, s::ACASState)
    if is_nmac(m, s) || s.d <= binedges(m.dbins)[2]
        return true
    else
        return false
    end
end
POMDPs.isterminal(m::UnresponsiveACASMDP, s::AbstractVector) = isterminal(m, convert(ACASState, s))

POMDPs.initialstate(::UnresponsiveACASMDP) = Deterministic(ACASState(sum(hlim)/2, 0.0, sum(hlim)/2, maxd))

function POMDPs.convert_s(::Type{ACASState}, s::Integer, m::UnresponsiveACASMDP)
    shape = map(nlabels, (m.hbins, m.hdotbins, m.hbins, m.dbins))
    bins = CartesianIndices(shape)[s]
    h_o = bincenter(m.hbins, bins[1])
    hdot_o = bincenter(m.hdotbins, bins[2])
    h_i = bincenter(m.hbins, bins[3])
    d = maxd-bincenter(m.dbins, bins[4])
    return ACASState(h_o, hdot_o, h_i, d)
end

function POMDPs.convert_s(::Type{Int}, s::ACASState, m::UnresponsiveACASMDP)
    return stateindex(m, s)
end
POMDPs.convert_s(::Type{Int}, s::AbstractVector, m::UnresponsiveACASMDP) = convert_s(Int, convert(ACASState, s), m)

function POMDPs.stateindex(m::UnresponsiveACASMDP, v)
    s = convert(ACASState, v)
    h_o_bin = encode(m.hbins, s.h_o)
    hdot_o_bin = encode(m.hdotbins, s.hdot_o)
    h_i_bin = encode(m.hbins, s.h_i)
    dbin = encode(m.dbins, maxd-s.d)
    shape = map(nlabels, (m.hbins, m.hdotbins, m.hbins, m.dbins))
    return LinearIndices(shape)[h_o_bin, hdot_o_bin, h_i_bin, dbin]
end

POMDPs.actionindex(m::UnresponsiveACASMDP, a) = sign(a)+2

# less than 100 ft vertically and 500 ft horizontally
function is_nmac(m::UnresponsiveACASMDP, s::ACASState)
    return (s.d < binedges(m.dbins)[2] || s.d <= 500) && abs(s.h_o - s.h_i) <= 100
end

function cached_cdf(hbins, stdev)
    nbins = nlabels(hbins)
    binsize = (hlim[2]-hlim[1])/nbins
    @assert binsize ≈ mean(binwidths(hbins))
    d = Normal(0.0, stdev)
    npoints = 2*nbins-1
    lower = -(nbins - 0.5)*binsize

    cdfvals = Vector{Float64}(undef, npoints)
    for p in 1:npoints
        pt = lower + p*binsize
        cdfvals[p] = cdf(d, pt)
    end

    # some sanity checks
    @assert lower ≈ hlim[1] - (hlim[2] - 0.5*binsize)
    @assert lower + binsize*npoints ≈ hlim[2] - (hlim[1] + 0.5*binsize)
    @assert all(0 .<= cdfvals .<= 1)

    return cdfvals
end

function bincenter(d::LinearDiscretizer, binindex)
    be = binedges(d)
    return (be[binindex+1] + be[binindex])/2
end

"""
    transition_matrices(m, [sparse=true])

Create a dictionary mapping actions to transition matrices for MDP m.

# Example
```julia
T = transition_matrices(m)
T[1][2,3] # probability of transitioning from state 2 to 3 when action 1 is taken
```

# Arguments
- `m`: `MDP` model

# Keyword Arguments
- `sparse::Bool`: if true, returns a sparse matrix representation (with significant memory savings!), if false, the matrices will be dense
"""
function transition_matrices(m::MDP; sparse::Bool=false)
    transmats = POMDPModelTools.transition_matrix_a_s_sp(m)
    if !sparse
        transmats = [convert(Matrix, t) for t in transmats]
    end
    mtype = typeof(first(transmats))
    oa = ordered_actions(m)
    return Dict{actiontype(m), mtype}(oa[ai]=>transmats[ai] for ai in 1:length(actions(m)))
end

function reward_vectors(m::MDP)
    d = Dict{actiontype(m), Vector{Float64}}()
    for a in actions(m)
        rv = POMDPModelTools.policy_reward_vector(m, FunctionPolicy(s->a))
        d[a] = rv
    end
    return d
end

include("hw2_vis.jl")

@binclude(".bin/hw2_eval")

end


text:
require("DataFrames")
using DataFrames

df = read_table("demo/senate112.csv")
senator_names = df[:, 425]

dm = DataArray(df[:, 1:424])

u, d, v = svd(dm, 2)

ideal_points = u * diagm(d)

df = DataFrame()
df["Senator"] = senator_names
df["X"] = DataArray(ideal_points[:, 1])
df["Y"] = DataArray(ideal_points[:, 2])

write_table("demo/ideal_points.tsv", df)


text:
@testset "Expression destructuring"  begin

    @testset "expr_bindings" begin

        @test expr_bindings(:(x = 5), :(x = _val_)) == Dict(:_val_ => 5)
        @test expr_bindings(:(x = 5), :(_sym_ = 5)) == Dict(:_sym_ => :x)
        @test expr_bindings(:(x = 5), :_ex_) == Dict(:_ex_ => :(x = 5))

        splatex = Expr(:let, :_body_, :_SPLAT_bindings_)
        ex = :(let x=1, y=2, z=3
                   x + y
                   y + z
               end)

        bindings = expr_bindings(ex, splatex)
        body = bindings[:_body_]
        b = bindings[:_bindings_]

        @test haskey(bindings, :_body_)
        @test haskey(bindings, :_bindings_)
        @test isa(body, Expr) && body.head == :block
        @test isa(b, Array) && length(b) == 3
    end

    @testset "expr_replace"  begin

        ex = quote
            let x, y, z
                bar
                x + y
                y + z
            end
        end

        template = quote
            let _SPLAT_bindings_
                _funname_
                _SPLAT_body_
            end
        end

        out = quote
            function _funname_( _UNSPLAT_bindings_)
                _UNSPLAT_body_
            end
        end

        foofun = expr_replace(ex, template, out)
        eval(foofun)

        @test foofun.head == :function
        @test eval(foofun)(1,2,3) == 5
        @test bar(1,2,3) == 5
    end

end


text:
using Test, YaoBlocks, YaoBase

rp = RepeatedBlock{5}(X, (1, 2, 3))
@test isreflexive(rp)
@test ishermitian(rp)
@test isunitary(rp)
@test (chsubblocks(rp, [Z]) |> subblocks .== [Z]) |> all
@test occupied_locs(rp) == (1, 2, 3)
@test rp |> copy == rp
@test YaoBlocks.PropertyTrait(rp) == YaoBlocks.PreserveAll()

@test repeat(10, H, 1:10) == repeat(10, H, Tuple(1:10))


text:
@testset "positional types" begin
    dct = Dict("foo" => 37, "bar" => 23)
    p = @inferred Plot barplot(dct)
    test_ref("references/barplot/default.txt", @print_col(p))
    test_ref("references/barplot/default.txt", @show_col(p))
    test_ref("references/barplot/nocolor.txt", @show_nocol(p))

    dct = Dict("foo" => 37., :bar => 23., 2.1 => 10)
    p = @inferred Plot barplot(dct)
    test_ref("references/barplot/default_mixed.txt", @print_col(p))

    dct = Dict(:foo => 37, :bar => 23)
    p = @inferred Plot barplot(dct)
    test_ref("references/barplot/default.txt", @print_col(p))

    p = @inferred barplot([:bar, :foo], Int16[23, 37])
    test_ref("references/barplot/default.txt", @print_col(p))

    p = @inferred barplot(["bar", "foo"], [23, 37])
    test_ref("references/barplot/default.txt", @print_col(p))
    @test_throws DimensionMismatch barplot!(p, ["zoom"], [90, 80])
    @test_throws DimensionMismatch barplot!(p, ["zoom", "boom"], [90])
    @test_throws MethodError barplot!(p, ["zoom"], [90.])
    @test_throws InexactError barplot!(p, "zoom", 90.1)
    @test @inferred(barplot!(p, ["zoom"], [90])) === p
    test_ref("references/barplot/default2.txt", @print_col(p))

    p = @inferred barplot(["bar", "foo"], [23, 37])
    @test @inferred(barplot!(p, :zoom, 90.)) === p
    test_ref("references/barplot/default2.txt", @print_col(p))

    p = @inferred barplot(["bar", "foo"], [23, 37])
    @test_throws MethodError barplot!(p, Dict("zoom" => 90.))
    @test @inferred(barplot!(p, Dict("zoom" => 90))) === p
    test_ref("references/barplot/default2.txt", @print_col(p))

    p = @inferred barplot(["bar", "foo"], [23, 37])
    @test @inferred(barplot!(p, Dict(:zoom => 90))) === p
    test_ref("references/barplot/default2.txt", @print_col(p))

    p = @inferred barplot(2:6, 11:15)
    test_ref("references/barplot/ranges.txt", @print_col(p))
    @test @inferred(barplot!(p, 9:10, 20:21)) === p
    test_ref("references/barplot/ranges2.txt", @print_col(p))
end

struct Scale{T} r::T end
(f::Scale)(x) = f.r * x  # functor

@testset "keyword arguments" begin
    p = @inferred barplot(
        [:a, :b, :c, :d, :e],
        [0, 1, 10, 100, 1000],
        title = "Logscale Plot",
        xscale = :log10,
    )
    test_ref("references/barplot/log10.txt", @print_col(p))

    p = @inferred barplot(
        [:a, :b, :c, :d, :e],
        [0, 1, 10, 100, 1000],
        title = "Logscale Plot",
        xscale = Scale(π),
    )
    test_ref("references/barplot/functor.txt", @print_col(p))

    p = @inferred barplot(
        [:a, :b, :c, :d, :e],
        [0, 1, 10, 100, 1000],
        title = "Logscale Plot",
        xlabel = "custom label",
        xscale = :log10,
    )
    test_ref("references/barplot/log10_label.txt", @print_col(p))

    p = @inferred barplot(
        ["Paris", "New York", "Moskau", "Madrid"],
        [2.244, 8.406, 11.92, 3.165],
        title = "Relative sizes of cities",
        xlabel = "population [in mil]",
        color = :blue,
        margin = 7,
        padding = 3,
    )
    test_ref("references/barplot/parameters1.txt", @print_col(p))

    p = @inferred barplot(
        ["Paris", "New York", "Moskau", "Madrid"],
        [2.244, 8.406, 11.92, 3.165],
        title = "Relative sizes of cities",
        xlabel = "population [in mil]",
        color = (0,0,255),
        margin = 7,
        padding = 3,
        labels = false,
    )
    test_ref("references/barplot/parameters1_nolabels.txt", @print_col(p))

    p = @inferred barplot(
        ["Paris", "New York", "Moskau", "Madrid"],
        [2.244, 8.406, 11.92, 3.165],
        title = "Relative sizes of cities",
        xlabel = "population [in mil]",
        color = :yellow,
        border = :solid,
        symbols = ["="],
        width = 60
    )
    test_ref("references/barplot/parameters2.txt", @print_col(p))

    # same but with Char as symbols
    p = @inferred barplot(
        ["Paris", "New York", "Moskau", "Madrid"],
        [2.244, 8.406, 11.92, 3.165],
        title = "Relative sizes of cities",
        xlabel = "population [in mil]",
        color = :yellow,
        border = :solid,
        symbols = ['='],
        width = 60
    )
    test_ref("references/barplot/parameters2.txt", @print_col(p))
end

@testset "edge cases" begin
    @test_throws ArgumentError barplot([:a, :b], [-1, 2])
    p = barplot([5,4,3,2,1], [0,0,0,0,0])
    test_ref("references/barplot/edgecase_zeros.txt", @print_col(p))
    p = barplot([:a,:b,:c,:d], [1,1,1,1000000])
    test_ref("references/barplot/edgecase_onelarge.txt", @print_col(p))
    barplot("one", 1)
end

@testset "colors" begin
    p = barplot(["B","A"], [2,1], color=9)
    test_ref("references/barplot/col1.txt", @print_col(p))
    p = barplot(["B","A"], [2,1], color=(200,50,0))
    test_ref("references/barplot/col2.txt", @print_col(p))
end


text:
"""
Loads in PowerModels network data given the name of a network case file,
then starts creating samples with distributed processing
"""
function dist_create_samples(net::String, K=Inf; U=0.0, S=0.0, V=0.0, max_iter=Inf, T=Inf, discard=false, variance=false,
							input_vars=DEFAULT_INPUTS, output_vars=DEFAULT_OUTPUTS, dual_vars=DEFAULT_DUALS,
							sampler=sample_polytope_cprnd, sampler_opts=Dict{Symbol,Any}()::Dict{Symbol}, A=[]::Array, b=[]::Array,
							pd_max=nothing, pd_min=nothing, pf_min=0.7071, pf_lagging=true, reset_level=0, nproc=nothing, replace_samples=false,
							save_max_load=false, save_certs=false,
							print_level=0, stat_track=false, save_while=false, save_infeasible=false, save_path="", net_path="",
							model_type=PM.QCLSPowerModel, r_solver=JuMP.optimizer_with_attributes(Ipopt.Optimizer, "tol" => TOL),
							opf_solver=JuMP.optimizer_with_attributes(Ipopt.Optimizer, "tol" => TOL))
	net, net_path = load_net(net, net_path, print_level)

	return dist_create_samples(net, K::Integer; U=U, S=S, V=V, max_iter=max_iter, T=T, discard=discard, variance=variance,
							  input_vars=input_vars, output_vars=output_vars, dual_vars=dual_vars,
							  sampler=sampler, sampler_opts=sampler_opts, A=A, b=b, reset_level=reset_level,
							  pd_max=pd_max, pd_min=pd_min, pf_min=pf_min, pf_lagging=pf_lagging, nproc=nproc, replace_samples=replace_samples,
							  save_certs=save_certs, save_max_load=save_max_load,
							  print_level=print_level, stat_track=stat_track, save_while=save_while,
							  save_infeasible=save_infeasible, save_path=save_path, net_path=net_path,
							  model_type=model_type, r_solver=r_solver, opf_solver=opf_solver)
end

"""
Creates an AC OPF dataset for the given PowerModels network dictionary. Generates samples until one of the given stopping criteria is met.
Takes options to determine how to sample points, what information to save, and what information is printed.

# Arguments
- 'net::Dict': network information stored in a PowerModels.jl format specified dictionary
- 'K::Integer': the maximum number of samples before stopping sampling
- 'U::Float': the minimum % of unique active sets sampled in the previous 1 / U samples to continue sampling
- 'S::Float': the minimum % of saved samples in the previous 1 / L samples to continue sampling
- 'V::Float': the minimum % of feasible samples that increase the variance of the dataset in the previous 1 / L samples to continue sampling
- 'T::Integer': the maximum time for the sampler to run in seconds.
- 'max_iter::Integer': maximum number of iterations for the sampler to run for.
- 'nproc::Integer': the number of processors for the sampler to run with. Defaults to the number reported by Distributed.nprocs().
- 'replace_samples::Bool': whether samples in the samples channel are replaced when a new infeasibility certificate is added. Found to sometimes block progress when turned on.
- 'sampler::Function': the sampling function to use. This function must take arguements A and b, and can take optional arguments.
- 'sampler_opts::Dict': a dictionary of optional arguments to pass to the sampler function.
- 'A::Array': defines the initial sampling space polytope Ax<=b. If not provided, initializes to a default.
- 'b::Array': defines the initial sampling space polytope Ax<=b. If not provided, initializes to a default.
- 'pd_max::Array': the maximum active load values to use when initializing the sampling space and constraining the loads. If nothing, finds the maximum load at each bus with the given relaxed model type.
- 'pd_min::Array': the minimum active load values to use when initializing the sampling space and constraining the loads. If nothing, this is set to 0 for all loads.
- 'pf_min::Array/Float:' the minimum power factor for all loads in the system (Number) or an array of minimum power factors for each load in the system.
- 'pf_lagging::Bool': indicating if load power factors can be only lagging (True), or both lagging or leading (False).
- 'reset_level::Integer': determines how to reset the load point to be inside the polytope before sampling. 2: Reset closer to nominal load & chebyshev center, 1: Reset closer to chebyshev center, 0: Reset at chebyshev center.
- 'save_certs::Bool': specifies whether the sampling space, Ax<=b (A & b matrices) are saved to the results dictionary.
- 'save\\_max_load::Bool': specifies whether the max active load demands used are saved to the results dictionary.
- 'model_type::Type': an abstract PowerModels type indicating the network model to use for the relaxed AC-OPF formulations (Max Load & Nearest Feasible)
- 'r_solver': an optimizer constructor used for solving the relaxed AC-OPF optimization problems.
- 'opf_solver': an optimizer constructor used to find the AC-OPF optimal solution for each sample.
- 'print_level::Integer': from 0 to 3 indicating the level of info to print to console, with 0 indicating minimum printing.
- 'stat_track::Integer': from 0 to 3 indicating the level of stats info saved during each iteration	0: No information saved, 1: Feasibility, New Certificate, Added Sample, Iteration Time, 2: Variance for all input & output variables
- 'save_while::Bool': indicates whether results and stats information is saved to a csv file during processing.
- 'save_infeasible::Bool': indicates if infeasible samples are saved. If true saves infeasible samples in a seperate file from feasible samples.
- 'save_path::String:' a string with the file path to the desired result save location.
- 'net_path::String': a string with the file path to the network file.
- 'variance::Bool': indicates if dataset variance information is tracked for each unique active set.
- 'discard::Bool': indicates if samples that do not increase the variance within a unique active set are discarded.

See 'OPF-Learn: An Open-Source Framework for Creating Representative AC Optimal Power Flow Datasets'
for more information on how the AC OPF datasets are created.

Modified from AgenerateACOPFsamples.m written by Ahmed S. Zamzam
"""
function dist_create_samples(net::Dict, K=Inf; U=0.0, S=0.0, V=0.0, max_iter=Inf, T=Inf, nproc=nothing, discard=false, variance=false,
								input_vars=DEFAULT_INPUTS, output_vars=DEFAULT_OUTPUTS, dual_vars=DEFAULT_DUALS,
								sampler=sample_polytope_cprnd, sampler_opts=Dict{Symbol,Any}()::Dict{Symbol}, A=[]::Array, b=[]::Array,
								pd_max=nothing, pd_min=nothing, pf_min=0.7071, pf_lagging=true, reset_level=0, replace_samples=false,
								save_max_load=false, save_certs=false,
								print_level=0, stat_track=false, save_while=false, save_infeasible=false, save_path="", net_path="",
								model_type=PM.QCLSPowerModel, r_solver=JuMP.optimizer_with_attributes(Ipopt.Optimizer, "tol" => TOL),
								opf_solver=JuMP.optimizer_with_attributes(Ipopt.Optimizer, "tol" => TOL))
	# Create channels for transfering data between processes
	isnothing(nproc) && (nproc = Distributed.nprocs())
	@assert nproc > 3 "Not enough processors available, nprocs:$(nproc). Need 4+ CPUs for improved runtimes."
	@assert nproc <= Distributed.nprocs() "Number of distributed processors added, $(Distributed.nprocs()), must be greater than specified nproc, $(nproc)"
	pid = Distributed.myid()
	sample_chnl = Distributed.Channel{Any}(4 * nproc)
	polytope_chnl = Distributed.Channel{Any}(4 * nproc)
	result_chnl = Distributed.Channel{Any}(4 * nproc)
	final_chnl = Distributed.Channel{Any}(1)
	sample_ch = Distributed.RemoteChannel(()->sample_chnl, pid)
	polytope_ch = Distributed.RemoteChannel(()->polytope_chnl, pid)
	result_ch = Distributed.RemoteChannel(()->result_chnl, pid)
	final_ch = Distributed.RemoteChannel(()->final_chnl, pid)

	num_procs = nproc - 2  #TASK: Determine why the producer gets stuck running on the main proc

	net_name = net["name"]
	save_order = vcat(input_vars, output_vars, dual_vars)
    if net_path == ""  # Set net path to save path if not specified for saving network max loads
		net_path = save_path
	end

	# Gather network information used during processing
	A, b, x, results, fnfp_model, base_load_feasible, net_r = initialize(net, pf_min, pf_lagging, pd_max, pd_min,
																		 input_vars, output_vars, dual_vars,
																		 save_certs, save_infeasible, save_while,
																		 stat_track, save_max_load, A, b,
																		 sampler, sampler_opts, net_name,
																		 net_path, model_type, r_solver,
																		 reset_level, print_level)

	(print_level > 0) && println("Starting sampling...")
	procs = Distributed.workers()
	producer = Distributed.remotecall(sample_producer, procs[1],
						  A, b, sampler, sampler_opts, base_load_feasible, K, U, S, V, max_iter, T, num_procs,
						  results, discard, variance, reset_level, save_while, save_infeasible,
						  stat_track, save_certs, net_name, dual_vars, save_order, replace_samples,
						  save_path, sample_ch, polytope_ch, result_ch, final_ch, print_level)
	for proc in procs[2:end]
		a = Distributed.remotecall(sample_processor, proc, net, net_r, r_solver, opf_solver,
							   sample_ch, final_ch, polytope_ch, result_ch,
							   print_level, model_type)
	end

	results = Distributed.fetch(producer)
	# Need to convert the unique active sets Set to and Array, due to a bug with PyJulia
	results["duals"]["unique_active_sets"] = collect(results["duals"]["unique_active_sets"])
	return results
end

"UNUSED"
function initialize_distributed(net, num_procs; pd_max=nothing, save_path="",  net_path="", save_while=false, A=[], b=[],
							  line_constraints=false, r_solver=JuMP.optimizer_with_attributes(Ipopt.Optimizer, "tol" => TOL),
							  model_type=PM.QCLSPowerModel, reset_level=0, print_level=0)
	net_name = net["name"]

    # Get system information for relaxed AC OPF runs
	net_info = gather_net_info(net)
	num_loads = net_info["num_loads"]
	pg_max = net_info["pg_max"]
	base_load = net_info["base_load"]

	# Set pf_min & base_load in relaxed network info dict
	net_r = create_net_r(net, pf_min, pf_lagging, base_load)

	# Find the maximum servable load for each load bus in the system
	pd_max = find_max_loads(pd_max, net_r, net_path, net_name, model_type, save_max_load, print_level)

    # Initialize polytope, Ax <= b, parameter data structures
	A, b = initialize_polytope(num_loads, pd_max, pg_max, A, b, sampler)

	# Create find nearest feasible point model without the objective
	pm = PM.instantiate_model(net_r, model_type, build_opf_var_load)
	fnfp_model = find_nearest_feasible_model(pm, print_level=print_level, solver=r_solver)

	# Initialize sampler origin
	x, base_load_feasible = initialize_sample_origin(A, b, base_load, num_loads, fnfp_model, reset_level, print_level)

	return net, net_r, net_info, A, b, x
end


text:
using MixedModels
import InteractiveUtils: versioninfo
import LinearAlgebra: BLAS

# there seem to be processor-specific issues and knowing this is helpful
println(versioninfo())
@show BLAS.vendor()
if startswith(string(BLAS.vendor()), "openblas")
    println(BLAS.openblas_get_config())
end

include("utilities.jl")
include("pivot.jl")
include("UniformBlockDiagonal.jl")
include("linalg.jl")
include("matrixterm.jl")
include("FactorReTerm.jl")
include("pls.jl")
include("pirls.jl")
include("gausshermite.jl")
include("fit.jl")
include("missing.jl")
include("likelihoodratiotest.jl")
include("grouping.jl")
include("bootstrap.jl")


text:
using Test
using Scruff
using Scruff.Operators
using Scruff.SFuncs
using Scruff.Utils

@testset "Score" begin

    @testset "Computing score" begin

        @testset "Hard score" begin
            s = HardScore(2)
            @test get_score(s, 1) == 0.0
            @test get_score(s, 2) == 1.0
            @test get_log_score(s, 1) == -Inf
            @test get_log_score(s, 2) == 0.0
        end

        @testset "Soft score" begin
            s = SoftScore([:a, :b], [0.1, 0.2])
            @test isapprox(get_score(s, :a), 0.1)
            @test isapprox(get_score(s, :b), 0.2)
            @test get_score(s, :c) == 0.0
            @test isapprox(get_log_score(s, :a), log(0.1))
            @test isapprox(get_log_score(s, :b), log(0.2))
            @test get_log_score(s, :c) == -Inf
        end

        @testset "Log score" begin
            s = LogScore([:a, :b], [-1.0, -2.0])
            @test isapprox(get_score(s, :a), exp(-1.0))
            @test isapprox(get_score(s, :b), exp(-2.0))
            @test get_log_score(s, :a) == -1.0
            @test get_log_score(s, :b) == -2.0
        end

        @testset "Functional score" begin
            f(x) = 1.0 / x
            s = FunctionalScore{Float64}(f)
            @test isapprox(get_score(s, 2.0), 0.5)
            @test isapprox(get_log_score(s, 2.0), log(0.5))
        end

        @testset "Normal score" begin
            s = NormalScore(1.0, 2.0)
            r = normal_density(3.0, 1.0, 2.0)
            @test isapprox(get_score(s, 3.0), r)
            @test isapprox(get_log_score(s, 3.0), log(r))
        end

        @testset "Parzen score" begin
            s = Parzen([-1.0, 1.0], 2.0)
            r1 = normal_density(0.5, -1.0, 2.0)
            r2 = normal_density(0.5, 1.0, 2.0)
            ans = 0.5 * (r1 + r2)
            @test isapprox(get_score(s, 0.5), ans)
            @test isapprox(get_log_score(s, 0.5), log(ans))
        end

    end
end


text:
using .SampleChains: SampleChains
using .SampleChains.TupleVectors: TupleVector

# TODO: supported samples with a nested structure

function namedtuple_of_arrays(x::TupleVector{<:NamedTuple{K}}) where {K}
    return NamedTuple(k => flatten(getproperty(x, k)) for k in K)
end
function namedtuple_of_arrays(chain::SampleChains.AbstractChain)
    return namedtuple_of_arrays(SampleChains.samples(chain))
end
function namedtuple_of_arrays(multichain::SampleChains.MultiChain)
    chains = SampleChains.getchains(multichain)
    return namedtuple_of_arrays(map(namedtuple_of_arrays, chains))
end

_maybe_multichain(x) = x
_maybe_multichain(x::SampleChains.MultiChain) = x
_maybe_multichain(x::SampleChains.AbstractChain) = SampleChains.MultiChain(x)
function _maybe_multichain(x::AbstractVector{<:SampleChains.AbstractChain})
    return SampleChains.MultiChain(x)
end
function _maybe_multichain(x::Tuple{Vararg{<:SampleChains.AbstractChain}})
    return SampleChains.MultiChain(x...)
end

# info(::AbstractChain) is only required to return an AbstractVector, which is not enough
# information for us to convert it
# see https://github.com/arviz-devs/ArviZ.jl/issues/124
_samplechains_info(::SampleChains.AbstractChain) = nothing
function _samplechains_info(multichain::SampleChains.MultiChain)
    stats = map(_samplechains_info, SampleChains.getchains(multichain))
    all(isnothing, stats) && return nothing
    return namedtuple_of_arrays(stats)
end

@require SampleChainsDynamicHMC = "6d9fd711-e8b2-4778-9c70-c1dfb499d4c4" begin
    using .SampleChainsDynamicHMC: SampleChainsDynamicHMC

    function _samplechains_info(chain::SampleChainsDynamicHMC.DynamicHMCChain)
        info = SampleChains.info(chain)
        termination = info.termination
        tree_stats = (
            energy=info.π,
            tree_depth=info.depth,
            acceptance_rate=info.acceptance_rate,
            n_steps=info.steps,
            diverging=map(t -> t.left == t.right, termination),
            turning=map(t -> t.left < t.right, termination),
        )
        used_info = (:π, :depth, :acceptance_rate, :steps, :termination)
        skipped_info = setdiff(propertynames(info), used_info)
        isempty(skipped_info) ||
            @debug "Skipped SampleChainsDynamicHMC info entries: $skipped_info."
        return tree_stats
    end
end

"""
    from_samplechains(
        posterior=nothing;
        prior=nothing,
        library=SampleChains,
        kwargs...,
    ) -> InferenceData

Convert SampleChains samples to an `InferenceData`.

Either `posterior` or `prior` may be a `SampleChains.AbstractChain` or
`SampleChains.MultiChain` object.

For descriptions of remaining `kwargs`, see [`from_namedtuple`](@ref).
"""
function from_samplechains(
    posterior=nothing;
    prior=nothing,
    sample_stats=nothing,
    sample_stats_prior=nothing,
    library=SampleChains,
    kwargs...,
)
    posterior_mc = _maybe_multichain(posterior)
    prior_mc = _maybe_multichain(prior)
    if sample_stats === nothing && posterior_mc isa SampleChains.MultiChain
        sample_stats = _samplechains_info(posterior_mc)
    end
    if sample_stats === nothing && prior_mc isa SampleChains.MultiChain
        sample_stats_prior = _samplechains_info(prior_mc)
    end
    return from_namedtuple(
        posterior_mc;
        prior=prior_mc,
        sample_stats=sample_stats,
        sample_stats_prior=sample_stats_prior,
        library=library,
        kwargs...,
    )
end

"""
    convert_to_inference_data(
        obj::SampleChains.AbstractChain;
        group=:posterior,
        kwargs...,
    ) -> InferenceData
    convert_to_inference_data(
        obj::SampleChains.AbstractChain;
        group=:posterior,
        kwargs...,
    ) -> InferenceData

Convert the chains `obj` to an [`InferenceData`](@ref) with the specified `group`.

Remaining `kwargs` are forwarded to [`from_samplechains`](@ref).
"""
function convert_to_inference_data(
    chain::T; group=:posterior, kwargs...
) where {
    T<:Union{
        SampleChains.AbstractChain,
        SampleChains.MultiChain,
        AbstractVector{<:SampleChains.AbstractChain},
        Tuple{Vararg{<:SampleChains.AbstractChain}},
    },
}
    group = Symbol(group)
    group === :posterior && return from_samplechains(chain; kwargs...)
    return from_samplechains(; group => chain, kwargs...)
end


text:
using Flux, Statistics, Trebuchet, Printf
using Flux.Tracker: data, grad, gradient
using Flux.Optimise: _update_params!
using DataStructures: CircularBuffer
using Distributions: sample

#using CuArrays

#=
# Description of the problem:
    There is Trebuchet, which throws a mass to a target. The mass is to be
    released at an angle, and at certain velocity so that it lands on the target.
    The velocity of release is determined by the counterweight of the Trebuchet.
    Given conditions of environment we are required to predict the angle of
    release and counterweight.

    The problem is gamified by introducing a threshold. The player gets 99
    attempts. The first attempt has threshold value of 1. Threshold determines
    the tolerable relative error between actual distance travelled by the Projectile
    and the target distance. This relative error has to be less than the threshold
    in order to continue playing the game. The game gets tougher at each attempts
    by reduction of threshold by 0.01.

# Input:  Wind speed,   Target distance
# Output: ReleaseAngle, Weight
=#

# ----------------------------- Parameters -------------------------------------

STATE_SIZE = 2
ACTION_SIZE = 2

DIST  = (2f1, 1f2)	# Maximum target distance
SPEED =   5f0 # Maximum wind speed

ACTION_BOUND = [DIST[2]-DIST[1], SPEED]
MAX_EP = 50_000
MAX_EP_LENGTH = 1000

BATCH_SIZE = 64
MEM_SIZE = 100_000
MIN_EXP_SIZE = 50_000

γ = 99f-2     # discount rate

τ = 1f-3 # for running average while updating target networks
η_act = 1f-4   # Learning rate
η_crit = 1f-3
L2_DECAY = 1f-2

# Ornstein-Uhlenbeck Noise params
μ = 0f0
θ = 15f-2
σ = 2f-1

# --------------------------------- Memory ------------------------------------

memory = CircularBuffer{Any}(MEM_SIZE)

function getData(batch_size = BATCH_SIZE)
  # Getting data in shape
  minibatch = sample(memory, batch_size)
  x = hcat(minibatch...)

  s      =   hcat(x[1, :]...) |> gpu
  a      =   hcat(x[2, :]...) |> gpu
  r      =   hcat(x[3, :]...) |> gpu
  s′     =   hcat(x[4, :]...) |> gpu
  s_mask = .!hcat(x[5, :]...) |> gpu

  return s, a, r, s′, s_mask
end

# -------------------------------- Action Noise --------------------------------

struct OUNoise
  μ
  θ
  σ
  X
end

ou = OUNoise(μ, θ, σ, zeros(Float32, ACTION_SIZE) |> gpu)

function sample_noise(ou::OUNoise)
  dx     = ou.θ * (ou.μ .- ou.X)
  dx   .+= ou.σ * randn(Float32, length(ou.X)) |> gpu
  ou.X .+= dx
end

# Noise scale
τ_ = 25
ϵ  = exp(-1f0 / τ_)
noise_scale = 1f0

# ----------------------------- Model Architecture -----------------------------

w_init(dims...) = 6f-3rand(Float32, dims...) .- 3f-3

actor = Chain(Dense(STATE_SIZE, 400, relu),
	      	  Dense(400, 300, relu),
              Dense(300, ACTION_SIZE, tanh, initW=w_init)) |> gpu
actor_target = deepcopy(actor)

struct crit
  state_crit
  act_crit
  sa_crit
end

Flux.@treelike crit

function (c::crit)(state, action)
  s = c.state_crit(state)
  a = c.act_crit(action)
  c.sa_crit(relu.(s .+ a))
end

Base.deepcopy(c::crit) = crit(deepcopy(c.state_crit),
                              deepcopy(c.act_crit),
			      			  deepcopy(c.sa_crit))

critic = crit(Chain(Dense(STATE_SIZE, 400, relu), Dense(400, 300)) |> gpu,
              	  	Dense(ACTION_SIZE, 300) |> gpu,
	      			Dense(300, 1, initW=w_init) |> gpu)
critic_target = deepcopy(critic)

# ------------------------------- Param Update Functions---------------------------------

function update_target!(target, model; τ = 1f0)
  for (p_t, p_m) in zip(params(target), params(model))
    p_t.data .= (1f0 - τ) * p_t.data .+ τ * p_m.data
  end
end

nullify_grad!(p) = p
nullify_grad!(p::TrackedArray) = (p.grad .= 0f0)

zero_grad!(model) = (model = mapleaves(nullify_grad!, model))

# ---------------------------------- Training ----------------------------------

function L2_loss(model)
  l2_loss = 0
  for p in params(model)
    l2_loss += sum(p .^ 2)
  end

  return L2_DECAY * l2_loss
end

opt_crit = ADAM(η_crit)
opt_act  = ADAM(η_act)

function replay()
  s, a, r, s′, s_mask = getData()

  # Update Critic
  a′ = data(actor_target(s′))
  v′ = data(critic_target(s′, a′))
  y = r .+ γ * v′ .* s_mask	# set v′ to 0 where s_ is terminal state
  @show size(s)
  v = critic(s, a)
  loss_crit = Flux.mse(y, v) + L2_loss(critic)

  zero_grad!(critic)
  Flux.back!(loss_crit)
  _update_params!(opt_crit, params(critic))

  # Update Actor
  actions = actor(s)
  crit_out = critic(s, actions)

  zero_grad!(actor)
  Flux.back!(-sum(crit_out))
  _update_params!(opt_act, params(actor))

  # Update Target models
  update_target!(actor_target, actor; τ = τ)
  update_target!(critic_target, critic; τ = τ)
end

# --------------------------- Helper Functions --------------------------------

lerp(x, lo, hi) = x*(hi-lo)+lo

target() = [randn() * SPEED, lerp(rand(), DIST...)]

# Stores tuple of state, action, reward, next_state, and done
remember(state, action, reward, next_state, done) =
  push!(memory, [data.((state, action, reward[1], next_state))..., done])

function action(state, train=true)
  state = reshape(state, size(state)..., 1)
  act_pred = data(actor(state |> gpu)) |> cpu
  if train
    act_pred .+= noise_scale * sample_noise(ou)
  end
  angle = clamp.(act_pred[1:1, :], -π, π) # returns action

  return vcat(angle, act_pred[2:2, :])
end

function shoot(wind, angle, weight)
  Trebuchet.shoot((wind, Trebuchet.deg2rad(angle), weight))[2]
end

function reward(wind, target_dist, angle, weight)
  actual_dist = shoot(wind, angle, weight)
  -(target_dist - actual_dist) ^ 2
end

function step!(state, action)
  r = reward(state..., action...)
  return r
end

function episode(train=true)
  wind, target_dist = target()
  s = [wind_speed, target_dist]
  a = action(s, train)
  r = step!(s, a)

  if train
    remember(s, a, zeros(Float32, 2), r, true)
    replay()
  end

  return √-r
end
# -------------------------------- Testing -------------------------------------

# Returns average score over 100 episodes
function test()
  score_mean = 0f0
  for _=1:100
    total_reward = episode(false)
    score_mean += total_reward / 100
  end
  return score_mean
end

# ------------------------------ Training --------------------------------------

# Populate memory with random actions

for e=1:MIN_EXP_SIZE
  s = target()
  a = 2rand(Float32, ACTION_SIZE) .* ACTION_BOUND .- ACTION_BOUND
  r = step!(s, a)
  remember(s, a, r, zeros(Float32, 2), true)
end

for e=1:MAX_EP
  global noise_scale
  total_reward = episode(true)
  total_reward = @sprintf "%9.3f" total_reward
  print("Episode: $e | Score: $total_reward | ")
  score_mean = test()
  score_mean = @sprintf "%9.3f" score_mean
  println("Mean score over 100 test episodes: $score_mean")
  noise_scale *= ϵ
end

max_stars_repo_name: c-p-murphy/model-zoo

max_issues_repo_name: c-p-murphy/model-zoo

max_forks_repo_name: c-p-murphy/model-zoo


text:
include("common.jl")

mmesh = Mesh(joinpath(PDATA,"crack.stl"))
c1=(0.50, 0.0)
c2=(1.5, 0.5)
fault_idx, fault_coords, fault_normal_vec = get_fault_info(mmesh, c1, c2)

close("all")
visualize_mesh(mmesh)
plot(mmesh.nodes[fault_idx,1], mmesh.nodes[fault_idx,2], ".")
savefig("mesh1.png")

mmesh, conn = fix_mesh(mmesh, fault_idx, fault_redblack)
cn1 = conn[:,1]
cn2 = conn[:,2]
close("all")
visualize_mesh(mmesh)
v = 0.01*(rand(length(cn1)).-0.5)
u = 0.01*(rand(length(cn1)).-0.5)
scatter(mmesh.nodes[cn1,1] + u, mmesh.nodes[cn1,2] + v, s = 2, color="red")
scatter(mmesh.nodes[cn2,1], mmesh.nodes[cn2,2], s = 2, color="blue")

savefig("mesh.png")

bdnode = bcnode((x, y)->(x<1e-5 || x>2-1e-5), mmesh)

Conn = compute_slip_boundary_condition( conn, mmesh)
NT = 100
Δt = 0.004
λ = 1.0
μ = 1.0
invη = constant(ones(get_ngauss(mmesh)))
sn = size(conn,1)
slipvec = constant(0.1*[2ones(NT, sn) ones(NT, sn)]) # acceleration
gravity = 0.0

U, Sigma, Varepsilon = subduction_solver(invη, gravity, Conn, slipvec, bdnode, mmesh)

sess = Session(); init(sess)
Uval, Stress, Strain = run(sess, [U, Sigma, Varepsilon])
p = visualize_displacement(Uval[1:5:end,:], mmesh)
saveanim(p, "subduction_rectangle.gif")

loss = sum(U^2)
close("all")
u = Uval[end,:]
visualize_scalar_on_fem_points(u[mmesh.ndof+1:2mmesh.ndof],mmesh)
savefig("y_displacement.png")

close("all")
p = visualize_scalar_on_gauss_points(Stress[end,:,1], mmesh)
savefig("fault_stress.png")


text:
export SimplexExact, transferoperator

"""
    SimplexExact() <: TriangulationBasedTransferOperator

A transfer operator estimator using a triangulation partition and exact
simplex intersections[^Diego2019].

To use this estimator, the Simplices.jl
package must be brought into scope by doing `using Simplices` after
`using PerronFrobenius`. *Note: due to computing the exact simplex intersection,
this estimator is slow compared to `SimplexApprox`.*

[^Diego2019]: Diego, David, Kristian Agasøster Haaga, and Bjarte Hannisdal. "Transfer entropy computation using the Perron-Frobenius operator." Physical Review E 99.4 (2019): 042212.
"""
struct SimplexExact <: TriangulationBasedTransferOperator
    bc::String

    function SimplexExact(bc::String = "circular")
        isboundarycondition(bc, "triangulation")  || error("Boundary condition '$bc' not valid.")
        new(bc)
    end
end
Base.show(io::IO, se::SimplexExact) = print(io, "SimplexExact{$(se.bc)}")

""" Generate a transopergenerator for an exact simplex estimator."""
function transopergenerator(pts, method::SimplexExact)
    # modified points, where the image of each point is guaranteed to lie within the convex hull of the previous points
    invariant_pts = invariantize(pts)

    # triangulation of the invariant points. The last point is excluded, so that
    # the last vertex also can be mapped forward one step in time.
    triang = DelaunayTriangulation(invariant_pts[1:end-1])

    init = (invariant_pts = invariant_pts, triang = triang)

    TransferOperatorGenerator(method, pts, init)
end

function (tog::TransferOperatorGenerator{<:SimplexExact})(; tol = 1e-8)
    invariant_pts, triang = getfield.(Ref(tog.init), (:invariant_pts, :triang))

    D = length(invariant_pts[1])
    N = length(triang)
    ϵ = tol / N

    # Pre-allocate simplex and its image
    image_simplex = SimplexTypes.MutableSimplex(zeros(Float64, D, D+1))
    simplex = SimplexTypes.MutableSimplex(zeros(Float64, D, D+1))
    M = zeros(Float64, N, N)

    for j in 1:N
        for k = 1:(D + 1)
            # Get the vertices of the image simplex
            image_simplex[k] .= invariant_pts[triang[j][k] + 1]
        end

        imvol = abs(SimplexTypes.orientation(image_simplex))

        for i = 1:N
            for k = 1:(D + 1)
                # Get the vertices of the source simplex
                simplex[k] .= invariant_pts[triang[i][k]]
            end

            # Only compute the entry of the transfer matrix
            # if simplices are of sufficient size.
            vol = abs(SimplexTypes.orientation(simplex))

            if vol * imvol > 0 && (vol/imvol) > ϵ
                M[j, i] = SimplexTypes.intersect(simplex, image_simplex) / imvol
            end
        end
    end
    params = (; tol = tol)
    TransferOperatorApproximation(tog, M, params)
end

function transferoperator(pts, method::SimplexExact; tol::Real = 1e-8)
    tog = transopergenerator(pts, method)
    tog(tol = tol)
end


text:
#----------------- Get the A matrix
function getDivGrad(n1,n2,n3)

        # the Divergence
        D1 = kron(sparse(I,n3,n3),kron(sparse(I,n2,n2),ddx(n1)))
        D2 = kron(sparse(I,n3,n3),kron(ddx(n2),sparse(I,n1,n1)))
        D3 = kron(ddx(n3),sparse(I,n1*n2,n1*n2))
        # DIV from faces to cell-centers
        Div = [D1 D2 D3]

        return Div*Div';
end
#----------------- 1D finite difference on staggered grid
function ddx(n)
# generate 1D derivatives
	II, JJ, VV = SparseArrays.spdiagm_internal(0 => -ones(n), 1 => ones(n));
	return sparse(II, JJ, VV, n, n+1)

end


text:
using Pkg
Pkg.activate(pwd())

# # Plots.jl

using Plots

x = range(0, 2π; length = 100)
y = sin.(x)

plot(x, y)

#+

y = hcat(sin.(x), cos.(x))

plot(x, y);
plot!(x, sin.(x .+ π/4))

#+

plt = plot(x, hcat(sin.(x), cos.(x)));
plot!(plt, x, sin.(x .+ π/4))

# ## Plot attributes

x = range(0, 2π; length = 100)
y = hcat(sin.(x), cos.(x))

plot(x, y;
    label = ["sine" "cosine"],
    xguide = "x",
    yguide = "y",
    legend = :bottomleft,
    title = "Trigonometric functions",
    xticks = (0:0.5π:2π, ["0", "0.5π", "π", "1.5π", "2π"]),
    color = [:red :blue],
    linestyle = [:dash :dot],
    linewidth = [2 4],
)

#+

n = 200
x = range(0, 2π; length = n)
linewidth = range(1, 50; length = n)
color = palette(:viridis, n)
xlims = (0, 7)
ylims = (-1.2, 1.2)
label = ""

plot(x, sin.(x); linewidth, color, label, xlims, ylims)

#+

plot(x, [sin.(x) cos.(x)]; linewidth, color = [color :red], label, xlims, ylims)

#+

plotattr(:Series)

#+

plotattr("title")

#+

plotattr("xguide")

#+

plotattr("guide")

# ### Exercise:
# Consider the following set of equations
#
# $$
# \begin{aligned}
# x(t) & = \cos(3t), \\
# y(t) & = \sin(2t),\\
# \end{aligned}
# $$
#
# where $t \in [0, 2\pi]$. Create a plot of the curve described by the equations above.
# Use plot attributes to set the following properties:
#
# 1. The line width should start at `1`, increase to `50` and then decrease back to `1`.
# 2. The line color should change with the changing line width.
#
# Use `:viridis` color scheme or any other
# [color scheme](http://docs.juliaplots.org/latest/generated/colorschemes/) supported by
# the Plots package. Use additional plot attributes to get a nice looking graph.
#
# **Hints:**
# - use the `pallete` function combined with the `collect` function to generate a vector of
# colors from the `:viridis` color scheme.
# - remove all decorators by using: `axis = nothing`, `border = :none`.
#
# ---
# ### Solution:

n = 1000
t = range(0, 2π; length = n)

fx(t) = cos(3t)
fy(t) = sin(2t)

#+

linewidth = vcat(
    range(1, 50; length = n ÷ 2),
    range(50, 1; length = n - n ÷ 2)
)

#+

c = palette(:viridis, 2);
typeof(c)

#+

c = collect(palette(:viridis, 2))

#+

color = vcat(
    collect(palette(:viridis, n ÷ 2)),
    collect(palette(:viridis, n - n ÷ 2; rev = true))
)

#+

plot(fx.(t), fy.(t);
    linewidth,
    color,
    lims = (-1.2, 1.2),
    legend = false,
    axis = nothing,
    border = :none,
)

# ---
#
# ## Function plotting

t = range(0, 2π; length = 100)

plot(t, [sin, cos]; label = ["sine" "cosine"])

#+

plot(sin, x -> sin(2x), t; linewidth = 2, label = "")

#+

plot(sin, x -> sin(2x), 0, 2π, 100; linewidth = 2, label = "")

# ### Exercise:
# Create a plot given by the following set of equations:
#
# $$
# \begin{aligned}
# x(t) & = (a + b)\cos(t) - b \cdot \cos \left( \left(\frac{a}{b} + 1 \right)t \right), \\
# y(t) & = (a + b)\sin(t) - b \cdot \sin \left( \left(\frac{a}{b} + 1 \right)t \right), \\
# \end{aligned}
# $$
# where $a = 4.23$, $b = 2.35$ and $t \in [-15, 20]$. Use additional plot attributes
# to get a nicely looking graph.
#
# ---
# ### Solution:

fx(t; a = 4.23, b = 2.35) = (a + b)*cos(t) - b*cos((a/b + 1)*t)
fy(t; a = 4.23, b = 2.35) = (a + b)*sin(t) - b*sin((a/b + 1)*t)

plot(fx, fy, -15, 20, 500;
    linewidth = 2,
    legend = false,
    axis = nothing,
    border = :none,
)

# ---
#
# Changing the plotting series

x = range(0, 2π; length = 100)
y = sin.(x)

plot(x, y; seriestype = :scatter)

#+

scatter(x, y)

# ### Exercise:
# Consider the following function:
#
# $$
# f(x, y) = \frac{x^2 \cdot y^2}{x^4 + y^4}.
# $$
#
# Draw this function for $x, y \in [-5, 5]$. Use the following three plot series
# `contourf`, `heatmap`, and `surface` with the following settings:
# - `:viridis` color scheme,
# - camera angle `(25, 65)`,
# - no legend, color bar, or decorators (`axis`, `frame` and `ticks`).
#
# ---
# ### Solution:

x = range(-5, 5; length = 200)

fz(x, y) = x^2*y^2/(x^4 + y^4)

#+

kwargs = (
    color = :viridis,
    legend = false,
    cbar = false,
    axis = nothing,
    border = :none,
)

#+

contourf(x, x, fz; kwargs...) # or plot(x, x, fz; seriestype = :contourf, kwargs...)
heatmap(x, x, fz; kwargs...)
surface(x, x, fz; camera = (25, 65), kwargs...)

# ---
#
# Subplots

x = range(0, 2π; length = 100)

plot(x, [sin, cos, tan, sinc];
    layout = 4,
    linewidth = 2,
    legend = false,
    title = ["1" "2" "3" "4"],
)

#+

plot(x, [sin, cos, tan, sinc];
    layout = grid(4, 1; heights = [0.1 ,0.4, 0.4, 0.1]),
    linewidth = 2,
    legend = false,
)

#+

l = @layout [a ; b{0.3w} c]
plot(x, [sin, cos, tan]; layout = l, linewidth = 2, legend = false)

#+

linewidth = range(1, 20; length = 100)

p1 = plot(x, sin; legend = false, line_z = 1:100, color = :viridis, linewidth);
p2 = plot(x, cos; legend = false, line_z = 1:100, color = :Blues_9, linewidth);
p3 = plot(x, tan; legend = false, line_z = 1:100, color = :hsv, linewidth);

l = @layout [a ; b{0.3w} c]
plot(p1, p2, p3; layout = l)

# ## Animations

n = 300

plt = plot(Float64[], [sin, cos];
    legend = false,
    xlims = (0, 6π),
    ylims = (-1.1, 1.1),
    linewidth = range(1, 20; length = n),
    color = palette(:viridis, n),
    axis = nothing,
    border = :none
)

anim = Animation()

for x in range(0, 6π; length = n)
    push!(plt, x, [sin(x), cos(x)])
    frame(anim)
end

gif(anim, joinpath(pwd(), "lecture_05", "animsincos.gif"), fps = 15)

#+

x = range(-5, 5; length = 400)
fz(x, y) = x^2*y^2/(x^4 + y^4)

plt = surface(x, x, fz;
    camera = (30, 65),
    color = :viridis,
    legend = false,
    axis = nothing,
    border = :none,
    cbar = false,
)

anim = @animate for i in vcat(30:60, 60:-1:30)
    plot!(plt, camera = (i, 65))
end

gif(anim, joinpath(pwd(), "lecture_05", "animsurf.gif"), fps = 15)


text:
"""
    publish_post!(state, post_list, agent_idx, tick_nr)
Publish a post to the network.
# Arguments
- `state`: a tuple of the current graph and agent_list
- `post_list`: List of all published posts in network
- `agent_idx`: agent index
- `tick_nr`: Number of current simulation tick
See also: [`Config`](@ref), [`Agent`](@ref)
"""
function publish_post!(
    state::Tuple{AbstractGraph, AbstractArray},
    post_list::AbstractArray,
    agent_idx::Integer,
    tick_nr::Integer=0
)

    graph, agent_list = state
    this_agent = agent_list[agent_idx]

    post_opinion = this_agent.opinion + 0.1 * (2 * rand() - 1)

    if post_opinion > 1
        post_opinion = 1.0
    elseif post_opinion < -1
        post_opinion = -1.0
    end

    post = Post(
        post_opinion,
        length(outneighbors(graph, agent_idx)),
        agent_idx,
        tick_nr
    )

    push!(post_list, post)

    for neighbor in outneighbors(graph, agent_idx)
        push!(agent_list[neighbor].feed, post)
    end

    return state, post_list
end

# suppress output of include()
;


text:
using QuDynamics.QuTiP
begin

# Defining the system - Hamiltonian, initial state, time steps
# Hamiltonian of the system
H = sigmax

# Initial state of the system
init_state = statevec(1, FiniteBasis(2))
init_state_dm = init_state*init_state'

# Time steps
tlist = 0.:0.1:2*pi

# Collapse operators
c_ops = [lowerop(2)]

# Operators whose expectation values are to computed
# in the evolved states
e_ops = [sqrt(2)*lowerop(2) + sqrt(3)*raiseop(2), raiseop(2), lowerop(2)]

# Test set involving Schrodinger like equations
se = QuSchrodingerEq(sigmax)

# Evaluating states using the general iterator version
evolved_states_se = Array(QuBase.AbstractQuArray, length(tlist)-1)
i = 1
for (t, psi) in QuStateEvolution(sigmax, init_state, tlist, QuExpmV())
    evolved_states_se[i] = psi
    i = i + 1
end

ind = rand(1:(length(tlist)-1))

# Evaluating using sesolve by passing Hamiltonian and initial state vector
schro_direct = sesolve(sigmax, init_state, tlist)

# Evaluating using sesolve by passing Schrodinger type equation
schro_eq = sesolve(se, init_state, tlist)

# Tests comparing the states from all the above constructs
@assert schro_direct[2][ind] == evolved_states_se[ind]
@assert schro_eq[2][ind] == evolved_states_se[ind]
@assert schro_direct[2][ind] == schro_eq[2][ind]

# Evaluating the expectation values
schro_direct_expect = sesolve(sigmax, init_state, tlist, e_ops)[2]
schro_eq_expect = sesolve(se, init_state, tlist, e_ops)[2]

ind_1 = rand(1:(length(e_ops)))
ind_2 = rand(1:length(tlist)-1)

@assert schro_direct_expect[ind_2, ind_1].im == schro_eq_expect[ind_2, ind_1].im

# Test set involving LiouvillevonNeumann like equations
lvne = QuLiouvillevonNeumannEq(sigmax)

# Evaluating states using the general iterator version
evolved_states_lvne = Array(QuBase.AbstractQuArray, length(tlist)-1)
i = 1
for (t, psi) in QuStateEvolution(sigmax, init_state_dm, tlist, QuExpmV())
    evolved_states_lvne[i] = psi
    i = i + 1
end

ind=rand(1:(length(tlist)-1))

# Evaluating using sesolve by passing Hamiltonian and initial density matrix
lvne_direct = sesolve(sigmax, init_state_dm, tlist)

# Evaluating using sesolve by passing LiouvillevonNeumann type equation
lvne_eq = sesolve(lvne, init_state_dm, tlist)

# Tests comparing the states from all the above constructs
@assert lvne_direct[2][ind] == evolved_states_lvne[ind]
@assert lvne_eq[2][ind] == evolved_states_lvne[ind]
@assert lvne_direct[2][ind] == lvne_eq[2][ind]

# Evaluating the expectation values
expectation_lvne = QuDynamics.QuTiP.eops_helper(e_ops, evolved_states_lvne)

# Evaluating using sesolve
lvne_direct_expect = sesolve(sigmax, init_state_dm, tlist, e_ops)[2]
lvne_eq_expect = sesolve(lvne, init_state_dm, tlist, e_ops)[2]

ind_1 = rand(1:(length(e_ops)))
ind_2 = rand(1:length(tlist)-1)

@assert lvne_direct_expect[ind_2, ind_1].im == expectation_lvne[ind_2, ind_1].im
@assert lvne_eq_expect[ind_2, ind_1].im  == expectation_lvne[ind_2, ind_1].im
@assert lvne_direct_expect[ind_2, ind_1].im == lvne_eq_expect[ind_2, ind_1].im

# Test set involving Lindblad Master equations
lme = QuLindbladMasterEq(sigmax, c_ops)

# Evaluating states using the general iterator version
evolved_states_lme = Array(QuBase.AbstractQuArray, length(tlist)-1)
i = 1
for (t, psi) in QuStateEvolution(lme, init_state_dm, tlist, QuExpmV())
    evolved_states_lme[i] = psi
    i = i + 1
end

ind=rand(1:(length(tlist)-1))

# Evaluating using mesolve by passing Hamiltonian and initial density matrix
lme_direct = mesolve(sigmax, init_state_dm, tlist, c_ops)

# Evaluating using mesolve by passing Lindblad Master equation
lme_eq = mesolve(lme, init_state_dm, tlist)

# Tests comparing the states from all the above constructs
@assert lme_direct[2][ind] == evolved_states_lme[ind]
@assert lme_eq[2][ind] == evolved_states_lme[ind]
@assert lme_direct[2][ind] == lme_eq[2][ind]

# Evaluating the expectation values
expectation_lme = QuDynamics.QuTiP.eops_helper(e_ops, evolved_states_lme)

# Evaluating using mesolve
lme_direct_expect = mesolve(sigmax, init_state_dm, tlist, c_ops, e_ops)[2]
lme_eq_expect = mesolve(lme, init_state_dm, tlist, e_ops)[2]

ind_1 = rand(1:(length(e_ops)))
ind_2 = rand(1:length(tlist)-1)

@assert lme_direct_expect[ind_2, ind_1].im == expectation_lme[ind_2, ind_1].im
@assert lme_eq_expect[ind_2, ind_1].im == expectation_lme[ind_2, ind_1].im
@assert lme_direct_expect[ind_2, ind_1].im == lme_eq_expect[ind_2, ind_1].im

# Testing the expectation using various mcsolve constructs
mcs = mcsolve(H, init_state, tlist, c_ops)
mcs_lme = mcsolve(lme, init_state, tlist)
mcs_lme_expect = mcsolve(lme, init_state_dm, tlist, e_ops)
mcs_expect = mcsolve(H, init_state, tlist, c_ops,  e_ops)

@test_approx_eq_eps mcs[2][ind_2] mcs_lme[2][ind_2] 1e-1
@test_approx_eq_eps trace(mcs[2][ind_2]*e_ops[ind_1]) mcs_expect[2][ind_2, ind_1] 1e-1
@test_approx_eq_eps trace(mcs[2][ind_2]*e_ops[ind_1]) mcs_lme_expect[2][ind_2, ind_1] 1e-1
@test_approx_eq_eps mcs_expect[2][ind_2, ind_1] mcs_lme_expect[2][ind_2, ind_1] 1e-1
end


text:
# Jonas August 2021

#=

This problem was recently asked by Google.

Given a list of numbers and a number k, return whether any two numbers from the list add up to k.

For example, given [10, 15, 3, 7] and k of 17, return true since 10 + 7 is 17.

Bonus: Can you do this in one pass?

=#

# Two pass method
function haspair2(list, sum)
   for i in 1:length(list)
      for j in 1:length(list)
        if j == 1
	    continue
	end
	if list[i] + list[j] == sum
	    return true
	end
      end
   end
   false
end

#=
TODO: single-pass
I can see quadratic-time method with only a single pass of list, by introducing a pass on a different array, so not really an improvement IMO.
=#

using Test

@test haspair2([10, 15, 3, 7], 17)
@test !haspair2([10, 15, 3, 7], 26)

max_stars_repo_path: Pair sum - Problem 1.jl

max_issues_repo_path: Pair sum - Problem 1.jl

max_forks_repo_path: Pair sum - Problem 1.jl


text:
# TODO: Outsource EPGs and rename MRIQuant -> MRIRelax
module ExtendedPhaseGraphs

	export required_kmax

	using LinearAlgebra
	using MRIQuant: γ

	@inline function required_kmax(num_TR::Integer)::Integer
		@assert num_TR >= 0
		(num_TR ÷ 2) - 1
	end

	# Code for computing which k-states need to be considered
	# +1 in the following are due to indices starting at 1
	# No argument checks
	@inline required_states(mode::Val{:minimal}, t::Integer, N::Integer, kmax::Integer)::Integer = min(t, N-t+1, kmax+1)
	@inline required_states(mode::Val{:full}, t::Integer, N::Integer, kmax::Integer)::Integer = kmax+1
	@inline required_states(mode::Val{:full_in}, t::Integer, N::Integer, kmax::Integer)::Integer = min(kmax+1, N-t+1)
	@inline required_states(mode::Val{:full_out}, t::Integer, N::Integer, kmax::Integer)::Integer = min(t, kmax+1)

	"""
		Basically b-values?

		gradient G in mT
		duration τ in ms
		k in [index]
		K0 in rad / m
		Apply function multiple times and set proper Δk

		Definition
		k: index in phase graph, multiples of the total shift
		produced by all gradients together, unit is index.
		K: The shift in physical units (radians/m), can be non-integer
		multiples of the shift produces by all gradients together.
		sum(τ) = TR

		TODO: Split into individual gradients, because rf pulses might come inbetween!

		returns b_longitudinal, b_transverse in units of (rad/m)^2 * s
	"""
	function diffusion_b_values(
		G::AbstractVector{<: Real},
		τ::AbstractVector{<: Real},
		kmax::Integer
	)::Tuple{Vector{Float64}, Vector{Float64}}

		@assert all(τ .> 0)
		@assert kmax .>= 0

		# Shift in K generated by individual gradients
		dK = γ .* G .* τ
		# Shift by all gradients together
		ΔK = sum(dK)

		# For longitudinal states, K remains constant even when gradients are non-zero.
		b_longitudinal = (0:kmax).^2 .* (ΔK^2 * sum(τ))

		# For transverse states, K changes while the gradients are non-zero.
		# Thus, individual parts of the gradient shape cannot be considered independently.
		# The previous shift K0 must be taken into account.

		# Precalculate some terms
		#=
			b_t	= b_l
				  + k * [ΔK * sum_{gradients} τ (2K0 + dK)]				# =: k * χ
				  + [sum_{gradients} τ (K0^2 + dK^2 / 3 + K0 * dK)]		# =: ξ
		=#
		local ξ, χ
		let
			K0 = Vector{Float64}(undef, length(G))
			tmp = Vector{Float64}(undef, length(G))
			# Compute K0, inital K produced by each gradient lobe
			K0[1] = 0
			K0[2:end] = cumsum(@view dK[1:end-1])
			# Compute elementwise ξ = sum_{gradients} τ * (K0^2 + 3 dK^2 + K0 * dK)
			@. tmp = K0^2 + dK^2 / 3 + K0 * dK
			tmp .*= τ # TODO: Needs to be split bc otherwise julia allocates memory, why?
			ξ = sum(tmp)
			# Reuse tmp to compute χ = ΔK * sum_{gradients} τ * (2K0 + dK)
			@. tmp = τ * (2K0 + dK)
			χ = ΔK * sum(tmp)
		end

		b_transverse = @. b_longitudinal + χ * (0:kmax) + ξ

		return b_longitudinal, b_transverse
	end

	@inline function diffusion_factor(b::Real, D::Real)
		# if b in units of (rad/m)^2 * s, then D in units of m^2/s (magnitude 1e-9)
		exp(-b * D);
	end

	function compute_combined_relaxation(
		TR::Real,
		kmax::Integer,
		G::AbstractVector{<: Real},
		τ::AbstractVector{<: Real},
		D::Real,
		R1::Real,
		R2::Real
	)
		# Compute combined effect of T1, T2 and diffusion
		E1 = exp(-TR * R1)
		E2 = exp(-TR * R2)
		ED_longitudinal, ED_transverse = diffusion_b_values(G, τ, kmax) # Misuse of names
		@. ED_longitudinal = E1 * diffusion_factor(ED_longitudinal, D)
		@. ED_transverse = E2 * diffusion_factor(ED_transverse, D)
		return 1.0 - E1, ED_longitudinal, ED_transverse
	end

	function compute_combined_relaxation(
		TR::AbstractVector{<: Real},
		kmax::Integer,
		G::AbstractVector{<: Real},
		τ::AbstractVector{<: Real},
		D::Real,
		R1::Real,
		R2::Real
	)
		# Compute exponentials for T1, T2 relaxation factors
		E1 = exp.(-TR * R1)
		E2 = exp.(-TR * R2)

		# Compute diffusion damping factors (depends on k)
		D_longitudinal, D_transverse = diffusion_b_values(G, τ, kmax) # Misuse of names
		D_longitudinal = diffusion_factor.(D_longitudinal, D)
		D_transverse = diffusion_factor.(D_transverse, D)

		return E1, E2, D_longitudinal, D_transverse
	end

	@inline function apply_combined_relaxation!(
		state::AbstractMatrix{<: Complex},
		E1::Real,
		E2::Real,
		D_longitudinal::AbstractVector{<: Real},
		D_transverse::AbstractVector{<: Real}
	)
		@inbounds begin
			state[:, 1] .*= E2 * D_transverse
			state[:, 2] .*= E2 * D_transverse
			state[:, 3] .*= E1 * D_longitudinal
			state[1, 3] .+= 1.0 - E1
		end
		return
	end

	@inline function apply_combined_relaxation!(
		state::AbstractMatrix{<: Complex},
		upper::Integer,
		one_minus_E1::Real,
		ED_longitudinal::AbstractVector{<: Real},
		ED_transverse::AbstractVector{<: Real}
	)
		@inbounds begin
			for (s, ED) in enumerate((ED_transverse, ED_transverse, ED_longitudinal))
				for i = 1:upper
					state[i, s] *= ED[i]
				end
			end
			state[1, 3] += one_minus_E1
		end
		return
	end

	function rf_pulse_matrix(
		α::Real,
		ϕ::Real;
		out=Matrix{ComplexF64}(undef, 3, 3)
	)::Matrix{ComplexF64}
		# This is the transposed version of the one from Weigel's paper
		phase_factor = exp(im * ϕ)
		phase_factor_squared = phase_factor^2
		im_phase_factor = im * phase_factor
		@inbounds begin
			out[1, 1] = cos(0.5α)^2
			out[2, 1] = sin(0.5α)^2 * phase_factor_squared;
			out[3, 1] = -sin(α) * im_phase_factor
			out[1, 2] = conj(out[2, 1])
			out[2, 2] = out[1, 1]
			out[3, 2] = conj(out[3, 1])
			out[1, 3] = -0.5 * out[3, 2]
			out[2, 3] = conj(out[1, 3])
			out[3, 3] = cos(α)
		end
		return out
	end

	"""
		TODO: Debug version where all states are computed (the ones in the second half that cannot contribute)
		Version where not everything is computed

		TODO: More complicated sequences with different gradient moments

		α, ϕ in radians
		TR in the unit determined by R1,2
		kmax in indices
		G in mT/m
		τ in ms
		D in m^2 / time unit matching TR (magnitude 10^-12)
	"""
	@generated function simulate(
		mode::Union{Val{:minimal}, Val{:full}, Val{:full_in}, Val{:full_out}};
		initial_state::Union{Nothing, AbstractMatrix{ComplexF64}} = nothing,
		α::AbstractVector{<: Real},
		ϕ::AbstractVector{<: Real},
		TR::Union{Real, AbstractVector{<: Real}},
		kmax::Integer,
		G::AbstractVector{<: Real},
		τ::AbstractVector{<: Real},
		D::Real,
		R1::Real,
		R2::Real,
		record::Union{Val{:signal}, Val{:all}} = Val(:signal)
	)
		# Check if initial state and mode are compatible
		use_initial_state = mode <: Union{Val{:full}, Val{:full_in}}
		if initial_state <: Nothing
			@assert !use_initial_state "No initial state provided"
			setup_state = quote
				# Use two allocated memory blocks
				allocated_states = (
					zeros(ComplexF64, kmax+1, 3), # +1 for k = 0
					zeros(ComplexF64, kmax+1, 3)
				)
				allocated_states[2][1, 3] = 1 # Set M0, start with second block
			end
		else
			@assert use_initial_state "Initial state ignored"
			setup_state = quote
				@assert size(initial_state) == (kmax+1, 3) # +1 for k = 0
				# Use two allocated memory blocks
				allocated_states = (
					copy(initial_state),
					copy(initial_state)
				)
			end
		end

		# Distinguish case of scalar respectively vector-valued TR
		if TR <: AbstractVector{<: Real}
			TR_check = :(N == length(TR) && all(TR .> 0))
			combined_relaxation = :((E1, E2, D_longitudinal, D_transverse))
			apply_combined_relaxation = quote
				@inbounds apply_combined_relaxation!(state, E1[t], E2[t], ED_longitudinal, ED_transverse)
			end
		else
			TR_check = :(TR > 0)
			combined_relaxation = :((one_minus_E1, ED_longitudinal, ED_transverse))
			apply_combined_relaxation = quote
				apply_combined_relaxation!(state, upper, one_minus_E1, ED_longitudinal, ED_transverse)
			end
		end

		# Code to set up variables to record the state
		if record <: Val{:signal}
			setup_recorded_states = :(Vector{ComplexF64}(undef, N))
			record_state = :(@inbounds recorded_states[t] = state[1, 2]) # Record F^-(k = 0)
		elseif record <: Val{:all}
			setup_recorded_states = :(Array{ComplexF64}(undef, kmax+1, 3, N))
			record_state = :(@inbounds recorded_states[:, :, t] = state)
		end

		# Code to apply the rf pulse matrix
		if mode <: Val{:full}
			apply_rf_pulse_matrix = quote
				mul!(state, allocated_states[3-j], T)
			end
		else
			apply_rf_pulse_matrix = quote
				@inbounds @views mul!(state[1:upper, :], allocated_states[3-j][1:upper, :], T)
			end
		end

		# Code to apply the return the final state if desired
		if mode <: Val{:full_out} || mode <: Val{:full}
			return_value = :((recorded_states, state))
		else
			return_value = :(recorded_states)
		end

		# Interpolate and return full code
		return quote

			# Check arguments
			N = length(α)
			@assert N == length(ϕ)
			@assert $TR_check
			@assert kmax >= 0

			# Precompute relaxation (T1, T2 and diffusion)
			$combined_relaxation = compute_combined_relaxation(TR, kmax, G, τ, D, R1, R2)

			# Array holding the current phase graph state
			# but indices start at 1, while k starts at zero
			# First index for k, second index for components, basically F^+(k), F^-(-k), Z(k) stacked
			$setup_state

			# Allocate memory
			recorded_states = $setup_recorded_states
			T = Matrix{ComplexF64}(undef, 3, 3)

			# Simulate
			local state
			for t = 1:N
				# How many states are involved?
				upper = required_states(mode, t, N, kmax)

				# Apply pulse
				rf_pulse_matrix(α[t], ϕ[t]; out=T)
				j = mod1(t, 2)
				state = allocated_states[j]
				$apply_rf_pulse_matrix

				# Store signal (F^-(k = 0)) or complete state
				$record_state
				# Note: No need to apply any relaxation before getting the signal, because it is only a scaling.

				# Apply relaxation (T1, T2 and diffusion)
				$apply_combined_relaxation

				# Shift (by gradients)
				@inbounds begin
					# F^+(k) goes to F^+(k+1)
					for i = min(upper+1, kmax+1):-1:2
						state[i, 1] = state[i-1, 1]
					end
					# F^+(k = 0) becomes the conjugate of F^-(k = -1)
					state[1, 1] = conj(state[2, 2])
					# F^-(k) also goes to F^-(k+1), but k is ordered the other way around in the array
					for i = 1:(upper-1)
						state[i, 2] = state[i+1, 2]
					end
					# F^-(k = -upper+1) becomes unpopulated
					state[upper, 2] = 0
					# Longitudinal states are not getting shifted
				end
			end

			return $return_value
		end
	end

end


text:
### Transformations

abstract type AbstractDataTransform end

# apply the transform
"""
    transform!(t::AbstractDataTransform, x)

Apply transformation `t` to vector or matrix `x` in place.
"""
transform!(t::AbstractDataTransform, x::AbstractMatrix{<:Real}) =
    transform!(x, t, x)
transform!(t::AbstractDataTransform, x::AbstractVector{<:Real}) =
    (transform!(t, reshape(x, :, 1)); x)

"""
    transform(t::AbstractDataTransform, x)

Return a standardized copy of vector or matrix `x` using transformation `t`.
"""
transform(t::AbstractDataTransform, x::AbstractMatrix{<:Real}) =
    transform!(similar(x), t, x)
transform(t::AbstractDataTransform, x::AbstractVector{<:Real}) =
    vec(transform(t, reshape(x, :, 1)))

# reconstruct the original data from transformed values
"""
    reconstruct!(t::AbstractDataTransform, y)

Perform an in-place reconstruction into an original data scale from a transformed
vector or matrix `y` using transformation `t`.
"""
reconstruct!(t::AbstractDataTransform, y::AbstractMatrix{<:Real}) =
    reconstruct!(y, t, y)
reconstruct!(t::AbstractDataTransform, y::AbstractVector{<:Real}) =
    (reconstruct!(t, reshape(y, :, 1)); y)

"""
    reconstruct(t::AbstractDataTransform, y)

Return a reconstruction of an originally scaled data from a transformed vector
or matrix `y` using transformation `t`.
"""
reconstruct(t::AbstractDataTransform, y::AbstractMatrix{<:Real}) =
    reconstruct!(similar(y), t, y)
reconstruct(t::AbstractDataTransform, y::AbstractVector{<:Real}) =
    vec(reconstruct(t, reshape(y, :, 1)))

"""
Standardization (Z-score transformation)
"""
struct ZScoreTransform{T<:Real} <: AbstractDataTransform
    len::Int
    dims::Int
    mean::Vector{T}
    scale::Vector{T}

    function ZScoreTransform(l::Int, dims::Int, m::Vector{T}, s::Vector{T}) where T
        lenm = length(m)
        lens = length(s)
        lenm == l || lenm == 0 || throw(DimensionMismatch("Inconsistent dimensions."))
        lens == l || lens == 0 || throw(DimensionMismatch("Inconsistent dimensions."))
        new{T}(l, dims, m, s)
    end
end

function Base.getproperty(t::ZScoreTransform, p::Symbol)
    if p === :indim || p === :outdim
        return t.len
    else
        return getfield(t, p)
    end
end

"""
    fit(ZScoreTransform, X; dims=nothing, center=true, scale=true)

Fit standardization parameters to vector or matrix `X`
and return a `ZScoreTransform` transformation object.

# Keyword arguments

* `dims`: if `1` fit standardization parameters in column-wise fashion;
  if `2` fit in row-wise fashion. The default is `nothing`, which is equivalent to `dims=2` with a deprecation warning.

* `center`: if `true` (the default) center data so that its mean is zero.

* `scale`: if `true` (the default) scale the data so that its variance is equal to one.

# Examples

```jldoctest
julia> using StatsBase

julia> X = [0.0 -0.5 0.5; 0.0 1.0 2.0]
2×3 Array{Float64,2}:
 0.0  -0.5  0.5
 0.0   1.0  2.0

julia> dt = fit(ZScoreTransform, X, dims=2)
ZScoreTransform{Float64}(2, 2, [0.0, 1.0], [0.5, 1.0])

julia> StatsBase.transform(dt, X)
2×3 Array{Float64,2}:
  0.0  -1.0  1.0
 -1.0   0.0  1.0
```
"""
function fit(::Type{ZScoreTransform}, X::AbstractMatrix{<:Real};
             dims::Union{Integer,Nothing}=nothing, center::Bool=true, scale::Bool=true)
    if dims === nothing
        Base.depwarn("fit(t, x) is deprecated: use fit(t, x, dims=2) instead", :fit)
        dims = 2
    end
    if dims == 1
        n, l = size(X)
        n >= 2 || error("X must contain at least two rows.")
        m, s = mean_and_std(X, 1)
    elseif dims == 2
        l, n = size(X)
        n >= 2 || error("X must contain at least two columns.")
        m, s = mean_and_std(X, 2)
    else
        throw(DomainError(dims, "fit only accept dims to be 1 or 2."))
    end
    T = eltype(X)
    return ZScoreTransform(l, dims, (center ? vec(m) : zeros(T, 0)),
                                    (scale ? vec(s) : zeros(T, 0)))
end

function fit(::Type{ZScoreTransform}, X::AbstractVector{<:Real};
             dims::Integer=1, center::Bool=true, scale::Bool=true)
    if dims != 1
        throw(DomainError(dims, "fit only accepts dims=1 over a vector. Try fit(t, x, dims=1)."))
    end

    T = eltype(X)
    m, s = mean_and_std(X)
    return ZScoreTransform(1, dims, (center ? [m] : zeros(T, 0)),
                                    (scale ? [s] : zeros(T, 0)))
end

function transform!(y::AbstractMatrix{<:Real}, t::ZScoreTransform, x::AbstractMatrix{<:Real})
    if t.dims == 1
        l = t.len
        size(x,2) == size(y,2) == l || throw(DimensionMismatch("Inconsistent dimensions."))
        n = size(y,1)
        size(x,1) == n || throw(DimensionMismatch("Inconsistent dimensions."))

        m = t.mean
        s = t.scale

        if isempty(m)
            if isempty(s)
                if x !== y
                    copyto!(y, x)
                end
            else
                broadcast!(/, y, x, s')
            end
        else
            if isempty(s)
                broadcast!(-, y, x, m')
            else
                broadcast!((x,m,s)->(x-m)/s, y, x, m', s')
            end
        end
    elseif t.dims == 2
        t_ = ZScoreTransform(t.len, 1, t.mean, t.scale)
        transform!(y', t_, x')
    end
    return y
end

function reconstruct!(x::AbstractMatrix{<:Real}, t::ZScoreTransform, y::AbstractMatrix{<:Real})
    if t.dims == 1
        l = t.len
        size(x,2) == size(y,2) == l || throw(DimensionMismatch("Inconsistent dimensions."))
        n = size(y,1)
        size(x,1) == n || throw(DimensionMismatch("Inconsistent dimensions."))

        m = t.mean
        s = t.scale

        if isempty(m)
            if isempty(s)
                if y !== x
                    copyto!(x, y)
                end
            else
                broadcast!(*, x, y, s')
            end
        else
            if isempty(s)
                broadcast!(+, x, y, m')
            else
                broadcast!((y,m,s)->y*s+m, x, y, m', s')
            end
        end
    elseif t.dims == 2
        t_ = ZScoreTransform(t.len, 1, t.mean, t.scale)
        reconstruct!(x', t_, y')
    end
    return x
end

"""
Unit range normalization
"""
struct UnitRangeTransform{T<:Real}  <: AbstractDataTransform
    len::Int
    dims::Int
    unit::Bool
    min::Vector{T}
    scale::Vector{T}

    function UnitRangeTransform(l::Int, dims::Int, unit::Bool, min::Vector{T}, max::Vector{T}) where {T}
        lenmin = length(min)
        lenmax = length(max)
        lenmin == l || lenmin == 0 || throw(DimensionMismatch("Inconsistent dimensions."))
        lenmax == l || lenmax == 0 || throw(DimensionMismatch("Inconsistent dimensions."))
        new{T}(l, dims, unit, min, max)
    end
end

function Base.getproperty(t::UnitRangeTransform, p::Symbol)
    if p === :indim || p === :outdim
        return t.len
    else
        return getfield(t, p)
    end
end

# fit a unit transform
"""
    fit(UnitRangeTransform, X; dims=nothing, unit=true)

Fit a scaling parameters to vector or matrix `X`
and return a `UnitRangeTransform` transformation object.

# Keyword arguments

* `dims`: if `1` fit standardization parameters in column-wise fashion;
 if `2` fit in row-wise fashion. The default is `nothing`.

* `unit`: if `true` (the default) shift the minimum data to zero.

# Examples

```jldoctest
julia> using StatsBase

julia> X = [0.0 -0.5 0.5; 0.0 1.0 2.0]
2×3 Array{Float64,2}:
 0.0  -0.5  0.5
 0.0   1.0  2.0

julia> dt = fit(UnitRangeTransform, X, dims=2)
UnitRangeTransform{Float64}(2, 2, true, [-0.5, 0.0], [1.0, 0.5])

julia> StatsBase.transform(dt, X)
2×3 Array{Float64,2}:
 0.5  0.0  1.0
 0.0  0.5  1.0
```
"""
function fit(::Type{UnitRangeTransform}, X::AbstractMatrix{<:Real};
             dims::Union{Integer,Nothing}=nothing, unit::Bool=true)
    if dims === nothing
        Base.depwarn("fit(t, x) is deprecated: use fit(t, x, dims=2) instead", :fit)
        dims = 2
    end
    if dims == 1
        l, tmin, tmax = _compute_extrema(X)
    elseif dims == 2
        l, tmin, tmax = _compute_extrema(X')
    else
        throw(DomainError(dims, "fit only accept dims to be 1 or 2."))
    end

    for i = 1:l
        @inbounds tmax[i] = 1 / (tmax[i] - tmin[i])
    end
    return UnitRangeTransform(l, dims, unit, tmin, tmax)
end

function _compute_extrema(X::AbstractMatrix{<:Real})
    n, l = size(X)
    tmin = X[1, :]
    tmax = X[1, :]
    for j = 1:l
        @inbounds for i = 2:n
            if X[i, j] < tmin[j]
                tmin[j] = X[i, j]
            elseif X[i, j] > tmax[j]
                tmax[j] = X[i, j]
            end
        end
    end
    return l, tmin, tmax
end

function fit(::Type{UnitRangeTransform}, X::AbstractVector{<:Real};
             dims::Integer=1, unit::Bool=true)
    if dims != 1
        throw(DomainError(dims, "fit only accept dims=1 over a vector. Try fit(t, x, dims=1)."))
    end

    l, tmin, tmax = _compute_extrema(reshape(X, :, 1))
    tmax = 1 / (tmax - tmin)
    return UnitRangeTransform(1, dims, unit, vec(tmin), vec(tmax))
end

function transform!(y::AbstractMatrix{<:Real}, t::UnitRangeTransform, x::AbstractMatrix{<:Real})
    if t.dims == 1
        l = t.len
        size(x,2) == size(y,2) == l || throw(DimensionMismatch("Inconsistent dimensions."))
        n = size(x,1)
        size(y,1) == n || throw(DimensionMismatch("Inconsistent dimensions."))

        tmin = t.min
        tscale = t.scale

        if t.unit
            broadcast!((x,s,m)->(x-m)*s, y, x, tscale', tmin')
        else
            broadcast!(*, y, x, tscale')
        end
    elseif t.dims == 2
        t_ = UnitRangeTransform(t.len, 1, t.unit, t.min, t.scale)
        transform!(y', t_, x')
    end
    return y
end

function reconstruct!(x::AbstractMatrix{<:Real}, t::UnitRangeTransform, y::AbstractMatrix{<:Real})
    if t.dims == 1
        l = t.len
        size(x,2) == size(y,2) == l || throw(DimensionMismatch("Inconsistent dimensions."))
        n = size(y,1)
        size(x,1) == n || throw(DimensionMismatch("Inconsistent dimensions."))

        tmin = t.min
        tscale = t.scale

        if t.unit
            broadcast!((y,s,m)->y/s+m, x, y, tscale', tmin')
        else
            broadcast!(/, x, y, tscale')
        end
    elseif t.dims == 2
        t_ = UnitRangeTransform(t.len, 1, t.unit, t.min, t.scale)
        reconstruct!(x', t_, y')
    end
    return x
end

"""
    standardize(DT, X; dims=nothing, kwargs...)

 Return a standardized copy of vector or matrix `X` along dimensions `dims`
 using transformation `DT` which is a subtype of `AbstractDataTransform`:

- `ZScoreTransform`
- `UnitRangeTransform`

# Example

```jldoctest
julia> using StatsBase

julia> standardize(ZScoreTransform, [0.0 -0.5 0.5; 0.0 1.0 2.0], dims=2)
2×3 Array{Float64,2}:
  0.0  -1.0  1.0
 -1.0   0.0  1.0

julia> standardize(UnitRangeTransform, [0.0 -0.5 0.5; 0.0 1.0 2.0], dims=2)
2×3 Array{Float64,2}:
 0.5  0.0  1.0
 0.0  0.5  1.0
```
"""
function standardize(::Type{DT}, X::AbstractVecOrMat{<:Real}; kwargs...) where {DT <: AbstractDataTransform}
    return transform(fit(DT, X; kwargs...), X)
end


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

"""
    KLIEP(σ=1.0, b=100)

Kullback-Leibler importance estimation procedure (KLIEP).

## Parameters

* `σ` - Bandwidth of Gaussian kernel (default to `2.0`)
* `b` - Maximum number of radial basis functions (default to `10`)

## References

* Sugiyama et al. 2008. Direct importance estimation for
  covariate shift adaptation.

### Author

* Júlio Hoffimann (julio.hoffimann@gmail.com)
"""
@with_kw struct KLIEP{T,RNG} <: DensityRatioEstimator
  σ::T=2.0
  b::Int=10
  rng::RNG=Random.GLOBAL_RNG
end

default_optlib(dre::Type{<:KLIEP}) = OptimLib

available_optlib(dre::Type{<:KLIEP}) = [OptimLib, ConvexLib]

function _densratio(x_nu, x_de, dre::KLIEP,
                    optlib::Type{<:OptimizationLibrary})
  K_nu, K_de, x_ba = _kliep_consts(x_nu, x_de, dre)
  α = _kliep_coeffs(K_nu, K_de, dre, optlib)
  K_de*α
end

function _densratiofunc(x_nu, x_de, dre::KLIEP,
                        optlib::Type{<:OptimizationLibrary})
  K_nu, K_de, x_ba = _kliep_consts(x_nu, x_de, dre)
  α = _kliep_coeffs(K_nu, K_de, dre, optlib)
  function r(x)
    K = gaussian_gramian([x], x_ba, σ=dre.σ)
    dot(K, α)
  end
end

# constants involved in KLIEP optimization
function _kliep_consts(x_nu, x_de, dre)
  x_ba = sample(dre.rng, x_nu, min(length(x_nu), dre.b), replace=false)
  K_nu = gaussian_gramian(x_nu, x_ba, σ=dre.σ)
  K_de = gaussian_gramian(x_de, x_ba, σ=dre.σ)

  K_nu, K_de, x_ba
end

"""
    _kliep_coeffs(K_nu, K_de, dre, optlib)

Return the coefficients of KLIEP basis expansion.
"""
_kliep_coeffs(K_nu, K_de, dre::KLIEP,
              optlib::Type{<:OptimizationLibrary}) =
  @error "not implemented"


text:
# ## Quiver plots

#md # ![FILE_NAME.png](images/FILE_NAME.png)

n = 7
f(x,y) = 1 / hypot(x, y)
x = repeat(-3:(2*3)/n:3, 1, n) |> vec
y = repeat(-3:(2*3)/n:3, 1, n)' |> vec
vx = f.(x,y) .* cos.(atan.(y,x)) |> vec
vy = f.(x,y) .* sin.(atan.(y,x)) |> vec

quiver(x, y, quiver=(vx, vy), aspect_ratio=:equal)


text:
import LanguageServer: LanguageServerInstance, Document
server = LanguageServerInstance(IOBuffer(), IOBuffer(), false)
server.runlinter = true
LanguageServer.process(LanguageServer.parse(LanguageServer.JSONRPC.Request, JSON.parse(init_request)), server)

function getresult(server)
    str = String(take!(server.pipe_out))
    JSON.parse(str[findfirst("{", str)[1]:end])["result"]["contents"]
end

testtext = """
module testmodule
struct testtype
    a
    b::Float64
    c::Vector{Float64}
end

function testfunction(a, b::Float64, c::testtype)
    return c
end
end
testmodule
"""

server.documents[LanguageServer.URI2("testdoc")] = Document("testdoc", testtext, true)
doc = server.documents[LanguageServer.URI2("testdoc")]
LanguageServer.parse_all(doc, server)

# clear init output
take!(server.pipe_out)

LanguageServer.process(LanguageServer.parse(LanguageServer.JSONRPC.Request, JSON.parse("""{"jsonrpc":"2.0","id":1,"method":"textDocument/hover","params":{"textDocument":{"uri":"testdoc"},"position":{"line":3,"character":11}}}""")), server)
res = getresult(server)
@test startswith(res[1]["value"], "\nFloat64")

LanguageServer.process(LanguageServer.parse(LanguageServer.JSONRPC.Request, JSON.parse("""{"jsonrpc":"2.0","id":1,"method":"textDocument/hover","params":{"textDocument":{"uri":"testdoc"},"position":{"line":8,"character":12}}}""")), server)
res = getresult(server)
@test res[1]["value"] == "testtype"

LanguageServer.process(LanguageServer.parse(LanguageServer.JSONRPC.Request, JSON.parse("""{"jsonrpc":"2.0","id":1,"method":"textDocument/hover","params":{"textDocument":{"uri":"testdoc"},"position":{"line":10,"character":1}}}""")), server)
res = getresult(server)
@test res[1]["value"] == "Closes `module` expression"


text:
module GraphPPLTest

using Test, Documenter, GraphPPL

doctest(GraphPPL)

@testset "GraphPPL" begin

    @testset "Detect ambiguities" begin
        @test length(Test.detect_ambiguities(GraphPPL)) == 0
    end

end

end