text:
# Prim's algorithm on minimum spanning tree

###################################################################
#
#   The type that capsulates the states of Prim's algorithm
#
###################################################################

type PrimStates{V,W,Heap,H}
    parents::Vector{V}
    colormap::Vector{Int}
    weightmap::Vector{W}

    heap::Heap
    hmap::Vector{H}
end

immutable PrimHEntry{V,E,W}
    vertex::V
    edge::E
    weight::W
end

<(e1::PrimHEntry, e2::PrimHEntry) = e1.weight < e2.weight

# create Prim states

function create_prim_states{V,E,W}(g::AbstractGraph{V,E}, ::Type{W})
    n = num_vertices(g)

    parents = Array(V, n)
    colormap = zeros(Int, n)
    weightmap = zeros(W, n)

    heap = mutable_binary_minheap(PrimHEntry{V,E,W})
    hmap = zeros(Int, n)

    PrimStates(parents, colormap, weightmap, heap, hmap)
end

###################################################################
#
#   visitors
#
###################################################################

abstract AbstractPrimVisitor

# invoked when a new vertex is first encountered
discover_vertex!(visitor::AbstractPrimVisitor, v, e, w) = nothing

# invoked when an edge is included into the tree
# it returns whether to continue
include_vertex!(visitor::AbstractPrimVisitor, v, e, w) = true

# invoked when the weight to a vertex is updated
update_vertex!(visitor::AbstractPrimVisitor, v, e, w) = nothing

# invoked when all neighbors of a vertex has been examined
close_vertex!(visitor::AbstractPrimVisitor, v) = nothing

# trivial visitor

type TrivialPrimVisitor <: AbstractPrimVisitor
end

# default visitor

type PrimVisitor{E,W} <: AbstractPrimVisitor
    edges::Vector{E}
    weights::Vector{W}
end

function default_prim_visitor{V,E,W}(g::AbstractGraph{V,E}, ::Type{W})
    edges = Array(E, 0)
    weights = Array(W, 0)
    n = num_vertices(g)
    if n > 1
        sizehint!(edges, n-1)
        sizehint!(weights, n-1)
    end
    PrimVisitor{E,W}(edges, weights)
end

function include_vertex!{V,E,W}(vis::PrimVisitor{E,W}, v::V, e::E, w::W)
    push!(vis.edges, e)
    push!(vis.weights, w)
    true
end

# log visitor

type LogPrimVisitor <: AbstractPrimVisitor
    io::IO
end

function discover_vertex!(visitor::LogPrimVisitor, v, e, w)
    println(visitor.io, "discover vertex $v (edge = $e, w = $w)")
end

function include_vertex!(visitor::LogPrimVisitor, v, e, w)
    println(visitor.io, "include vertex $v (edge = $e, w = $w)")
    true
end

function update_vertex!(visitor::LogPrimVisitor, v, e, w)
    println(visitor.io, "update vertex $v (edge = $e, w = $w)")
end

function close_vertex!(visitor::LogPrimVisitor, v)
    println(visitor.io, "close vertex $v")
end

###################################################################
#
#   core algorithm implementation
#
###################################################################

function process_neighbors!{V,E,W,Heap,H}(
    graph::AbstractGraph{V,E},          # the graph
    edge_weights::AbstractEdgePropertyInspector{W},            # weights associated with edges
    visitor::AbstractPrimVisitor,       # visitor object
    u::V,                               # the vertex whose neighbor to be examined
    state::PrimStates{V,W,Heap,H})      # the states (created)

    parents::Vector{V} = state.parents
    colormap::Vector{Int} = state.colormap
    weightmap::Vector{W} = state.weightmap
    heap::Heap = state.heap
    hmap::Vector{H} = state.hmap

    ew::W = zero(W)

    for e in out_edges(u, graph)
        v = target(e, graph)
        vi = vertex_index(v, graph)
        v_color = colormap[vi]

        if v_color == 0
            ew = edge_property(edge_weights, e, graph)
            colormap[vi] = 1
            weightmap[vi] = ew
            discover_vertex!(visitor, v, e, ew)

            hmap[vi] = push!(heap, PrimHEntry(v, e, ew))

        elseif v_color == 1
            ew = edge_property(edge_weights, e, graph)
            if ew < weightmap[vi]
                weightmap[vi] = ew

                update_vertex!(visitor, v, e, ew)
                update!(heap, hmap[vi], PrimHEntry(v, e, ew))
            end
        end
    end
end

function prim_minimum_spantree!{V,E,W,Heap,H}(
    graph::AbstractGraph{V,E},          # the graph
    edge_weights::AbstractEdgePropertyInspector{W},            # weights associated with edges
    root::V,                            # the root vertex
    visitor::AbstractPrimVisitor,       # visitor object
    state::PrimStates{V,W,Heap,H})        # the states (created)

    @graph_requires graph vertex_map incidence_list

    if is_directed(graph)
        throw(ArgumentError("graph must be undirected."))
    end

    # initialize

    root_idx = vertex_index(root, graph)
    state.parents[root_idx] = root
    state.colormap[root_idx] = 2

    process_neighbors!(graph, edge_weights, visitor, root, state)
    close_vertex!(visitor, root)

    # main loop

    heap::Heap = state.heap

    while !isempty(heap)
        entry = pop!(heap)
        v::V = entry.vertex
        e::E = entry.edge
        w::W = entry.weight

        state.colormap[vertex_index(v, graph)] = 2
        if !include_vertex!(visitor, v, e, w)
            return
        end
        process_neighbors!(graph, edge_weights, visitor, v, state)
        close_vertex!(visitor, v)
    end
end

# convenient functions

function prim_minimum_spantree{V,E,W}(
    graph::AbstractGraph{V,E},
    edge_weight_vec::Vector{W},
    root::V)

    state = create_prim_states(graph, W)
    visitor = default_prim_visitor(graph, W)
    edge_weights = VectorEdgePropertyInspector(edge_weight_vec)
    prim_minimum_spantree!(graph, edge_weights, root, visitor, state)
    return (visitor.edges, visitor.weights)
end

function prim_minimum_spantree{V,E,W}(
    graph::AbstractGraph{V,E},
    edge_weights::AbstractEdgePropertyInspector{W},
    root::V)

    state = create_prim_states(graph, W)
    visitor = default_prim_visitor(graph, W)
    prim_minimum_spantree!(graph, edge_weights, root, visitor, state)
    return (visitor.edges, visitor.weights)
end

function prim_minimum_spantree_withlog{V,E,W}(
    graph::AbstractGraph{V,E},
    edge_weight_vec::Vector{W},
    root::V)

    state = create_prim_states(graph, W)
    visitor = LogPrimVisitor(STDOUT)
    edge_weights = VectorEdgePropertyInspector(edge_weight_vec)
    prim_minimum_spantree!(graph, edge_weights, root, visitor, state)
end


text:
using Iris
using Test

@testset "Dispersions" begin
    @testset "No Dispersion" begin
        @test_nowarn NoDispersion()
        d = NoDispersion()
        @test susceptability(d,1,2,3) ≈ 0
    end
    @testset "Two Level" begin
        @test_nowarn TwoLevelSystem()
        d = TwoLevelSystem()
        @test_nowarn susceptability(d,3)
        @test susceptability(d,3) ≈ 0
        ω = 3
        ωs = 2.5 .+ 1rand(3)
        @testset "Scalar" begin
            r = rand(ComplexF64,20,length(ωs))
            pos = Vector{Point{1}}(undef,20)
            for i ∈ 1:20 pos[i] = Point(randn(Float64,1)...) end
            e = ScalarField(pos,r)
            @test_nowarn susceptability(d,ω,ωs,e)
            d.D₀ = .1
            d.ωa = ω
            @test_nowarn susceptability(d,ω,ωs,e)
        end
        @testset "Vector" begin
            N = 1 # dimension
            M = 3 # number of components
            n = 25 # number of sites
            nmodes = 3
            r = rand(ComplexF64,M*n,nmodes)
            pos = Vector{Point{N}}(undef,n)
            for i ∈ 1:n pos[i] = Point(randn(Float64,N)...) end
            e = ElectricField(pos,r)
            @test_nowarn susceptability(d,ω,ωs,e)
            d.D₀ = .1
            d.ωa = ω
            @test_nowarn susceptability(d,ω,ωs,e)
        end
    end
    # @testset "Kerr" begin
    #     d = Kerr()
    #     @test susceptability(d,1,2,3) ≈ 0.0
    # end
end

ext: jl


text:
# convenience function to calculate the mean-squared error
function mse(arr1::AbstractArray, arr2::AbstractArray)
    @assert length(arr1) == length(arr2)
    N = length(arr1)
    err = 0.0
    for i in 1:N
        err += (arr2[i] - arr1[i])^2
    end
    err /= N
end


text:
struct IJuliaCell

    window::ViewerWindow
    embed::Bool

    IJuliaCell(window, embed=false) = new(window, embed)
end

IJuliaCell(v::Visualizer) = IJuliaCell(v.window)

url(c::IJuliaCell) = url(c.window)

function Base.show(io::IO, ::MIME"text/html", frame::IJuliaCell)
    if frame.embed
        show_embed(io, frame)
    else
        show_inline(io, frame)
    end
end

const iframe_attrs = "height=\"100%\" width=\"100%\" style=\"min-height: 500px;\""

function show_inline(io::IO, frame::IJuliaCell)
    print(io, """

""")
end

srcdoc_escape(x) = replace(replace(x, "&", "&"), "\"", """)

function show_embed(io::IO, frame::IJuliaCell)
    id = Base.Random.uuid1()
    print(io, """

    """)
end

# struct Snapshot
#   json::String

#     Snapshot(fname::AbstractString) = new(open(readstring, fname))
#     Snapshot(io::IO) = new(readstring(io))
# end

# function Base.show(io::IO, ::MIME"text/html", snap::Snapshot)
#   content = readstring(open(joinpath(viewer_root, "build", "inline.html")))
#   # TODO: there has to be a better way than doing a replace() on the html.
#   script = """
#
#
#   """
#   html = replace(content, "", script)
#   print(io, """
#
#
#   """)
# end

# function save(fname::String, snap::Snapshot)
#     content = readstring(open(joinpath(viewer_root, "build", "inline.html")))
#     # TODO: there has to be a better way than doing a replace() on the html.
#     script = """
#
#
#     """
#     html = replace(content, "", script)
#     open(fname, "w") do file
#         write(file, html)
#     end
# end


text:
#-------------------------------------------------------------------------------------------------
#--
#-- Ususal libraries
#--
using CSV, Dates
using DataFrames, DataFramesMeta

#-------------------------------------------------------------------------------------------------
#--
#-- DATA
#--
# Update latest data from Neherlab
function updateData()

    ### COUNTRY CODES
    filename = "https://github.com/neherlab/covid19_scenarios_data/raw/master/country_codes.csv"
    download(filename, "data/country_codes.csv")

    ### LIST OF CASES
    filename = "https://github.com/neherlab/covid19_scenarios_data/raw/master/case-counts/World.tsv"
    download(filename, "data/cases_world.tsv")

    ### ICU beds
    filename = "https://github.com/neherlab/covid19_scenarios_data/raw/master/hospital-data/ICU_capacity.tsv"
    download(filename, "data/ICU_capacity.tsv")

    ### Hospital beds
    filename = "https://github.com/neherlab/covid19_scenarios_data/raw/master/hospital-data/hospital_capacity.csv"
    download(filename, "data/hospital_capacity.csv")
end

function loadData()
    ### Population Data
    popData = CSV.read("data/populationData.tsv", delim = "\t")

    ### COUNTRY CODES
    country_codes = select(CSV.read("data/country_codes.csv"), :name, Symbol("alpha-3"))

    ### LIST OF CASES
    cases = DataFrame(CSV.read("data/cases_world.tsv", delim = "\t", header = 4))

    # Add a time column in the same format as the other dataframes
    cases = hcat(DataFrame(t = date2days.(cases[:, :time])), cases)

    # Remove any row with no recorded death
    cases = cases[cases.deaths.>0, :]

    ### ICU beds
    ICU_capacity = select(CSV.read("data/ICU_capacity.tsv", delim = "\t"), :country, :CriticalCare)

    ### Hospital beds
    hospital_capacity = select(
        CSV.read("data/hospital_capacity.csv", types = Dict(:COUNTRY => String), limit = 1267),
        :COUNTRY,
        :YEAR,
        :VALUE,
    )
    hospital_capacity = @where(hospital_capacity, Not(ismissing.(:COUNTRY)))

    ### Age pyramids
    age_distribution = CSV.read("data/country_age_distribution.csv")

    return Dict(
        :codes => country_codes,
        :popData => popData,
        :cases => cases,
        :ICU => ICU_capacity,
        :Beds => hospital_capacity,
        :age_distribution => age_distribution,
    )
end

# Let's get the Neherlab repo's data (the files can be updated with `updateData()`).
COUNTRY_DB = loadData();

#-------------------------------------------------------------------------------------------------
#--
#-- Create a dictionary populated with the information of a given country
#--
function populateCountryDate(country, hemisphere)

    # This is the start of a country specific data structure. A dictionary is good enough for that purpose.
    countryData = Dict()

    countryData = Dict([(:name, country), (:R₀, BaseR₀), (:hemisphere, hemisphere)])

    countryData[:peak] = peakDate[countryData[:hemisphere]]
    countryData[:ϵ] = ϵ[countryData[:hemisphere]]
    countryData[:mitigation] = DEFAULT_MITIGATION

    # Clean up and extract country-specific information
    country_codes = @where(COUNTRY_DB[:codes], occursin.(country, :name))
    countryShort = country_codes[:, Symbol("alpha-3")][1]

    cases = @where(COUNTRY_DB[:cases], occursin.(country, :location))
    sort!(cases, :time)

    ICU_capacity = try
        @where(COUNTRY_DB[:popData], occursin.(country, :name))[!, :ICUBeds][1]
    catch
        0
    end
    ICU_capacity = convert(Float64, ICU_capacity)

    hospital_capacity = try
        @where(COUNTRY_DB[:popData], occursin.(country, :name))[!, :hospitalBeds][1]
    catch
        0
    end
    hospital_capacity = convert(Float64, hospital_capacity)

    age_distribution = @where(COUNTRY_DB[:age_distribution], occursin.(country, :_key))[!, 2:10]

    countryData[:country_code] = countryShort
    countryData[:cases] = cases
    countryData[:ICU_capacity] = ICU_capacity
    countryData[:hospital_capacity] = hospital_capacity
    countryData[:age_distribution] = age_distribution

    p, r = createDefaultParameters()
    countryData[:params] = p
    countryData[:range] = r

    return countryData
end


text:
using LinearAlgebra
BLAS.set_num_threads(1)

using Zygote, Flux, BenchmarkTools, Test, DataFrames, GLM, Statistics, StaticArrays
using Plots

to_smatrix(x) = x
function to_smatrix(x::Matrix)
    return SMatrix{size(x)...}(x)
end

function build_model(num_layers::Int, size::Int)
    layers = Any[]
    for layer_idx in 1:num_layers
        push!(layers, LSTM(size, size))
    end
    return mapleaves(to_smatrix, mapleaves(Flux.data, Chain(layers...)))
end

@info "Reporting (num_layers, seq_len, feature_size, batch_size, time)"
function sweep_batchsizes(batch_sizes, num_layers = 1, seq_len = 4, feature_size = 4)
    timings = Float64[]
    for batch_size in batch_sizes
        model = build_model(num_layers, feature_size)
        x = SMatrix{feature_size,batch_size}(randn(feature_size, batch_size))
        test(m, x) = gradient(m, x) do m, x
            x_t = x
            for idx in 1:seq_len
                x_t = m(x_t)
            end
            return sum(x_t)
        end
        push!(timings, minimum(@benchmark $test($model, $x)).time)
        @info num_layers, seq_len, feature_size, batch_size, timings[end]
    end
    return timings
end

batch_sizes = 1:8
layer_sizes = 1:3
timings = Dict()
for num_layers in layer_sizes,
    seq_len in (4,),
    feature_size in (4,)

    timings[(num_layers, seq_len, feature_size, batch_sizes)] = sweep_batchsizes(batch_sizes, num_layers, seq_len, feature_size)
end

# These gotten via count_ops/count_lstm_ops.jl
num_ops = Dict(
    1 => 255,
    2 => 491,
    3 => 727,
    4 => 963,
)

overhead_estimates = Dict()
workload_estimates = Dict()
for k in keys(timings)
    num_layers, seq_len, feature_size, batch_sizes = k
    data = DataFrame(X=batch_sizes, Y=timings[k])
    ols = lm(@formula(Y ~ X), data)
    overhead_estimates[k] = coeftable(ols).cols[1][1]/num_ops[num_layers]
    workload_estimates[k] = coeftable(ols).cols[1][2]
end

gr()
ENV["GKSwstype"]="100"
p = plot()
for key in sort(collect(keys(timings)))
    num_layers, seq_len, feature_size, batch_sizes = key
    plot!(p, batch_sizes, timings[key]; title="LSTM Runtime", ylabel="Absolute Runtime (ns)", xlabel="Batch Size", label="$(num_layers) layers")
    extrapolate(x) = overhead_estimates[key]*num_ops[num_layers] + workload_estimates[key]*x
    x_points = [0, batch_sizes...]
    plot!(p, x_points, extrapolate.(x_points); style=:dashdot, label="$(num_layers) layers extrapolation")
end
savefig(p, joinpath(@__DIR__, "lstm_runtime.pdf"))

max_stars_repo_name: MikeInnes/zygote-paper

max_issues_repo_name: MikeInnes/zygote-paper

max_forks_repo_name: MikeInnes/zygote-paper


text:
# Reexport GridTools as they are useful outside too
@reexport module GridTools

using LinearAlgebra, SparseArrays

#===================================
Off-diaongal matrices
===================================#

"""
    buildIbelow(wet3D, iwet)

Build the shifted-diagonal sparse matrix of the indices of below neighbours.

Ibelow[i,j] = 1 if the box represented by the linear index i
lies directly below the box represented by the linear index j.
"""
function buildIbelow(wet3D, iwet)
    nlat, nlon, ndepth = size(wet3D)
    n = nlon * nlat * (ndepth + 1)
    In = sparse(I, n, n)
    idx = zeros(Int, nlat, nlon, ndepth + 1)
    idx[:] = 1:n
    idx .= idx[:, :, [2:ndepth + 1; 1]]      # downward shift
    return In[idx[:], :][iwet, iwet]
end
export buildIbelow

"""
    buildIabove(wet3D, iwet)

Build the shifted-diagonal sparse matrix of the indices of above neighbours.

Iabove[i,j] = 1 if the box represented by the linear index i
lies directly above the box represented by the linear index j.
"""
buildIabove(wet3D, iwet) = copy(transpose(buildIbelow(wet3D, iwet)))
export buildIabove

#===================================
Functions returning useful arrays,
vectors, and constants
===================================#

"""
    array_of_volumes(grid)

Returns the 3D array of volumes of grid boxes.
"""
array_of_volumes(grid) = grid.volume_3D
export array_of_volumes

"""
    indices_of_wet_boxes(wet3D)

Returns the vector of the indices of wet grid boxes.
"""
indices_of_wet_boxes(wet3D) = (LinearIndices(wet3D))[findall(!iszero, wet3D)]
export indices_of_wet_boxes

"""
    number_of_wet_boxes(wet3D)

Returns the number of wet grid boxes.
"""
number_of_wet_boxes(wet3D) = length(indices_of_wet_boxes(wet3D))
export number_of_wet_boxes

"""
    vector_of_volumes(wet3D, grid)

Returns the vector of volumes of wet boxes.
"""
vector_of_volumes(wet3D, grid) = array_of_volumes(grid)[indices_of_wet_boxes(wet3D)]
export vector_of_volumes

"""
    vector_of_depths(wet3D, grid)

Returns the vector of depths of the center of wet boxes.
"""
vector_of_depths(wet3D, grid) = grid.depth_3D[indices_of_wet_boxes(wet3D)]
export vector_of_depths

"""
    vector_of_depths(wet3D, grid)

Returns the vector of depths of the top of wet boxes.
"""
vector_of_top_depths(wet3D, grid) = grid.depth_top_3D[indices_of_wet_boxes(wet3D)]
export vector_of_top_depths

#===================================
Functions return weighted norms
===================================#

"""
    weighted_norm²(x, w)

Returns the square of the weighted norm of `x` using weights `w`.
"""
weighted_norm²(x, w) = transpose(x) * Diagonal(w) * x

"""
    weighted_norm(x, v)

Returns the weighted norm of `x` using weights `w`.
"""
weighted_norm(x, w) = sqrt(weighted_norm²(x, w))

"""
    weighted_mean(x, w)

Returns the weighted mean of `x` using weights `w`.
"""
weighted_mean(x, w) = transpose(w) * x / sum(w)

"""
    vmean(x, w, I)

Returns the weighted mean of `x` using weights `w`, but only over indices `I`.
This is useful to get the observed mean.
(Because there are some grid boxes without observations.)
"""
weighted_mean(x, w, I) = transpose(w[I]) * x[I] / sum(w[I])
export weighted_norm²

#===================================
1D Vector <-> 3D array conversions
===================================#

"""
    rearrange_into_3Darray(x, wet3D)

Returns a 3D array of `x` rearranged to the `true` entries of `wet3D`.
Entries where `wet3D` is `false` are filled with `NaN`s.
"""
function rearrange_into_3Darray(x, wet3D)
    iwet = indices_of_wet_boxes(wet3D)
    x3d = fill(NaN, size(wet3D))
    x3d[iwet] .= x
    return x3d
end

"""
    rearrange_into_1Dvector(x3d, wet3D)

Returns a 1D vector of `x3d` from the linear indices where `wet3D` is `true`.
"""
function rearrange_into_1Dvector(x3d, wet3D)
    iwet = indices_of_wet_boxes(wet3D)
    return x3d[iwet]
end
export rearrange_into_3Darray, rearrange_into_1Dvector

#=============================================
unpacking of multi-tracers
=============================================#

state_to_tracers(x, nb, nt) = ntuple(i -> state_to_tracer(x, nb, nt, i), nt)
state_to_tracer(x, nb, nt, i) = x[tracer_indices(nb, nt, i)]
tracer_indices(nb, nt, i) = (i-1)*nb+1 : i*nb
tracers_to_state(xs) = reduce(vcat, xs)
export state_to_tracers, state_to_tracer, tracers_to_state, tracer_indices

end


text:
using Distributed, Logging, Random, Test

function safe_addprocs(n)
    try
        length(addprocs(n))
    catch
        @warn "problem calling addprocs, nworkers=$(nworkers())"
    end
end

@testset "pmap, stable cluster test" begin
    safe_addprocs(5)
    @everywhere using Distributed, Schedulers
    s = randstring(6)
    @everywhere function foo1(tsk, s)
        write(joinpath(tempdir(), "task-$s-$tsk.txt"), "$tsk,$(myid())")
        sleep(10)
    end
    epmap(foo1, 1:10, s)

    h = Dict()
    for w in workers()
        h[w] = 0
    end

    @test nworkers() == 5
    rmprocs(workers())

    for tsk = 1:10
        r = read(joinpath(tempdir(), "task-$s-$tsk.txt"), String)
        r_tsk, r_pid = split(r, ",")
        h[parse(Int,r_pid)] += 1
        @test r_tsk == "$tsk"
        rm(joinpath(tempdir(), "task-$s-$tsk.txt"))
    end

    for (key,value) in h
        @test value ∈ 0:4
    end
end

@testset "pmap, growing cluster test" begin
    safe_addprocs(2)
    @everywhere using Distributed, Schedulers
    s = randstring(6)
    @everywhere function foo2(tsk, s)
        write(joinpath(tempdir(), "task-$s-$tsk.txt"), "$tsk, $(myid())")
        sleep(10)
    end

    epmap(foo2, 1:100, s; epmap_maxworkers=10, epmap_minworkers=10)

    h = Dict()
    for w in workers()
        h[w] = 0
    end

    @test nworkers() == 10
    rmprocs(workers())

    for tsk = 1:100
        r = read(joinpath(tempdir(), "task-$s-$tsk.txt"), String)
        r_tsk, r_pid = split(r, ",")
        h[parse(Int,r_pid)] += 1
        @test r_tsk == "$tsk"
        rm(joinpath(tempdir(), "task-$s-$tsk.txt"))
    end

    for (key,value) in h
        @test value ∈ 5:15
    end
end

@testset "pmap, elastic cluster with faults" begin
    safe_addprocs(10)
    wrkrs = workers()
    @everywhere using Distributed, Schedulers
    s = randstring(6)
    @everywhere function foo3(tsk, s)
        write(joinpath(tempdir(), "task-$s-$tsk.txt"), "$tsk, $(myid())")
        sleep(10)
    end
    tsk = @async epmap(foo3, 1:100, s; epmap_maxworkers=10)

    sleep(15)
    faulty_pids = workers()[randperm(length(workers()))[1:2]]
    rmprocs(faulty_pids)

    wait(tsk)
    @test nworkers() == 10
    _wrkrs = workers()
    @test wrkrs != _wrkrs

    h = Dict()

    rmprocs(workers())

    for tsk = 1:100
        r = read(joinpath(tempdir(), "task-$s-$tsk.txt"), String)
        r_tsk, r_pid = split(r, ",")
        if haskey(h, parse(Int, r_pid))
            h[parse(Int,r_pid)] += 1
        else
            h[parse(Int,r_pid)] = 1
        end
        @test r_tsk == "$tsk"
        rm(joinpath(tempdir(), "task-$s-$tsk.txt"))
    end

    for (key,value) in h
        if key ∈ _wrkrs
            @test value ∈ 5:15
        end
    end
end

@testset "pmap with shrinking cluster" begin
    safe_addprocs(10)
    @everywhere using Distributed, Schedulers
    s = randstring(6)
    @everywhere function foo4(tsk, s)
        write(joinpath(tempdir(), "task-$s-$tsk.txt"), "$tsk, $(myid())")
        sleep(10)
    end

    h = Dict()
    for w in workers()
        h[w] = 0
    end

    epmap(foo4, 1:105, s; epmap_maxworkers=10, epmap_minworkers=4)

    @test nworkers() == 4
    rmprocs(workers())

    for tsk = 1:105
        r = read(joinpath(tempdir(), "task-$s-$tsk.txt"), String)
        r_tsk, r_pid = split(r, ",")
        h[parse(Int,r_pid)] += 1
        @test r_tsk == "$tsk"
        rm(joinpath(tempdir(), "task-$s-$tsk.txt"))
    end

    for (key,value) in h
        @test value ∈ 5:15
    end
end

@testset "pmap with interactive growing cluster" begin
    safe_addprocs(2)
    @everywhere using Distributed, Schedulers
    s = randstring(6)
    @everywhere function foo5(tsk, s)
        write(joinpath(tempdir(), "task-$s-$tsk.txt"), "$tsk, $(myid())")
        sleep(10)
    end

    _nworkers = 5

    tsk = @async epmap(foo5, 1:100, s; epmap_maxworkers=()->_nworkers, epmap_minworkers=()->_nworkers)

    sleep(15)
    _nworkers = 10

    wait(tsk)
    @test nworkers() == 10

    h = Dict()
    for w in workers()
        h[w] = 0
    end

    rmprocs(workers())

    for tsk = 1:100
        r = read(joinpath(tempdir(), "task-$s-$tsk.txt"), String)
        r_tsk, r_pid = split(r, ",")
        h[parse(Int,r_pid)] += 1
        @test r_tsk == "$tsk"
        rm(joinpath(tempdir(), "task-$s-$tsk.txt"))
    end

    for (key,value) in h
        @test value ∈ 5:15
    end
end

@testset "pmapreduce, stable cluster test" begin
    safe_addprocs(5)
    @everywhere using Distributed, Schedulers, Random
    @everywhere function foo6(x, tsk, a; b)
        fetch(x)::Vector{Float32} .+= a*b*tsk
        nothing
    end

    tmpdir = mktempdir(;cleanup=false)

    a,b = 2,3
    x = epmapreduce!(zeros(Float32,10), foo6, 1:100, a; b=b, epmap_scratch=tmpdir)
    rmprocs(workers())
    @test x ≈ sum(a*b*[1:100;]) * ones(10)
    @test mapreduce(file->startswith("checkpoint", file), +, ["x";readdir(tmpdir)]) == 0
    rm(tmpdir; recursive=true, force=true)
end

@testset "pmapreduce, cluster with ProcessExitedException during tasks" begin
    # TODO, how do we excersie the differnt fault mechanisms in the task loop
    #       1. fault during a f eval
    #       2. fault during checkpoint
    #       3. fault when removing old checkpoint
    #       4. fault when reducing-in an ophan
    #       5. fault when checkpoiting the reduced-in result
    #       6. fault when removing old orphaned checkpoints
    safe_addprocs(5)
    wrkrs = workers()
    @everywhere using Distributed, Schedulers, Random
    s = randstring(6)
    @everywhere function foo7(x, tsk, a, b)
        fetch(x)::Vector{Float32} .+= a*b*tsk
        sleep(1)
        nothing
    end

    a,b = 2,3

    tmpdir = mktempdir(;cleanup=false)

    tsk = @async epmapreduce!(zeros(Float32,10), foo7, 1:100, a, b; epmap_maxworkers=5, epmap_scratch=tmpdir)

    sleep(10)
    rmprocs(workers()[randperm(nworkers())[1]])

    x = fetch(tsk)
    @test nworkers() == 5
    _wrkrs = workers()
    @test wrkrs != _wrkrs

    rmprocs(workers())

    @test x ≈ sum(a*b*[1:100;]) * ones(10)
    @test mapreduce(file->startswith(file, "checkpoint"), +, ["x";readdir(tmpdir)]) == 0
    rm(tmpdir; recursive=true, force=true)
end

@testset "pmapreduce, cluster with RemoteException during tasks" begin
    safe_addprocs(5)
    @everywhere using Distributed, Schedulers, Random
    s = randstring(6)
    @everywhere function foo8(x, tsk, a, b, toggle, fault_id)
        _toggle = fetch(toggle)
        if myid() == fault_id && _toggle[1]
            _toggle[1] = false
            error("throwing an error")
        end

        fetch(x)::Vector{Float32} .+= a*b*tsk
        sleep(1)
        nothing
    end

    a,b = 2,3

    tmpdir = mktempdir(;cleanup=false)

    _pid = workers()[randperm(nworkers())[1]]
    toggle = remotecall_wait(()->[true], _pid)
    x = epmapreduce!(zeros(Float32,10), foo8, 1:10, a, b, toggle, _pid; epmap_maxworkers=5, epmap_scratch=tmpdir, epmap_retries=1, epmap_maxerrors=Inf)

    @test nworkers() == 5
    rmprocs(workers())

    @test x ≈ sum(a*b*[1:10;]) * ones(10)
    @test mapreduce(file->startswith(file, "checkpoint"), +, ["x";readdir(tmpdir)]) == 0
    rm(tmpdir; recursive=true, force=true)
end

@testset "pmapreduce, cluster with RemoteException during tasks, and max errors triggered" begin
    safe_addprocs(5)
    @everywhere using Distributed, Schedulers, Random
    s = randstring(6)
    @everywhere function foo9(x, tsk, a, b, toggle, fault_id)
        _toggle = fetch(toggle)
        if myid() == fault_id && _toggle[1]
            _toggle[1] = false
            error("throwing an error")
        end

        fetch(x)::Vector{Float32} .+= a*b*tsk
        sleep(1)
        nothing
    end

    a,b = 2,3

    tmpdir = mktempdir(;cleanup=false)

    _pid = workers()[randperm(nworkers())[1]]
    toggle = remotecall_wait(()->[true], _pid)
    @test_throws Exception epmapreduce!(zeros(Float32,10), foo9, 1:10, a, b, toggle, _pid; epmap_maxworkers=5, epmap_scratch=tmpdir, epmap_retries=1, epmap_maxerrors=1)
    rmprocs(workers())

    rm(tmpdir; recursive=true, force=true)
end

@testset "pmapreduce, cluster with ProcessExitedException during reduce" begin ### whoop!!!
    # TODO, how do we excersie the differnt fault mechanisms in the task loop
    #       1. fault during reduce
    #       2. fault during removal of checkpoint1
    #       3. fault during removal of checkpoint2
    safe_addprocs(5)
    __workers = workers()
    @everywhere using Distributed, Schedulers, Random
    s = randstring(6)
    @everywhere function foo10(x, tsk, a, b)
        fetch(x)::Vector{Float32} .+= a*b*tsk
        sleep(1)
        nothing
    end

    a,b = 2,3
    id = randstring(6)

    tmpdir = mktempdir(;cleanup=false)

    epmap_zeros = ()->zeros(Float32,10)

    result = epmap_zeros()

    journal = Schedulers.journal_init(1:100, tsks->nothing; reduce=true)

    eloop = Schedulers.ElasticLoop(;
        epmap_init = Schedulers.epmap_default_init,
        epmap_addprocs = Schedulers.epmap_default_addprocs,
        epmap_quantum = ()->32,
        epmap_minworkers = nworkers(),
        epmap_maxworkers = nworkers(),
        epmap_nworkers = nworkers,
        epmap_usemaster = false,
        tasks = 1:100,
        exit_on_empty = false)

    _elastic_loop = @async Schedulers.loop(eloop)

    checkpoints = Schedulers.epmapreduce_map(foo10, result, eloop, journal, a, b;
        epmapreduce_id = id,
        epmap_reducer! = Schedulers.default_reducer!,
        epmap_zeros = epmap_zeros,
        epmap_preempted = Schedulers.epmap_default_preempted,
        epmap_scratches = [tmpdir],
        epmap_reporttasks = true,
        epmap_maxerrors = Inf,
        epmap_retries = 0,
        epmap_keepcheckpoints = false,
        epmap_journal_task_callback = tsk->nothing)

    empty!(eloop, [1:length(checkpoints)-1;])
    eloop.exit_on_empty = true

    tsk = @async Schedulers.epmapreduce_reduce!(result, checkpoints, eloop, journal;
        epmapreduce_id = id,
        epmap_reducer! = Schedulers.default_reducer!,
        epmap_preempted = Schedulers.epmap_default_preempted,
        epmap_scratches = [tmpdir],
        epmap_reporttasks = true,
        epmap_maxerrors = Inf,
        epmap_retries = 0,
        epmap_keepcheckpoints = false)

    rmprocs(workers()[randperm(nworkers())[1]])

    fetch(_elastic_loop)
    _workers = workers()
    @test _workers != __workers
    1 ∈ _workers && popfirst!(_workers)
    rmprocs(_workers[1:(length(_workers) - eloop.epmap_minworkers())])

    x = fetch(tsk)

    Schedulers.journal_final(journal)

    rmprocs(workers())

    @test x ≈ sum(a*b*[1:100;])*ones(10)
    @test mapreduce(file->startswith(file, "checkpoint"), +, ["x";readdir(tmpdir)]) == 0
    rm(tmpdir; recursive=true, force=true)
end

@testset "pmapreduce, cluster with RemoteException during reduce" begin
    # TODO, how do we excersie the different fault mechanisms in the task loop
    #       1. fault during reduce
    #       2. fault during removal of checkpoint1
    #       3. fault during removal of checkpoint2

    safe_addprocs(5)
    @everywhere using Distributed, Schedulers, Random
    s = randstring(6)
    @everywhere function foo11(x, tsk, a, b)
        fetch(x)::Vector{Float32} .+= a*b*tsk
        sleep(1)
        nothing
    end

    a,b = 2,3
    id = randstring(6)

    tmpdir = mktempdir(;cleanup=false)

    _pid = workers()[randperm(nworkers())[1]]
    toggle = remotecall_wait(()->[true], _pid)

    @everywhere function myreducer(x, y, toggle, _pid)
        _toggle = fetch(toggle)
        if _toggle[1] && _pid == myid()
            _toggle[1] = false
            error("this is an error")
        end
        x .+= y
        nothing
    end

    x = epmapreduce!(zeros(Float32,10), foo11, 1:20, a, b;
        epmap_reducer! = (x,y)->myreducer(x,y,toggle,_pid),
        epmap_maxworkers = 5,
        epmap_scratch = tmpdir,
        epmap_retries = 1,
        epmap_maxerrors = Inf)

    rmprocs(workers())

    @test x ≈ sum(a*b*[1:20;])*ones(10)
    @test mapreduce(file->startswith(file, "checkpoint"), +, ["x";readdir(tmpdir)]) == 0
    rm(tmpdir; recursive=true, force=true)
end

@testset "pmapreduce, growing cluster test" begin
    safe_addprocs(5)
    @everywhere using Distributed, Schedulers, Random
    s = randstring(6)
    @everywhere function foo12(x, tsk, a, b)
        fetch(x)::Vector{Float32} .+= a*b*tsk
        sleep(5)
        nothing
    end

    a,b = 2,3

    tmpdir = mktempdir(;cleanup=false)

    x = epmapreduce!(zeros(Float32,10), foo12, 1:100, a, b; epmap_maxworkers=10, epmap_scratch=tmpdir)
    rmprocs(workers())
    @test x ≈ sum(a*b*[1:100;]) * ones(10)

    @test mapreduce(file->startswith("checkpoint", file), +, ["x";readdir(tmpdir)]) == 0
    rm(tmpdir; recursive=true, force=true)
end

@testset "pmapreduce, interactive growing cluster test" begin
    safe_addprocs(5)
    @everywhere using Distributed, Schedulers, Random
    s = randstring(6)
    @everywhere function foo13(x, tsk, a, b)
        fetch(x)::Vector{Float32} .+= a*b*tsk
        sleep(5)
        nothing
    end

    a,b = 2,3

    tmpdir = mktempdir(;cleanup=false)

    _nworkers = 5

    local x
    tsk = @async epmapreduce!(zeros(Float32,10), foo13, 1:100, a, b; epmap_maxworkers=()->_nworkers, epmap_scratch=tmpdir)

    sleep(20)
    _nworkers = 10
    sleep(20)
    @test nworkers() == 10

    x = fetch(tsk)

    rmprocs(workers())
    @test x ≈ sum(a*b*[1:100;]) * ones(10)

    @test mapreduce(file->startswith("checkpoint", file), +, ["x";readdir(tmpdir)]) == 0
    rm(tmpdir; recursive=true, force=true)
end

@testset "pmapreduce, structured data test" begin
    safe_addprocs(5)
    @everywhere using Distributed, Schedulers, Random
    s = randstring(6)
    @everywhere function foo14(x, tsk, a, b)
        fetch(x).y::Vector{Float32} .+= a*tsk
        fetch(x).z::Vector{Float32} .+= b*tsk
        sleep(1)
        nothing
    end

    a,b = 2,3

    tmpdir = mktempdir(;cleanup=false)

    my_zeros() = (y=zeros(Float32,10),z=zeros(Float32,10))
    x = my_zeros()

    epmapreduce!(x, foo14, 1:100, a, b;
        epmap_maxworkers = 10,
        epmap_scratch = tmpdir,
        epmap_zeros = my_zeros,
        epmap_reducer! = (x,y)->(x.y .+= y.y; x.z .+= y.z; nothing))

    rmprocs(workers())

    @test x.y ≈ sum([a*tsk for tsk in 1:100])*ones(Float32,10)
    @test x.z ≈ sum([b*tsk for tsk in 1:100])*ones(Float32,10)
end

@testset "pmapreduce, multiple scratch locations" begin
    safe_addprocs(5)
    @everywhere using Distributed, Schedulers, Random
    @everywhere function foo15(x, tsk, a; b)
        fetch(x)::Vector{Float32} .+= a*b*tsk
        nothing
    end

    tmpdirs = [mktempdir(;cleanup=false) for i=1:3]

    a,b = 2,3
    x = epmapreduce!(zeros(Float32,10), foo15, 1:100, a; b=b, epmap_scratch=tmpdirs, epmap_keepcheckpoints=true)

    ncheckpoints = [length(readdir(tmpdir)) for tmpdir in tmpdirs]
    ncheckpoints_average = sum(ncheckpoints) / 3
    for i = 1:3
        @test ncheckpoints[i] > 0
        @test (ncheckpoints[i] - ncheckpoints_average) < .1*ncheckpoints_average
    end

    rmprocs(workers())
    @test x ≈ sum(a*b*[1:100;]) * ones(10)
    rm.(tmpdirs; recursive=true, force=true)
end

@testset "logerror" begin
    try
        notafunction()
    catch e
        io = IOBuffer()
        with_logger(ConsoleLogger(io, Logging.Info)) do
            Schedulers.logerror(e)
        end
        s = String(take!(io))
        @test contains(s, "notafunction")
    end
end


text:
using StagedFilters
using Test

using PyCall, BenchmarkTools;

N = 1_000_000

data = convert.(Float64,collect(range(1,N,length=N)));
smoothed = zeros(eltype(data),length(data)); # <--- wholesome, type stable code.
savgol = pyimport("scipy.signal")."savgol_filter";
x = PyObject(data);

@info "Julia f32x$N"
@btime smooth!(SavitzkyGolayFilter{2,2}, $data, $smoothed);
@info "SciPy f32x$N"
@btime pycall($savgol, PyObject, $x,5,2,mode="wrap");

data = convert.(Float32,collect(range(1,N,length=N)));
smoothed = zeros(eltype(data),length(data)); # <--- wholesome, type stable$
savgol = pyimport("scipy.signal")."savgol_filter";
x = PyObject(data);

@info "Julia f64x$N"
@btime smooth!(SavitzkyGolayFilter{2,2}, $data, $smoothed);
@info "SciPy f64x$N"
@btime pycall($savgol, PyObject, $x,5,2,mode="wrap");
nothing;


text:
@testset "whittledata" begin
    @test_throws ArgumentError WhittleLikelihoodInference.WhittleData(CorrelatedOU, ones(10,3), 1.0)
    @test_throws ArgumentError WhittleLikelihoodInference.DebiasedWhittleData(CorrelatedOU, ones(10,3), 1.0)
    @test_throws ArgumentError WhittleLikelihoodInference.WhittleData(CorrelatedOU, ones(10,2), 1.0, taper = ones(10,2))
    @test_throws ArgumentError WhittleLikelihoodInference.DebiasedWhittleData(CorrelatedOU, ones(10,2), 1.0, taper = ones(10,2))
    @test_throws ArgumentError WhittleLikelihoodInference.WhittleData(CorrelatedOU, ones(10,2), 1.0, taper = ones(ComplexF64,10))
    @test_throws ArgumentError WhittleLikelihoodInference.DebiasedWhittleData(CorrelatedOU, ones(10,2), 1.0, taper = ones(ComplexF64,10))
end


text:
# ----------------------
# -- read             --
# ----------------------
function read_nodenum(skipnum)
    """
    xmax : 仮想セルも含めたnodeのxの数
    ymax : 仮想セルも含めたnodeのyの数
    """
    fff=[]
    open("grid/nodesnum", "r") do f
        fff=read(f,String)
    end
    fff=split(fff,"\n",keepempty=false)
    num_nodes=length(fff)-skipnum

    for i in 1+skipnum:length(fff)
        fff[i]=replace(fff[i]," \r" => "")
    end

    temp = split(fff[2]," ")
    xmax = parse(Int64,temp[1])
    ymax = parse(Int64,temp[2])

    return xmax, ymax
end

function read_nodes(skipnum,xmax,ymax)
    """
    nodes[i][j][k]
    i : x点の番号
    j : y点の番号
    k=1 : 点のx座標
    k=2 : 点のy座標
    """

    fff=[]
    open("grid/nodes", "r") do f
        fff=read(f,String)
    end
    fff=split(fff,"\n",keepempty=false)
    num_nodes=length(fff)-skipnum

    for i in 1+skipnum:length(fff)
        fff[i]=replace(fff[i]," \r" => "")
    end

    nodes=zeros(xmax,ymax,2)
    for i in 1:num_nodes
        temp=split(fff[i+skipnum]," ")

        x = parse(Int64,temp[1])
        y = parse(Int64,temp[2])
        nodes[x,y,1]=parse(Float64,temp[3])
        nodes[x,y,2]=parse(Float64,temp[4])
    end
    return nodes
end

function read_vecA(skipnum,xmax,ymax)
    """ vecAx[i][j][k]
        i : x番号
        j : y番号
        k=1 : x vec
        k=2 : y vec
    """
    fff=[]
    open("grid/vecAx", "r") do f
        fff=read(f,String)
    end
    fff=split(fff,"\n",keepempty=false)
    num_cell=length(fff)-skipnum

    for i in 1+skipnum:length(fff)
        fff[i]=replace(fff[i]," \r" => "")
    end

    vecAx=zeros(xmax,ymax-1,2)
    for i in 1:num_cell
        temp=split(fff[i+skipnum]," ")

        x = parse(Int64,temp[1])
        y = parse(Int64,temp[2])
        vecAx[x,y,1]=parse(Float64,temp[3])
        vecAx[x,y,2]=parse(Float64,temp[4])
    end

    fff=[]
    open("grid/vecAy", "r") do f
        fff=read(f,String)
    end
    fff=split(fff,"\n",keepempty=false)
    num_cell=length(fff)-skipnum

    for i in 1+skipnum:length(fff)
        fff[i]=replace(fff[i]," \r" => "")
    end

    vecAy=zeros(xmax-1,ymax,2)
    for i in 1:num_cell
        temp=split(fff[i+skipnum]," ")

        x = parse(Int64,temp[1])
        y = parse(Int64,temp[2])
        vecAy[x,y,1]=parse(Float64,temp[3])
        vecAy[x,y,2]=parse(Float64,temp[4])
    end
    return vecAx, vecAy
end

function read_allgrid()
    skip=1
    xmax,ymax=read_nodenum(skip)
    nodes=read_nodes(skip,xmax,ymax)
    vecAx, vecAy=read_vecA(skip,xmax,ymax)
    println("fin read grid")
    return xmax,ymax,nodes,vecAx,vecAy
end


text:
#
# This file is part of the DiscreteEvents.jl Julia package, MIT license
#
# Paul Bayer, 2020
#
# This is a Julia package for discrete event simulation
#

println("... basic tests: resources  ...")

R = Resource{Int}()
@test R.items isa Deque{Int}
@test R.capacity == typemax(Int)
R = Resource{Int}(5)
@test R.capacity == 5
@test capacity(R) == 5
@test length(R) == 0
@test isempty(R)
@test !isready(R)

@test !islocked(R)
lock(R)
@test islocked(R)
unlock(R)
@test !islocked(R)
@test trylock(R)
@test islocked(R)
t = @async trylock(R)
@test !fetch(t)
unlock(R)
@test !islocked(R)

for i in 1:5
    push!(R, i)
end
@test isready(R)
@test !isempty(R)
@test isfull(R)
@test length(R) == 5
@test_throws ArgumentError push!(R, 6)
@test_throws ArgumentError pushfirst!(R, 0)
@test first(R) == 1
@test last(R) == 5
@test pop!(R) == 5
@test popfirst!(R) == 1
@test first(R) == 2
pushfirst!(R, 1)
@test first(R) == 1
empty!(R)
@test isempty(R)

println("... basic tests: extended channel API  ...")
ch = Channel{Int}(5)
@test capacity(ch) == 5
@test isempty(ch)
for i in 1:5
    put!(ch, i)
end
@test first(ch) == 1
@test length(ch) == 5
@test isfull(ch)
@test take!(ch) == 1
@test length(ch) == 4
empty!(ch)
@test isempty(ch)


text:
const GMS_BLOCK_HEADER = ["Variables",
                      "Positive", # Variables
                      "Negative", # Variables
                      "Binary", # Variables
                      "Semicont", #Variables
                      "Integer",
                      "Equations"]

const PROBTYPE = ["LP",      # Linear Programming
                  "QCP",     # Quadratic Constraint Programming
                  "NLP",     # Nonlinear Programming
                  "DNLP",    # Nonlinear Programming with Discontinuous derivatives
                  "MIP",     # Mixed Integer Programming
                  "RMIP",    # Relaxed Mixed Integer Programming
                  "MIQCP",   # Mixed Integer Quadratic Constraint Programming
                  "MINLP",   # Mixed Integer Nonlinear Programming
                  "RMIQCP",  # Relaxed Mixed Integer Quadratic Constraint Programming
                  "RMINLP",  # Relaxed Mixed Integer Nonlinear Programming
                  "MCP",     # Mixed Complementarity Problem
                  "MPEC",	 # Mathematical Program iwth equilibrium Constraints
                  "CNS"		 # Constrained Nonlinear System
				  ]


text:
function makeArray(studentCount, studentDF, courseCount, coursesDF, weights)
	array = zeros((studentCount, courseCount)) .+ weights[1]
	for i = 1:courseCount
		for j = 1:studentCount
			if coursesDF.CourseCode[i] == studentDF.Choice1[j]
				array[j, i] = weights[2]
			elseif coursesDF.CourseCode[i] == studentDF.Choice2[j]
				array[j, i] = weights[3]
			elseif coursesDF.CourseCode[i] == studentDF.Choice3[j]
				array[j, i] = weights[4]
			elseif coursesDF.CourseCode[i] == studentDF.Choice4[j]
				array[j, i] = weights[5]
			elseif coursesDF.CourseCode[i] == studentDF.Choice5[j]
				array[j, i] = weights[6]
			end
		end
	end
	return array
end

mutable struct allotment
	studentCount::Int
	courseCount::Int
	times
	data
end

function squeezeAllotment(anAllotment::allotment)
	dataHolder = zeros(anAllotment.studentCount, anAllotment.courseCount)
	start = 1
	for i = 1:anAllotment.courseCount
		stop = start + anAllotment.times[i]
		dataHolder[:, i] = collect(Iterators.flatten(sum(anAllotment.data[:,start:stop-1], dims=2)))
		start = stop
	end
	anAllotment.data = dataHolder
	return anAllotment
end

function calcGoodness(anAllotment::allotment, choiceWeights, cpiArray)
	a = anAllotment.data .* choiceWeights
	choiceGoodness = sum(a, dims=2) ./ sum(a .!= 0, dims=2)
	choiceGoodness = replace!(choiceGoodness, NaN=>0)
	cpiRepeated = repeat(cpiArray, outer=[1, anAllotment.courseCount])
	SvC_cpi = anAllotment.data .* cpiRepeated
	cpiGoodness = sum(SvC_cpi, dims=1) ./ sum(SvC_cpi .!= 0, dims=1)
	cpiGoodness = replace!(cpiGoodness, NaN=>0)
	return choiceGoodness, cpiGoodness
end

function swapRows(anAllotment::allotment)
	swapID1 = rand(1:anAllotment.studentCount)
	swapID2 = rand(1:anAllotment.studentCount)
	temp = copy(anAllotment.data)
	tempVar = copy(temp[swapID1, :])
	temp[swapID1,:] = temp[swapID2,:]
	temp[swapID2,:] = tempVar
	return temp
end

# Running with no errors untill here

function calcChoiceCost(anAllotmentData, SvCArray)
	cost1 = sum(anAllotmentData .* SvCArray)
	return cost1
end

function calculateVariance(cpiArray, anAllotment, courseCount, times)
	cpiRepeated = repeat(cpiArray, outer=[1, courseCount])
	SvC_cpi = anAllotment .* cpiRepeated
	a = sum(SvC_cpi, dims=1) ./ sum(SvC_cpi .!= 0, dims=1)
	a = replace!(a, NaN=>0)
	variance = var(a)
	return variance
end

function allottedCourseGrade(studentCount, courseCount,studentDF, anAllotment, SvCArray)
	allotmentMatrix = anAllotment .* SvCArray
	temp = zeros((studentCount, 1))
	for i = 1:studentCount
		for j = 1:courseCount
			if allotmentMatrix[i,j] != 0
				if allotmentMatrix[i,j] == 10
					temp[i] = 0
					break
				else
					temp[i] = studentDF[i, Symbol("Grade"*string(Int(allotmentMatrix[i,j])))]
					# temp[i] = eval(Meta.parse("studentDF["*string(i)*",:].Grade"*string(Int(allotmentMatrix[i,j]))))
					break
				end
			end
		end
	end
	return sum(temp)
end

function calculateUtility(anAllotmentData, costWeights, choiceIdx, choiceWeights, courseCount, studentCount, studentsDF, cpiArray, times)
	choiceSum = calcChoiceCost(anAllotmentData, choiceWeights)
	cost1 = costWeights[1]*choiceSum
	cost2 = costWeights[2]*allottedCourseGrade(studentCount, courseCount, studentsDF, anAllotmentData, choiceWeights)
	var = calculateVariance(cpiArray, anAllotmentData, courseCount, times)
	cost3 = costWeights[3]/var
	totalUtility = cost1 + cost2 + cost3
	return totalUtility
end

function runMCMC(currentAllotment, nIters, studentCount, courseCount, costWeights, choiceIdx, choiceWeights, studentsDF, cpiArray, times)
	beta = 10*log10(1 + studentCount)
	utility = zeros((nIters, 1))
	for i = 1:nIters
		println(i)
		newAllotmentData = swapRows(currentAllotment)

		u1 = calculateUtility(newAllotmentData, costWeights, choiceIdx, choiceWeights, courseCount, studentCount, studentsDF, cpiArray, times)
		u2 = calculateUtility(currentAllotment.data, costWeights, choiceIdx, choiceWeights, courseCount, studentCount, studentsDF, cpiArray, times)
		utility[i] = u1
		println(u1)
		if u1 >= u2
			currentAllotment.data = newAllotmentData
		else
			beta = 10*log10(1+(i-1))
			probVar = exp(beta*(u1 - u2))
			randNum = rand()
			if randNum < probVar
				currentAllotment.data = newAllotmentData
			else
				currentAllotment.data = currentAllotment.data
			end
		end
	end
	return currentAllotment, utility
end

function writePerformance(anAllotment, choiceGoodnessOld, cpiGoodnessOld, choiceGoodnessNew, cpiGoodnessNew, utility)
	save("results/myfile.jld", "allotment", allotment, "choiceGoodnessOld", choiceGoodnessOld, "cpiGoodnessOld", cpiGoodnessOld, "choiceGoodnessNew", choiceGoodnessNew, "cpiGoodnessNew", cpiGoodnessNew, "utility", utility)
end

max_stars_repo_name: saurabhkm/juliaAllotter

max_issues_repo_name: saurabhkm/juliaAllotter

max_forks_repo_name: saurabhkm/juliaAllotter


text:
export InDims, OutDims, MovingWindow
using ..Cubes.Axes: get_descriptor, findAxis, Axes
using ...YAXArrays: YAXDefaults
using DataFrames: DataFrame
using YAXArrayBase: yaxcreate

"""
    MovingWindow(desc, pre, after)

Constructs a `MovingWindow` object to be passed to an `InDims` constructor to define that
the axis in `desc` shall participate in the inner function (i.e. shall be looped over), but inside
the inner function `pre` values before and `after` values after the center value
will be passed as well.

For example passing `MovingWindow("Time", 2, 0)` will loop over the time axis and
always pass the current time step plus the 2 previous steps. So in the inner function
the array will have an additional dimension of size 3.
"""
struct MovingWindow
    desc::Any
    pre::Int
    after::Int
end
Axes.get_descriptor(m::MovingWindow) = MovingWindow(get_descriptor(m.desc), m.pre, m.after)
Axes.findAxis(m::MovingWindow, c) = findAxis(m.desc, c)

wrapWorkArray(::Type{Array}, a, axes) = a
wrapWorkArray(T, a, axes) =
    yaxcreate(T, a, map(axsym, axes), map(i -> i.values, axes), Dict{String, Any}())

abstract type ProcFilter end
struct AllMissing <: ProcFilter end
struct NValid <: ProcFilter
    n::Int
end
struct AnyMissing <: ProcFilter end
struct AnyOcean <: ProcFilter end
struct NoFilter <: ProcFilter end
struct StdZero <: ProcFilter end
struct UserFilter{F} <: ProcFilter
    f::F
end

checkskip(::NoFilter, x) = false
checkskip(::AllMissing, x) = all(ismissing, x)
checkskip(::AnyMissing, x) = any(ismissing, x)
checkskip(nv::NValid, x) = count(!ismissing, x) <= nv.n
checkskip(uf::UserFilter, x) = uf.f(x)
checkskip(::StdZero, x) = all(i -> i == x[1], x)
docheck(pf::ProcFilter, x)::Bool = checkskip(pf, YAXArrayBase.getdata(x))
docheck(pf::Tuple, x) = reduce(|, map(i -> docheck(i, x), pf))

getprocfilter(f::Function) = (UserFilter(f),)
getprocfilter(pf::ProcFilter) = (pf,)
getprocfilter(pf::NTuple{N,<:ProcFilter}) where {N} = pf

"""
    InDims(axisdesc...;...)

Creates a description of an Input Data Cube for cube operations. Takes a single
  or multiple axis descriptions as first arguments. Alternatively a MovingWindow(@ref) struct can be passed to include
  neighbour slices of one or more axes in the computation.
  Axes can be specified by their
  name (String), through an Axis type, or by passing a concrete axis.

### Keyword arguments

* `artype` how shall the array be represented in the inner function. Defaults to `Array`, alternatives are `DataFrame` or `AsAxisArray`
* `filter` define some filter to skip the computation, e.g. when all values are missing. Defaults to
    `AllMissing()`, possible values are `AnyMissing()`, `AnyOcean()`, `StdZero()`, `NValid(n)`
    (for at least n non-missing elements). It is also possible to provide a custom one-argument function
    that takes the array and returns `true` if the compuation shall be skipped and `false` otherwise.
* `window_oob_value` if one of the input dimensions is a MowingWindow, this value will be used to fill out-of-bounds areas
"""
mutable struct InDims
    axisdesc::Tuple
    artype::Any
    procfilter::Tuple
    window_oob_value::Any
end
function InDims(
    axisdesc::Union{String,CubeAxis,Symbol,MovingWindow}...;
    artype = Array,
    filter = AllMissing(),
    window_oob_value = missing,
)
    descs = get_descriptor.(axisdesc)
    InDims(descs, artype, getprocfilter(filter), window_oob_value)
end

struct OutDims
    axisdesc::Any
    backend::Symbol
    backendargs::Any
    update::Bool
    artype::Any
    chunksize::Any
    outtype::Union{Int,DataType}
end
"""
    OutDims(axisdesc;...)

Creates a description of an Output Data Cube for cube operations. Takes a single
  or a Vector/Tuple of axes as first argument. Axes can be specified by their
  name (String), through an Axis type, or by passing a concrete axis.

- `axisdesc`: List of input axis names
- `backend` : specifies the dataset backend to write data to, must be either :auto or a key in `YAXArrayBase.backendlist`
- `update` : specifies wether the function operates inplace or if an output is returned
- `artype` : specifies the Array type inside the inner function that is mapped over
- `chunksize`: A Dict specifying the chunksizes for the output dimensions of the cube, or `:input` to copy chunksizes from input cube axes or `:max` to not chunk the inner dimensions
- `outtype`: force the output type to a specific type, defaults to `Any` which means that the element type of the first input cube is used
"""
function OutDims(
    axisdesc...;
    backend = :auto,
    update = false,
    artype = Array,
    chunksize = YAXDefaults.chunksize[],
    outtype = 1,
    backendargs...,
)
    descs = get_descriptor.(axisdesc)
    OutDims(descs, backend, backendargs, update, artype, chunksize, outtype)
end

registerDATFunction(a...; kwargs...) =
    @warn("Registration does not exist anymore, ignoring....")

max_stars_repo_name: meggart/YAXArrays

max_forks_repo_name: meggart/YAXArrays


text:
using JLD2; using PyPlot; using Statistics; using DelimitedFiles

# Load in the posteriors:
@load "../../../data/T1_hmc_total_02212020.jld2" state_total

clf()
# Make a confidence plot of eccentricity vectors, and
# show the longitudes of periastron.

nplanet = 7

e2 = 0.015; e1 = -0.015
ncos = 30
nsin = 30
evecgrid = zeros(nplanet,ncos,nsin)
ecosgrid = zeros(ncos)
esingrid = zeros(nsin)
for i=1:ncos
  ecosgrid[i] = (i-0.5)*(e2-e1)/ncos+e1
  esingrid[i] = (i-0.5)*(e2-e1)/nsin+e1
end

nsamp = 224000

conf = [0.95, 0.683]
nconf = length(conf)
level_ecc = zeros(nplanet,nconf)
planet_name=["b","c","d","e","f","g","h"]

cname = ["C0","C1","C2","C3","C4","C5","C6","C7"]

for j=1:nplanet
  ncosp = 10000
  # Set up vectors for samples:
  ecossamp= zeros(nsamp)
  esinsamp = zeros(nsamp)
  # Loop over number of samples:
  for isamp=1:nsamp
    # Draw from posterior distribution:
    # Compute the mass of planet in Earth mass units:
    ecossamp[isamp] = state_total[(j-1)*5+4,isamp]
    esinsamp[isamp] = state_total[(j-1)*5+5,isamp]
    # Now compute density distribution:
    if ecossamp[isamp] > e1 && ecossamp[isamp] < e2 && esinsamp[isamp] > e1 && esinsamp[isamp] < e2
      ir = ceil(Int64,(esinsamp[isamp]-e1)/(e2-e1)*ncos)
      im = ceil(Int64,(ecossamp[isamp]-e1)/(e2-e1)*nsin)
      evecgrid[j,ir,im] +=1.0
    end
  end
  # Set levels:
  evecgrid[j,:,:] /= float(sum(evecgrid[j,:,:]))
  evecgrid_sort = sort(vec(evecgrid[j,:,:]))
  i0 = nsin*ncos
  evecgrid_cum = evecgrid_sort[i0]
  while evecgrid_cum < conf[1] && i0 > 1
    i0 -= 1
    for k=1:nconf
      if evecgrid_cum < conf[k] && (evecgrid_cum+evecgrid_sort[i0]) > conf[k]
        level_ecc[j,k] = 0.5*(evecgrid_sort[i0]+evecgrid_sort[i0+1])
      end
    end
    evecgrid_cum += evecgrid_sort[i0]
  end
  println(j," ",mean(ecossamp)," ",std(ecossamp)," fractional error: ",std(ecossamp)/mean(ecossamp)*100.0,"%")
end
mopt = readdlm("../../../data/optimum_values.txt")
figure(figsize=(6,6))
# Plot maximum likelihood values:
for j=1:nplanet
  cs = contour(ecosgrid,esingrid,evecgrid[j,:,:],levels = level_ecc[j,:],colors=cname[j],linewidth=2.0)
  plot(mopt[(j-1)*5+4],mopt[(j-1)*5+5],"o",label=planet_name[j],color=cname[j])
  text(mopt[(j-1)*5+4]+0.0005,mopt[(j-1)*5+5]-0.0005,planet_name[j],color=cname[j],fontsize=15)
end
#legend(fontsize=15)
xlabel(L"$e \cos{\omega}$",fontsize=15)
ylabel(L"$e \sin{\omega}$",fontsize=15)
axis([-0.0125,0.0125,-0.0125,0.0125])
plot([-0.015,0.015],[0,0],linestyle=":",color="k")
plot([0,0],[-0.015,0.015],linestyle=":",color="k")
tight_layout()
savefig("../esin_vs_ecos.pdf",bbox_inches="tight")


text:
n = 8 # rand(1:10)
@test matrixdepot("frank", n) == matrixdepot("frank", Float64, n, 0)

A = zeros(Float64, n, n)
for i = 1:n, j = 1:n
    A[i,j] = i<=j ? n + 1 -j : j == i - 1 ? n - j : zero(Float64)
end
@test matrixdepot("frank", n) ≈ A

println("'frank' passed test...")


text:
mutable struct PWGrid
    Ns::Array{Int64,1}
    LatVecs::Array{Float64,2}
    RecVecs::Array{Float64,2}
    Npoints::Int
    Ω::Float64
    R::Array{Float64,2}
    G::Array{Float64,2}
    G2::Array{Float64}
end

function PWGrid( Ns::Array{Int,1},LatVecs::Array{Float64,2} )
    Npoints = prod(Ns)
    RecVecs = 2*pi*inv(LatVecs')
    Ω = det(LatVecs)
    R,G,G2 = init_grids( Ns, LatVecs, RecVecs )
    return PWGrid( Ns, LatVecs, RecVecs, Npoints, Ω, R, G, G2 )
end

function mm_to_nn(mm::Int,S::Int)
    if mm > S/2
        return mm - S
    else
        return mm
    end
end

function init_grids( Ns, LatVecs, RecVecs )
    #
    Npoints = prod(Ns)
    #
    R = Array(Float64,3,Npoints)
    ip = 0
    for k in 0:Ns[3]-1
    for j in 0:Ns[2]-1
    for i in 0:Ns[1]-1
        ip = ip + 1
        R[1,ip] = LatVecs[1,1]*i/Ns[1] + LatVecs[1,2]*j/Ns[2] + LatVecs[1,3]*k/Ns[3]
        R[2,ip] = LatVecs[2,1]*i/Ns[1] + LatVecs[2,2]*j/Ns[2] + LatVecs[2,3]*k/Ns[3]
        R[3,ip] = LatVecs[3,1]*i/Ns[1] + LatVecs[3,2]*j/Ns[2] + LatVecs[3,3]*k/Ns[3]
    end
    end
    end
    #
    G  = Array(Float64,3,Npoints)
    G2 = Array(Float64,Npoints)
    ip    = 0
    for k in 0:Ns[3]-1
    for j in 0:Ns[2]-1
    for i in 0:Ns[1]-1
        gi = mm_to_nn( i, Ns[1] )
        gj = mm_to_nn( j, Ns[2] )
        gk = mm_to_nn( k, Ns[3] )
        ip = ip + 1
        G[1,ip] = RecVecs[1,1]*gi + RecVecs[1,2]*gj + RecVecs[1,3]*gk
        G[2,ip] = RecVecs[2,1]*gi + RecVecs[2,2]*gj + RecVecs[2,3]*gk
        G[3,ip] = RecVecs[3,1]*gi + RecVecs[3,2]*gj + RecVecs[3,3]*gk
        G2[ip] = G[1,ip]^2 + G[2,ip]^2 + G[3,ip]^2
    end
    end
    end
    return R,G,G2
end


text:
using MPCCModelManager, Test, SafeTestsets
@time begin
@time @safetestset "Model Procesing" begin include("test_proc_model.jl") end
end

# using MPCCModelManager
# using Test

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


text:
# ------------------------------------------------------------------------------------------
# # `Table` Usage
#
# In this notebook we'll show common `Table` operations using the stock market data
# introduced in the previous notebook.  `NDSparse` operations are *nearly* identical, so we
# will focus on `Table`.  The functions we cover here are:
#
# 1. `select`
# 1. `filter`
# 1. `map`
# 1. `reduce`
# 1. `groupreduce`
# 1. `groupby`
# 1. `summarize`
# 1. `columns`/`rows`
# 1. `join`
# 1. `merge`
#
# Each of the above functions has detailed inline documentation, accessed from a Julia REPL
# with `?select`, for example.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# # Begin by Loading Data
#
# - Let's load the data we saved in the previous notebook:
# ------------------------------------------------------------------------------------------

using JuliaDB

# Print table rather than column summary
IndexedTables.set_show_compact!(false)

# loadtable("stocksample"; filenamecol = :Ticker, indexcols = [:Ticker, :Date]);
t = load("stocks.jdb")

# ------------------------------------------------------------------------------------------
# # Return a Subset of Columns:
#
# - We can use `select` to return a selector (introduced in the previous notebook) applied
# to a table.
#
# ## `select(table, selection)`
#
# - When multiple selectors are involved, rows are "passed around" as a `NamedTuple`.
# - A function paired with multiple selections must then accept a `NamedTuple`.
#
# - For example, to calculate the range of stock prices for each day we can:
#   1. Select `:High` and `:Low`
#   1. Pair it with the anonymous function `row -> row.High - row.Low`
# ------------------------------------------------------------------------------------------

select(t, (:High, :Low) => row -> row.High - row.Low)

# ------------------------------------------------------------------------------------------
# # Return a Subset of Rows:
#
# - We can get the rows that satisfy some condition (when a function returns true) with the
# syntax:
#
# ## `filter(function, table; selection)`
#
# - Here we retrieve the data for AMZN (Amazon) by getting the rows for which `Ticker ==
# "amzn.us.txt"`.
# ------------------------------------------------------------------------------------------

filter(x -> x == "amzn.us.txt", t; select = :Ticker)

# ------------------------------------------------------------------------------------------
# # Apply a Function to a Selection:
#
# - We can use `map` to apply a function on a selection of a table with the syntax below:
#
# ## `map(function, table; select)`
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# - If `select` is not provided, each full row will be passed to the function.
# - Here we return the first item in each row:
# ------------------------------------------------------------------------------------------

map(first, t)

# ------------------------------------------------------------------------------------------
# - Note that `map` and `select` can often be used to produce the same result since
# selections can be paired with a function.
# - For example, we previously used
#
#     ```julia
#     select(t, (:High, :Low) => row -> row.High - row.Low)
#     ```
#
#     to calculate stock price ranges.  Equivalently, we can use:
# ------------------------------------------------------------------------------------------

map(r -> r.High - r.Close, t)

# ------------------------------------------------------------------------------------------
# # `reduce`
#
# - `reduce` applies a function (`reducer`) pair-wise to a selection through the syntax:
#
# ## `reduce(reducer, table; select)`
#
# - For example, if a table is four rows long, `reduce(reducer, t)` is equivalent to
#
# ```julia
# out = reducer(row1, row2)
# out = reducer(out, row3)
# out = reducer(out, row4)
# ```
#
# - In order to be meaningful, the `reducer` must have the associative property:
#
# $$(A + B) + C = A + (B + C)$$
#
#
# ------------------------------------------------------------------------------------------

reduce(+, t; select = :Volume)

# ------------------------------------------------------------------------------------------
# You can also `reduce` with estimators from **OnlineStats** (more on this later):
# ------------------------------------------------------------------------------------------

using OnlineStats

reduce(Sum(Int), t; select = :Volume)

# ------------------------------------------------------------------------------------------
# # `groupreduce`
#
# - Like `reduce`, `groupreduce` applies a reducer pair-wise to table elements.
# - However, the reducer is applied separately across groups (unique values of another
# selection).
# - The syntax is:
#
# ## `groupreduce(reducer, table, by; selection)`
#
# - For example, we can find the total number of trades for each stock by calculating the
# sum of `:Volume`, grouped by `:Ticker`:
# ------------------------------------------------------------------------------------------

groupreduce(+, t, :Ticker; select = :Volume)

# ------------------------------------------------------------------------------------------
# # `groupby`
#
# - `groupby` applies a function to each group subset (not pair-wise like `reduce`) through
# the syntax:
#
# ## `groupby(function, table [, by]; select)`
#
# - Here we get the mean and standard deviation of closing price for each stock:
# ------------------------------------------------------------------------------------------

groupby((mean, std), t, :Ticker; select = :Close)

# ------------------------------------------------------------------------------------------
# # `summarize`
#
# - `summarize` applies a function (or functions) column-wise.  The syntax is:
#
# ## `summarize(function, table, by; select)`
# ------------------------------------------------------------------------------------------

summarize((mean, std), t, :Ticker; select = (:Open, :Close))

# ------------------------------------------------------------------------------------------
# # AoS and SoA
#
# - We can retrieve the table as a "struct of arrays" (`NamedTuple` of `Vector`s) or as an
# "array of structs" (`Vector` of `NamedTuple`s) via `columns` and `rows`, respectively.
#
# ## `columns(t; selection)`
#
# ## `rows(t; selection)`
# ------------------------------------------------------------------------------------------

# NamedTuple of Vectors
columns(t)[1]

# Vector of NamedTuples
rows(t)[1]

# ------------------------------------------------------------------------------------------
# # Joins
#
# ## `join(left, right; how, <options>)`
#
# Join tables together based on matching keys.
#
# - `how` can be one of `:inner`, `:left`, `:outer`, or`:anti`
# - `<options>`: `rkey`, `lkey` (default to indexed variable), `rselect`, `lselect`
# ------------------------------------------------------------------------------------------

t1 = table(@NT(x=1:5, y = rand(5)); pkey = :x)

t2 = table(@NT(x=3:7, z = rand(5)); pkey = :x)

# try :inner, :outer, :left
tjoin = join(t1, t2; how = :inner)

# ------------------------------------------------------------------------------------------
# # Merging
#
# - A `merge` results in a table that is still ordered by the primary key(s).
# ------------------------------------------------------------------------------------------

t3 = table(@NT(x=11:15, y = randn(5)), pkey = :x)

merge(t1, t3)

max_stars_repo_name: grenkoca/JuliaTutorials

max_issues_repo_name: grenkoca/JuliaTutorials

max_forks_repo_name: grenkoca/JuliaTutorials


text:
### A Pluto.jl notebook ###
# v0.19.0

using Markdown
using InteractiveUtils

# ╔═╡ 891c5658-b5e3-4cb8-8a21-9a1c28e122a5
md"""
为提高科学研究可重复性， 本站汇总了项目组近期开展项目的原始数据及数据分析过程。
"""

# ╔═╡ 95c53824-7376-4f0a-90c2-8863a4ddd0de
const NOTEBOOKS = [
    ("后悔情绪理解",  "Regret/Regret.jl"),
    ("脑电简单统计",  "Disjunctions/DisjunctionStatistics.jl"),
    ("脑电事件提取",  "Disjunctions/DisjunctionEvents.jl"),
    ("脑电事件转化",  "Disjunctions/RefactorTriggers.jl"),
    ("视觉观点采择",  "VisualPerspective/VisualPerspective.jl"),
    ("自闭症情绪理解", "AuditoryEmotion/AuditoryEmotion.jl"),
    ("否定加工过程",  "Negation/NegationWithoutLanguage.jl"),
    ("内隐心理理论",  "ImplicitTheoryOfMind/ImplicitTheoryOfMind.jl")
];

# ╔═╡ 82b7a596-1b6c-11ec-1731-5b8c4d67cc2c
let
function _linkname(path, nb, basedir)
    if haskey(ENV, "html_export") && ENV["html_export"] == "true"
        joinpath(basedir, "$(splitext(nb)[1]).html")
    else
        "open?path=" * joinpath(path, nb)
    end
end
function list_notebooks(file, basedir = "")
    path = joinpath(@__DIR__, "notebooks")
    sp = splitpath(file)
	filename = split(sp[end], "#")[1]

    list = join(["1. [$name]($(_linkname(path, nb, basedir)))" for (name, nb) in NOTEBOOKS], "\n")

    Markdown.parse("""## 项目目录 \n $list """)
end
list_notebooks(@__FILE__, "notebooks")
end

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
"""

# ╔═╡ 00000000-0000-0000-0000-000000000002
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised

julia_version = "1.7.2"
manifest_format = "2.0"

[deps]
"""

# ╔═╡ Cell order:
# ╟─891c5658-b5e3-4cb8-8a21-9a1c28e122a5
# ╟─82b7a596-1b6c-11ec-1731-5b8c4d67cc2c
# ╟─95c53824-7376-4f0a-90c2-8863a4ddd0de
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

using BAT
using Test

@testset "error_log" begin
    @test BAT.enable_error_log(false) === nothing
    @test BAT.error_log() === nothing
    @test_throws ArgumentError BAT.@throw_logged ArgumentError("foo")
    @test BAT.error_log() === nothing

    @test BAT.enable_error_log(true) === nothing
    @test BAT.error_log() isa AbstractVector{<:BAT.ErrLogEntry}
    @test_throws ArgumentError BAT.@throw_logged ArgumentError("foo")
    @test last(BAT.error_log()).error == ArgumentError("foo")

    @test_throws ErrorException try
        BAT.@throw_logged ArgumentError("bar")
    catch
        BAT.@rethrow_logged ErrorException("Some rethrow")
    end

    @test BAT.error_log()[end-1].error == ArgumentError("bar")
    @test BAT.error_log()[end].error == ErrorException("Some rethrow")

    BAT.enable_error_log(false)
    @test BAT.error_log() === nothing
end


text:
module TestInterface

using Phylo
using DataFrames
using Compat.Test

@testset "Build and tear down trees" begin
    @testset "For $TreeType" for TreeType in
        [NamedTree,
         BinaryTree{DataFrame, Vector{Float64}},
         BinaryTree{DataFrame, Vector{String}}]

        species = ["Dog", "Cat", "Human", "Potato", "Apple"]
        tree = TreeType(species)
        othernodes = ["Some 1", "Some 2"]
        @test othernodes == addnodes!(tree, othernodes)
        extra = addnodes!(tree, 1)
        @test isa(extra[1], String)
        append!(othernodes, extra)
        also = addnode!(tree)
        @test also ∉ othernodes
        @test_throws ErrorException addnode!(tree, also)
        @test_throws ErrorException addnodes!(tree, [also])
        @test_throws ErrorException deletenode!(tree, "really not")
        @test also == deletenode!(tree, also)
        also = addnode!(tree)
        push!(othernodes, also)
        @test Set(othernodes) ⊆ Set(nodenameiter(tree))
        innodes = append!(copy(species), othernodes)
        allnodes = reverse(innodes)
        pop!(innodes)
        branches =
            map(innodes) do node
                itr = Iterators.filter(name ->
                                       hasoutboundspace(tree, name) &&
                                       name != node, allnodes)
                addbranch!(tree, first(itr), node)
            end
        @test Set(branches) == Set(branchnameiter(tree))
        @test validate(tree)
        @test_throws ErrorException addbranch!(tree, allnodes[1], allnodes[2])
        who = getparent(tree, species[1])
        b = getinbound(tree, species[1])
        node1 = getnode(tree, species[1])
        @test !isunattached(node1) &&
            hasoutboundspace(node1) && !hasinboundspace(node1) &&
            outdegree(node1) == 0 && indegree(node1) == 1
        @test hasbranch(tree, b)
        @test dst(getbranch(tree, b)) == species[1]
        @test src(getbranch(tree, b)) == who
        @test b == deletebranch!(tree, b)
        @test !hasbranch(tree, b)
        @test isunattached(tree, species[1]) &&
            hasoutboundspace(tree, species[1]) &&
            hasinboundspace(tree, species[1]) &&
            outdegree(tree, species[1]) == 0 &&
            indegree(tree, species[1]) == 0
        @test_throws ErrorException getbranch(tree, b)
        branches = collect(Iterators.filter(name -> name != b, branches))
        @test Set(branches) == Set(branchnameiter(tree))
        b3 = getinbound(tree, species[2])
        source = src(tree, b3)
        destination = dst(tree, b3)
        @test b3 == changesrc!(tree, b3, who)
        @test b3 == changedst!(tree, b3, species[1])
        b2 = addbranch!(tree, source, destination)
        push!(branches, b2)
        @test Set(branches) == Set(branchnameiter(tree))
        @test species[1] == deletenode!(tree, species[1])
        @test !hasnode(tree, species[1])
        @test species[1] == branch!(tree, who, destination=species[1])
        @test hasnode(tree, species[1])
        @test validate(tree)
        @test all(map(==, Pair(first(branchiter(tree))),
                      Tuple(tree, first(branchnameiter(tree)))))
        @test all(map(==, Pair(tree, first(branchnameiter(tree))),
                      Tuple(first(branchiter(tree)))))
        @test all(map(node -> isleaf(tree, node), species))
        @test all(map(node -> !isroot(tree, node) &&
                      !isunattached(tree, node) &&
                      !isinternal(tree, node), species))
        @test Set(map(pair -> src(pair[2]), collect(getbranches(tree)))) ∪
            Set(map(pair -> dst(pair[2]), collect(getbranches(tree)))) ==
            Set(allnodes)
    end
end

end


text:
"""
    calculatesignal_gbloch_ide(α, TRF, TR, ω0, B1, m0s, R1f, R2f, Rx, R1s, T2s[; grad_list, Ncyc=2, output=:complexsignal])

Calculate the signal or magnetixation evolution with the full generalized Bloch model assuming a super-Lorentzian lineshape (slow).

The simulation assumes a sequence of rectangluar RF-pulses with varying flip angles α and RF-pulse durations TRF, but a fixed repetition time TR. Further, it assumes balanced gradient moments.

# Arguments
- `α::Vector{<:Number}`: Array of flip angles in radians
- `TRF::Vector{<:Number}`: Array of the RF-pulse durations in seconds
- `TR::Number`: Repetition time in seconds
- `ω0::Number`: Off-resonance frequency in rad/s
- `B1::Number`: Normalized transmit B1 field, i.e. B1 = 1 corresponds to a well-calibrated B1 field
- `m0s::Number`: Fractional size of the semi-solid pool; should be in range of 0 to 1
- `R1f::Number`: Longitudinal relaxation rate of the free pool in 1/seconds
- `R2f::Number`: Transversal relaxation rate of the free pool in 1/seconds
- `Rx::Number`: Exchange rate between the two spin pools in 1/seconds
- `R1s::Number`: Longitudinal relaxation rate of the semi-solid pool in 1/seconds
- `T2s::Number`: Transversal relaxationt time of the semi-solid pool in seconds

Optional:
- `grad_list=[]`: Vector to indicate which gradients should be calculated; the vector can either be empty `[]` for no gradient, or contain any subset/order of `grad_list=[grad_m0s(), grad_R1s(), grad_R2f(), grad_Rx(), grad_R1s(), grad_T2s(), grad_ω0(), grad_B1()]`; the derivative wrt. to apparent `R1a = R1f = R1s` can be calculated with `grad_R1a()`
- `Ncyc=2`: The magnetization is initialized with thermal equilibrium and then performed Ncyc times and only the last cycle is stored. The default value is usually a good approximation for antiperiodic boundary conditions. Increase the number for higher precision at the cost of computation time.
- `output=:complexsignal`: The defaul keywords triggers the function to output a complex-valued signal (`xf + 1im yf`); the keyword `output=:realmagnetization` triggers an output of the entire (real valued) vector `[xf, yf, zf, xs, zs]`
- `greens=(greens_superlorentzian, dG_o_dT2s_x_T2s_superlorentzian)`: Tuple of a Greens function `G(κ) = G((t-τ)/T2s)` and its partial derivative wrt. T2s, multiplied by T2s `∂G((t-τ)/T2s)/∂T2s * T2s`. This package supplies the three Greens functions `greens=(greens_superlorentzian, dG_o_dT2s_x_T2s_superlorentzian)` (default), `greens=(greens_lorentzian, dG_o_dT2s_x_T2s_lorentzian)`, and `greens=(greens_gaussian, dG_o_dT2s_x_T2s_gaussian)`

# Examples
```jldoctest
julia> calculatesignal_gbloch_ide(ones(100)*π/2, ones(100)*5e-4, 4e-3, 0, 1, 0.2, 0.3, 15, 20, 2, 10e-6)
100×1 Matrix{ComplexF64}:
   -0.00796631644501679 + 0.0im
  0.0012590590419314775 - 0.0im
  -0.006088855588122639 + 0.0im
  0.0024187389404174958 - 0.0im
  -0.004361339395232399 + 0.0im
   0.003489135821004988 - 0.0im
 -0.0027633710614651317 + 0.0im
   0.004483217941394422 - 0.0im
 -0.0012812573517352534 + 0.0im
    0.00540885403467719 - 0.0im
                        ⋮
   0.017760808273049045 - 0.0im
   0.017576118974646872 + 0.0im
   0.017813950945910605 - 0.0im
    0.01764385633550697 + 0.0im
   0.017863575855931853 - 0.0im
    0.01770692603385291 + 0.0im
   0.017909914934082043 - 0.0im
    0.01776565037329648 + 0.0im
   0.017953184893717555 - 0.0im

julia> calculatesignal_gbloch_ide(ones(100)*π/2, ones(100)*5e-4, 4e-3, 0, 1, 0.2, 0.3, 15, 20, 2, 10e-6; grad_list=[grad_R1f(), grad_T2s()], output=:realmagnetization)
100×15 transpose(::Matrix{Float64}) with eltype Float64:
 -0.00796627   0.0   0.000637773  …   0.0   -10.8757  -335.26   0.0
  0.00125903  -0.0  -0.00700671      -0.0   125.882   -326.977  0.0
 -0.00608882   0.0   0.00185086       0.0   -30.4187  -317.56   0.0
  0.00241873  -0.0  -0.00520622      -0.0    96.1776  -309.906  0.0
 -0.00436133   0.0   0.00296471       0.0   -47.5803  -302.948  0.0
  0.003489    -0.0  -0.00354518   …  -0.0    69.5148  -298.697  0.0
 -0.00276366   0.0   0.00399588       0.0   -62.8453  -294.886  0.0
  0.00448273  -0.0  -0.00200673      -0.0    45.3179  -292.783  0.0
 -0.00128187   0.0   0.00495478       0.0   -76.6573  -290.321  0.0
  0.00540814  -0.0  -0.000578836     -0.0    23.1756  -289.245  0.0
  ⋮                               ⋱
  0.0177563   -0.0   0.0175372       -0.0  -290.779   -349.855  0.0
  0.0175716    0.0   0.0177845        0.0  -295.347   -350.002  0.0
  0.0178094   -0.0   0.0176073       -0.0  -292.44    -350.163  0.0
  0.0176393    0.0   0.0178359        0.0  -296.668   -350.3    0.0
  0.017859    -0.0   0.0176727    …  -0.0  -294.001   -350.451  0.0
  0.0177024    0.0   0.0178838        0.0  -297.914   -350.579  0.0
  0.0179053   -0.0   0.0177335       -0.0  -295.467   -350.72   0.0
  0.0177611    0.0   0.0179286        0.0  -299.09    -350.84   0.0
  0.0179486   -0.0   0.0177902       -0.0  -296.845   -350.972  0.0
```
"""
function calculatesignal_gbloch_ide(α, TRF, TR, ω0, B1, m0s, R1f, R2f, Rx, R1s, T2s; grad_list=[], Ncyc=2, output=:complexsignal, greens=(greens_superlorentzian, dG_o_dT2s_x_T2s_superlorentzian))
    ω1 = α ./ TRF

    # Define G(τ) and its derivative as ApproxFuns
    G = interpolate_greens_function(greens[1], 0, maximum(TRF) / T2s)
    if any(isa.(grad_list, grad_T2s))
        dG_o_dT2s_x_T2s = interpolate_greens_function(greens[2], 0, maximum(TRF) / T2s)
    else
        dG_o_dT2s_x_T2s = []
    end

    # initialization and memory allocation
    N_s = 5 * (1 + length(grad_list))
    alg = MethodOfSteps(DP8())
    s = zeros(N_s, length(TRF))
    m0 = zeros(N_s)
    m0[3] = (1 - m0s)
    m0[4] = m0s
    m0[5] = 1
    mfun(p, t; idxs=nothing) = typeof(idxs) <: Number ? m0[idxs] : m0

    # prep pulse
    sol = solve(DDEProblem(apply_hamiltonian_gbloch!, m0, mfun, (0, TRF[2]), (-ω1[2] / 2, B1, ω0, m0s, R1f, R2f, Rx, R1s, T2s, G, dG_o_dT2s_x_T2s, grad_list)), alg)
    m0 = sol[end]

    T_FP = (TR - TRF[2]) / 2 - TRF[1] / 2
    sol = solve(ODEProblem(apply_hamiltonian_freeprecession!, m0, (0, T_FP), (ω0, m0s, R1f, R2f, Rx, R1s, grad_list)), Tsit5())
    m0 = sol[end]

    for ic = 0:(Ncyc - 1)
        # free precession for TRF/2
        sol = solve(ODEProblem(apply_hamiltonian_freeprecession!, m0, (0, TRF[1] / 2), (ω0, m0s, R1f, R2f, Rx, R1s, grad_list)), Tsit5())
        m0 = sol[end]

        # inversion pulse with crusher gradients (assumed to be instantanious)
        u00 = m0[1:3]
        m0[1:5:end] .*= -sin(B1 * ω1[1] * TRF[1] / 2)^2
        m0[2:5:end] .*= sin(B1 * ω1[1] * TRF[1] / 2)^2
        m0[3:5:end] .*= cos(B1 * ω1[1] * TRF[1])

        # calculate saturation of RF pulse
        sol = solve(DDEProblem(apply_hamiltonian_gbloch_inversion!, m0, mfun, (0, TRF[1]), ((-1)^(1 + ic) * ω1[1], B1, ω0, m0s, 0, 0, 0, 0, T2s, G, dG_o_dT2s_x_T2s, grad_list)), alg)
        m0[4:5:end] = sol[end][4:5:end]

        for i in eachindex(grad_list)
            if isa(grad_list[i], grad_B1)
                m0[5i + 1] -= u00[1] * sin(B1 * ω1[1] * TRF[1] / 2) * cos(B1 * ω1[1] * TRF[1] / 2) * ω1[1] * TRF[1]
                m0[5i + 2] += u00[2] * sin(B1 * ω1[1] * TRF[1] / 2) * cos(B1 * ω1[1] * TRF[1] / 2) * ω1[1] * TRF[1]
                m0[5i + 3] -= u00[3] * sin(B1 * ω1[1] * TRF[1]) * ω1[1] * TRF[1]
            end
        end

        # free precession
        T_FP = TR - TRF[2] / 2
        TE = TR / 2
        sol = solve(ODEProblem(apply_hamiltonian_freeprecession!, m0, (0.0, T_FP), (ω0, m0s, R1f, R2f, Rx, R1s, grad_list)), Tsit5(), save_everystep=false, saveat=TE)
        s[:,1] = sol[2]
        s[1:5:end,1] .*= (-1)^(1 + ic)
        s[2:5:end,1] .*= (-1)^(1 + ic)
        m0 = sol[end]

        for ip = 2:length(TRF)
            sol = solve(DDEProblem(apply_hamiltonian_gbloch!, m0, mfun, (0.0, TRF[ip]), ((-1)^(ip + ic) * ω1[ip], B1, ω0, m0s, R1f, R2f, Rx, R1s, T2s, G, dG_o_dT2s_x_T2s, grad_list)), alg)
            m0 = sol[end]

            T_FP = TR - TRF[ip] / 2 - TRF[mod(ip, length(TRF)) + 1] / 2
            TE = TR / 2 - TRF[ip] / 2
            sol = solve(ODEProblem(apply_hamiltonian_freeprecession!, m0, (0.0, T_FP), (ω0, m0s, R1f, R2f, Rx, R1s, grad_list)), Tsit5(), save_everystep=false, saveat=TE)
            if sol.t[2] / TE - 1 > 1e-10
                throw(DimensionMismatch("sol.t[2] is not equal to TE"))
            end
            s[:,ip] = sol[2]
            s[1:5:end,ip] .*= (-1)^(ip + ic)
            s[2:5:end,ip] .*= (-1)^(ip + ic)
            m0 = sol[end]
        end
    end
    s = transpose(s)
    if output == :complexsignal
        s = s[:, 1:5:end] + 1im * s[:, 2:5:end]
    end
    return s
end

"""
    calculatesignal_graham_ode(α, TRF, TR, ω0, B1, m0s, R1f, R2f, Rx, R1s, T2s[; grad_list, Ncyc=2, output=:complexsignal])

Calculate the signal or magnetixation evolution with Graham's spectral model assuming a super-Lorentzian lineshape.

The simulation assumes a sequence of rectangluar RF-pulses with varying flip angles α and RF-pulse durations TRF, but a fixed repetition time TR. Further, it assumes balanced gradient moments.

# Arguments
- `α::Vector{<:Number}`: Array of flip angles in radians
- `TRF::Vector{<:Number}`: Array of the RF-pulse durations in seconds
- `TR::Number`: Repetition time in seconds
- `ω0::Number`: Off-resonance frequency in rad/s
- `B1::Number`: Normalized transmit B1 field, i.e. B1 = 1 corresponds to a well-calibrated B1 field
- `m0s::Number`: Fractional size of the semi-solid pool; should be in range of 0 to 1
- `R1f::Number`: Longitudinal relaxation rate of the free pool in 1/seconds
- `R2f::Number`: Transversal relaxation rate of the free pool in 1/seconds
- `Rx::Number`: Exchange rate between the two spin pools in 1/seconds
- `R1s::Number`: Longitudinal relaxation rate of the semi-solid pool in 1/seconds
- `T2s::Number`: Transversal relaxationt time of the semi-solid pool in seconds

Optional:
- `grad_list=[]`: Vector to indicate which gradients should be calculated; the vector can either be empty `[]` for no gradient, or contain any subset/order of `grad_list=[grad_m0s(), grad_R1f(), grad_R2f(), grad_Rx(), grad_R1s(), grad_T2s(), grad_ω0(), grad_B1(), grad_R1a()]`
- `Ncyc=2`: The magnetization is initialized with thermal equilibrium and then performed Ncyc times and only the last cycle is stored. The default value is usually a good approximation for antiperiodic boundary conditions. Increase the number for higher precision at the cost of computation time.
- `output=:complexsignal`: The defaul keywords triggers the function to output a complex-valued signal (`xf + 1im yf`); the keyword `output=:realmagnetization` triggers an output of the entire (real valued) vector `[xf, yf, zf, xs, zs]`

# Examples
```jldoctest
julia> calculatesignal_graham_ode(ones(100)*π/2, ones(100)*5e-4, 4e-3, 0, 1, 0.2, 0.3, 15, 20, 2, 10e-6)
100×1 Matrix{ComplexF64}:
  -0.00807345119202598 + 0.0im
 0.0012686432904785246 - 0.0im
 -0.006178694437440372 + 0.0im
  0.002435865817868477 - 0.0im
 -0.004437476277388273 + 0.0im
  0.003516464650880384 - 0.0im
 -0.002831554256573095 + 0.0im
  0.004523902251349811 - 0.0im
 -0.001342229952563084 + 0.0im
  0.005454562035727423 - 0.0im
                       ⋮
  0.018148222055411015 - 0.0im
  0.017957696642062964 + 0.0im
  0.018204860894110032 - 0.0im
  0.018029363128597056 + 0.0im
  0.018257820158882408 - 0.0im
   0.01809616895266688 + 0.0im
  0.018307337749218777 - 0.0im
   0.01815844455461564 + 0.0im
  0.018353636251669653 - 0.0im

julia> calculatesignal_graham_ode(ones(100)*π/2, ones(100)*5e-4, 4e-3, 0, 1, 0.2, 0.3, 15, 20, 2, 10e-6; grad_list=[grad_R1f(), grad_T2s()], output=:realmagnetization)
100×15 transpose(::Matrix{Float64}) with eltype Float64:
 -0.0080756    0.0   0.000643162  …   0.0   -10.4986  -323.634  0.0
  0.00126867  -0.0  -0.00710692      -0.0   123.078   -316.358  0.0
 -0.00618067   0.0   0.00186482       0.0   -29.4458  -307.862  0.0
  0.00243634  -0.0  -0.0052899       -0.0    94.1692  -300.821  0.0
 -0.00443746   0.0   0.00298646       0.0   -46.1422  -294.116  0.0
  0.00351386  -0.0  -0.00361718   …  -0.0    68.262   -289.421  0.0
 -0.00282882   0.0   0.00402793       0.0   -61.0236  -285.242  0.0
  0.00451808  -0.0  -0.00206741      -0.0    44.7888  -282.867  0.0
 -0.00133608   0.0   0.00499316       0.0   -74.4075  -280.79   0.0
  0.00544896  -0.0  -0.00062663      -0.0    23.2536  -280.148  0.0
  ⋮                               ⋱
  0.0181755   -0.0   0.0179431       -0.0  -285.071   -337.48   0.0
  0.0179911    0.0   0.0181963        0.0  -289.591   -337.634  0.0
  0.0182343   -0.0   0.018019        -0.0  -286.777   -337.804  0.0
  0.0180645    0.0   0.0182531        0.0  -290.964   -337.948  0.0
  0.0182893   -0.0   0.0180898    …  -0.0  -288.383   -338.106  0.0
  0.0181329    0.0   0.0183062        0.0  -292.261   -338.241  0.0
  0.0183407   -0.0   0.0181559       -0.0  -289.895   -338.39   0.0
  0.0181968    0.0   0.0183559        0.0  -293.488   -338.516  0.0
  0.0183888   -0.0   0.0182175       -0.0  -291.318   -338.656  0.0
```
"""
function calculatesignal_graham_ode(α, TRF, TR, ω0, B1, m0s, R1f, R2f, Rx, R1s, T2s; grad_list=[], Ncyc=2, output=:complexsignal)
    ω1 = α ./ TRF

    # initialization and memory allocation
    N_s = 5 * (1 + length(grad_list))
    s = zeros(N_s, length(TRF))
    m0 = zeros(N_s)
    m0[3] = (1 - m0s)
    m0[4] = m0s
    m0[5] = 1.0

    # prep pulse
    sol = solve(ODEProblem(apply_hamiltonian_graham_superlorentzian!, m0, (0.0, TRF[2]), (-ω1[2] / 2, B1, ω0, TRF[2], m0s, R1f, R2f, Rx, R1s, T2s, grad_list)), Tsit5(), save_everystep=false)
    m0 = sol[end]

    T_FP = TR / 2 - TRF[2] / 2 - TRF[1] / 2
    sol = solve(ODEProblem(apply_hamiltonian_freeprecession!, m0, (0.0, T_FP), (ω0, m0s, R1f, R2f, Rx, R1s, grad_list)), Tsit5(), save_everystep=false)
    m0 = sol[end]

    for ic = 0:(Ncyc - 1)
        # free precession for TRF/2
        sol = solve(ODEProblem(apply_hamiltonian_freeprecession!, m0, (0.0, TRF[1] / 2), (ω0, m0s, R1f, R2f, Rx, R1s, grad_list)), Tsit5(), save_everystep=false)
        m0 = sol[end]

        # inversion pulse with crusher gradients (assumed to be instantanious)
        u00 = m0[1:3]
        m0[1:5:end] .*= -sin(B1 * ω1[1] * TRF[1] / 2)^2
        m0[2:5:end] .*= sin(B1 * ω1[1] * TRF[1] / 2)^2
        m0[3:5:end] .*= cos(B1 * ω1[1] * TRF[1])

        # calculate saturation of RF pulse
        sol = solve(ODEProblem(apply_hamiltonian_graham_superlorentzian_inversionpulse!, m0, (0, TRF[1]), ((-1)^(1 + ic) * ω1[1], B1, ω0, TRF[1], m0s, 0, 0, 0, 0, T2s, grad_list)), Tsit5(), save_everystep=false)
        m0[4:5:end] = sol[end][4:5:end]

        for i in eachindex(grad_list)
            if isa(grad_list[i], grad_B1)
                m0[5i + 1] -= u00[1] * sin(B1 * ω1[1] * TRF[1] / 2) * cos(B1 * ω1[1] * TRF[1] / 2) * ω1[1] * TRF[1]
                m0[5i + 2] += u00[2] * sin(B1 * ω1[1] * TRF[1] / 2) * cos(B1 * ω1[1] * TRF[1] / 2) * ω1[1] * TRF[1]
                m0[5i + 3] -= u00[3] * sin(B1 * ω1[1] * TRF[1]) * ω1[1] * TRF[1]
            end
        end

        # free precession
        T_FP = TR - TRF[2] / 2
        TE = TR / 2
        sol = solve(ODEProblem(apply_hamiltonian_freeprecession!, m0, (0.0, T_FP), (ω0, m0s, R1f, R2f, Rx, R1s, grad_list)), Tsit5(), save_everystep=false, saveat=TE)
        s[:,1] = sol[2]
        s[1:5:end,1] .*= (-1)^(1 + ic)
        s[2:5:end,1] .*= (-1)^(1 + ic)
        m0 = sol[end]

        for ip = 2:length(TRF)
            sol = solve(ODEProblem(apply_hamiltonian_graham_superlorentzian!, m0, (0.0, TRF[ip]), ((-1)^(ip + ic) * ω1[ip], B1, ω0, TRF[ip], m0s, R1f, R2f, Rx, R1s, T2s, grad_list)), Tsit5(), save_everystep=false)
            m0 = sol[end]

            T_FP = TR - TRF[ip] / 2 - TRF[mod(ip, length(TRF)) + 1] / 2
            TE = TR / 2 - TRF[ip] / 2
            sol = solve(ODEProblem(apply_hamiltonian_freeprecession!, m0, (0.0, T_FP), (ω0, m0s, R1f, R2f, Rx, R1s, grad_list)), Tsit5(), save_everystep=false, saveat=TE)
            if sol.t[2] / TE - 1 > 1e-10
                throw(DimensionMismatch("sol.t[2] is not equal to TE"))
            end
            s[:,ip] = sol[2]
            s[1:5:end,ip] .*= (-1)^(ip + ic)
            s[2:5:end,ip] .*= (-1)^(ip + ic)
            m0 = sol[end]
        end
    end
    s = transpose(s)
    if output == :complexsignal
        s = s[:, 1:5:end] + 1im * s[:, 2:5:end]
    end
    return s
end


text:
#=
# 2D advection using DG
=#

### If the Finch package has already been added, use this line #########
using Finch # Note: to add the package, first do: ]add "https://github.com/paralab/Finch.git"

### If not, use these four lines (working from the examples directory) ###
# if !@isdefined(Finch)
#     include("../Finch.jl");
#     using .Finch
# end
##########################################################################

init_finch("advection2d");

useLog("advection2dlog")

domain(2)
solverType(DG)
functionSpace(order=3)
timeStepper(CRANK_NICHOLSON, cfl=0.1)

mesh(QUADMESH, elsperdim=8)

u = variable("u")
testSymbol("v")

T = 2;
timeInterval(T)
initial(u, "0")

### Add boundary regions and set BCs ###############
addBoundaryID(2, (x,y) -> (x <= 0 || y <= 0))
addBoundaryID(3, (x,y) -> (x < 0.2 && y < 0.2))
boundary(u, 1, NO_BC)
boundary(u, 2, DIRICHLET, 0)
boundary(u, 3, DIRICHLET, "sin(pi*(t - 0.7*(x+y)))^2")
####################################################

coefficient("a", type=VECTOR, [0.5, 0.5])
coefficient("beta", 0)
weakForm(u, "Dt(u*v) + a*grad(u)*v + surface(-dot(a, normal())*u*v + dot(a, normal())*ave(u)*v - 0.5*(1-beta)*dot(a, jump(u))*v)")

solve(u);

# Uncomment to plot
# using Plots
# pyplot();
# display(plot(Finch.grid_data.allnodes[1,:], Finch.grid_data.allnodes[2,:], u.values[:], st=:surface, reuse=false))

finalize_finch()

max_stars_repo_name: paralab/Finch

max_issues_repo_name: paralab/Finch

max_forks_repo_name: paralab/Finch


text:
module VarianceTest

using Test
using DiffPrecTest
using ProximalBase
using Random, LinearAlgebra, SparseArrays, Statistics, Distributions

function varianceSymAlt(
    Sx,
    Sy,
    X,
    Y,
    ω::Vector,
    indω::Vector{Int64},
    Δ::Vector,
    indΔ::Vector{Int64}
    )

    nx, p = size(X)
    ny    = size(Y, 1)

    q = zeros(nx)
    r = zeros(ny)

    # compute statistics
    H = (kron(Sx, Sy) + kron(Sy, Sx)) / 2.
    t = dot(ω, H[indω, indΔ] * Δ) - dot(ω, Sy[indω] - Sx[indω])

    # compute qk
    for k=1:nx
        H = (kron(X[k, :]*X[k, :]', Sy) + kron(Sy, X[k, :]*X[k, :]')) / 2.
        q[k] = dot(ω, H[indω, indΔ] * Δ) - dot(ω, Sy[indω] - (X[k, :]*X[k, :]')[indω])
    end

    # compute rk
    for k=1:ny
        H = (kron(Y[k, :]*Y[k, :]', Sx) + kron(Sx, Y[k, :]*Y[k, :]')) / 2.
        r[k] = dot(ω, H[indω, indΔ] * Δ) - dot(ω, (Y[k, :]*Y[k, :]')[indω] - Sx[indω])
    end

    σ1 = sum(abs2, q) / (nx - 1) - nx / (nx - 1) * t^2
    σ2 = sum(abs2, r) / (ny - 1) - ny / (ny - 1) * t^2

    return σ1/nx + σ2/ny
end

@testset "sym variance" begin
    Random.seed!(134)
    n = 100

    Ωx = [1. 0.4; 0.4 1]
    Ωy = [1. 0.3; 0.3 1]
    Σx = inv(Symmetric(Ωx))
    Σy = inv(Symmetric(Ωy))
    dist_X = MvNormal(convert(Matrix, Σx))
    dist_Y = MvNormal(convert(Matrix, Σy))

    X = rand(dist_X, n)'
    Y = rand(dist_Y, n)'

    Sx = Symmetric( X'X / n )
    Sy = Symmetric( Y'Y / n )

    ind = [1, 2, 3, 4]
    H = (kron(Sx, Sy) + kron(Sy, Sx)) / 2.
    tmp = zeros(length(ind))
    tmp[1] = 1.
    ω = H \ tmp
    Δ = H \ (Sy[ind] - Sx[ind])

    @test variance(SymmetricOracleNormal(), Sx, Sy, X, Y, ω, ind, Δ, ind) ≈ varianceSymAlt(Sx, Sy, X, Y, ω, ind, Δ, ind)
end

end


text:
baremodule Folds

export DistributedEx, SequentialEx, ThreadedEx

using Transducers: DistributedEx, SequentialEx, ThreadedEx

module Core
abstract type ReducerFunctionAndFoldable <: Function end
end

import DefineSingletons

function reduce end
function mapreduce end
function foreach end
function collect end
function copy end
function map end
function map! end

DefineSingletons.@def_singleton all isa Core.ReducerFunctionAndFoldable
DefineSingletons.@def_singleton any isa Core.ReducerFunctionAndFoldable
DefineSingletons.@def_singleton count isa Core.ReducerFunctionAndFoldable
DefineSingletons.@def_singleton extrema isa Core.ReducerFunctionAndFoldable
DefineSingletons.@def_singleton findall isa Core.ReducerFunctionAndFoldable
DefineSingletons.@def_singleton findfirst isa Core.ReducerFunctionAndFoldable
DefineSingletons.@def_singleton findlast isa Core.ReducerFunctionAndFoldable
DefineSingletons.@def_singleton maximum isa Core.ReducerFunctionAndFoldable
DefineSingletons.@def_singleton minimum isa Core.ReducerFunctionAndFoldable
DefineSingletons.@def_singleton prod isa Core.ReducerFunctionAndFoldable
DefineSingletons.@def_singleton sum isa Core.ReducerFunctionAndFoldable
DefineSingletons.@def_singleton unique isa Core.ReducerFunctionAndFoldable

function issorted end

function set end
function dict end

module Implementations

using Accessors: @set
using BangBang:
    Empty,
    SingletonVector,
    append!!,
    collector,
    finish!,
    push!!,
    union!!,
    merge!!,
    setdiff!!
using Base: HasShape, IteratorSize, add_sum, mapreduce_empty, mul_prod, reduce_empty
using Baselet: Baselet
using DefineSingletons: @def_singleton
using InitialValues: @def_monoid, SpecificInitialValue, asmonoid
using MicroCollections: emptyshim, singletonshim
using Referenceables: isreferenceable, referenceable
using Requires: @require
using Transducers:
    Consecutive,
    DefaultInit,
    DistributedEx,
    Executor,
    Filter,
    IdentityTransducer,
    InitOf,
    Map,
    MapSplat,
    OnInit,
    PreferParallel,
    ReduceIf,
    SequentialEx,
    ThreadedEx,
    Transducers,
    dcollect,
    executor_for,
    extract_transducer,
    maybe_set_simd,
    reducingfunction,
    right,
    tcollect,
    transduce,
    unreduced,
    wheninit

using ..Folds.Core: ReducerFunctionAndFoldable
using ..Folds: Folds

include("utils.jl")
include("reduce.jl")
include("collect.jl")
include("misc.jl")

function __init__()
    @require OnlineStatsBase="925886fa-5bf2-5e8e-b522-a9147a512338" begin
        include("interop/onlinestats.jl")
    end
end

end

module Testing
include("testing.jl")
end

Implementations.define_docstrings()

end


text:
# Script that executes the Kabsch algorithm to align two paired sets of points
# by minimising the RMSD between them.
# Takes 2 .xyz files as input
# Program executed as
# 'julia kabsch_optimiser.jl geomP.xyz geomQ.xyz'
# Darya Shchepanovska
# 15th July 2019

using DelimitedFiles

using LinearAlgebra
using Statistics

function xyz2matrix(input_xyz)
    # Imports .xyz file format as data frame, and converts it to an N x 3 Matrix
    raw_xyz=readdlm(input_xyz)
    natoms=size(raw_xyz)[1]
    just_coords=Array{Float64}(raw_xyz[2:natoms,2:4])
    return just_coords
end

function translate_to_centroid(coord_matrix)
    # Normalises the molecular coordinates by centering them.
    center = [mean(coord_matrix[:,1]);mean(coord_matrix[:,2]);mean(coord_matrix[:,3])]
    centroid = transpose(center)
    translated_geom = broadcast(-,coord_matrix,centroid)
    return translated_geom
end

function cross_covariance_matrix(Pmatrix,Qmatrix)
    # Cross covariance matrix gives measure of variability between two matrices
    CCmatrix = transpose(Pmatrix) * Qmatrix
    return CCmatrix
end

function optimal_rotation_matrix(CCmatrix)
    # Returns 3x3 matrix that can be applied to P to get Q
    ORmatrix = sqrt(transpose(CCmatrix)*CCmatrix)*inv(CCmatrix)
    return ORmatrix
end

Pgeom = xyz2matrix(ARGS[1])
Qgeom = xyz2matrix(ARGS[2])

normalisedP = (translate_to_centroid(Pgeom))
normalisedQ = (translate_to_centroid(Qgeom))

xcov = cross_covariance_matrix(normalisedP,normalisedQ)
orot = optimal_rotation_matrix(xcov)

num_atoms = size(normalisedP)[1]
rotated = zeros(Float64,num_atoms,3)

for i=1:num_atoms
    rotated[i,:] = orot*normalisedP[i,:]
end

# display(rotated)

RMSD_value = norm(rotated-normalisedQ)

println("RMSD between these two structures is ",RMSD_value)

max_stars_repo_name: dasha-shchep/kabsch

max_issues_repo_name: dasha-shchep/kabsch

max_forks_repo_name: dasha-shchep/kabsch


text:
function validate(obj::ParameterValidator, arg0::RealVector)
    return jcall(obj, "validate", RealVector, (RealVector,), arg0)
end


text:
function solve_ip(input::Input)

    init_log_ip(input)
    n = input.n
    m = input.m
    rho = input.rho
    alpha = input.alpha
    tol = input.tol
    A = input.A
    b = input.b
    c = input.c

    nx = (1:n)
    np = (n+1:m+n)
    ns = (m+n+1:m+2n)

    x = ones(n)
    s = ones(n)
    p = ones(m)
    X = x
    S = s
    P = p
    gap = []
    mu = rho*x'*s/n

    for iter in 1:input.max_iter
        F = [A zeros(m,m) zeros(m,n);
            zeros(n,n) A' I(n);
            Diagonal(s) zeros(n,m) Diagonal(x)]

        g = [A*x - b; A'p + s + c; x.*s .- mu]

        d = - F \ g

        d_x = d[nx]
        d_p = d[np]
        d_s = d[ns]

        filter_x = d_x .< 0.0
        filter_p = d_p .< 0.0
        filter_s = d_s .< 0.0

        if sum(filter_x) == 0
            beta_x = 1.0
        else
            beta_x = min(1.0, alpha*minimum((-x./d_x)[filter_x]))
        end

        if sum(filter_p) == 0
            beta_p = 1.0
        else
            beta_p = min(1.0, alpha*minimum((-p./d_p)[filter_p]))
        end

        if sum(filter_s) == 0
            beta_s = 1.0
        else
            beta_s = min(1.0, alpha*minimum((-s./d_s)[filter_s]))
        end

        x += beta_x * d_x
        p += beta_s * d_p
        s += beta_s * d_s
        z = c'x
        w = b'p
        dual_inf = x'*s

        if z > 1/tol
            output = OutputIP(x,s,p,X,S,P,gap,mu,Inf,2,input.max_iter)
            last_log(input, output)
            return output
        end

        if w > 1/tol
            output = OutputIP(x,s,p,X,S,P,gap,mu,-Inf,3,input.max_iter)
            last_log(input, output)
            return output
        end
        X = [X x]
        S = [S s]
        P = [P p]
        push!(gap,dual_inf)
        iteration_log(input,iter,z,dual_inf)
        mu = rho*dual_inf/n
        if mu < tol
            output = OutputIP(x,s,p,X,S,P,gap,mu,z,1,iter)
            last_log(input, output)
            return output
        end
    end
    z = c'x
    output = OutputIP(x,s,p,X,S,P,gap,mu,z,0,input.max_iter)
    last_log(input, output)
    return output
end

function crossover(input::Input, output::OutputIP)
    m = input.m
    n = input.n
    crossover_tol = input.crossover_tol
    x = output.x
    midterm = MidTerm()

    base = []
    nbase = []

    for i in 1:n
        if x[i] > 1E-3
            push!(base,i)
        else
            push!(nbase,i)
        end
    end

    midterm.base = base
    midterm.nbase = nbase

    B = view(input.A,:,midterm.base)
    for (k,w) in enumerate(midterm.nbase)
        sol = B \ input.A[:,w]
        if norm(B*sol - input.A[:,w]) > crossover_tol
            push!(midterm.base, w)
            deleteat!(midterm.nbase, k)
            if length(midterm.base) == m
                break
            end
        end
    end

    init_log_simplex2(input)
    while midterm.termination_status == 0 && midterm.iter < input.max_iter
        iterate!(input, midterm)
    end
    update_midterm!(input, midterm)
    output = write_output(input, midterm)
    return output
end