text:
# now uses ParticleCollection from ParticleFilters.jl
# uncomment below if ParticleFilters ever dies or anything
#=
"""
Belief represented by an unweighted collection of particles
"""
type ParticleCollection{S}
    particles::Vector{S}
    ParticleCollection(particles) = new(particles)
    ParticleCollection() = new(S[])
end
ParticleCollection{S}(particles::AbstractVector{S}) = ParticleCollection{S}(particles)

function rand(rng::AbstractRNG, b::ParticleCollection, sample=nothing)
    return b.particles[rand(rng, 1:length(b.particles))]
end
mean(pc::ParticleCollection) = mean(pc.particles)
=#

"""
    uses_states_from_planner(belief)

Return true if the belief should use states simulated by the planner for belief estimation.

If this returns true, push!(belief, state) will be called whenever a new state for the belief node is generated by the planner. This is intended to allow users to be able to implement their own, more advanced particle filters that take advantage of the simulations run by the planner.
"""
uses_states_from_planner(::Any) = false
uses_states_from_planner(::ParticleCollection) = true

Base.push!{S}(b::ParticleCollection{S}, state::S) = push!(b.particles, state)

"""
Abstract base for a domain specific device to reinvigorate the particle collection when it has become depleted.

For use with POMCP, a subtype of this should implement the functions reinvigorate! and handle_unseen_observation
"""
abstract ParticleReinvigorator{S} <: POMDPs.Updater{ParticleCollection{S}}

"""
    reinvigorate!(pc::ParticleCollection, r::ParticleReinvigorator, old_node::BeliefNode, a, o)

Add states to pc to prevent particle depletion.

These states shoulc be consistent with the observation-action history accessible through old_node and a and o.
"""
function reinvigorate!(pc::ParticleCollection, r::ParticleReinvigorator, old_node::BeliefNode, a, o)
    error("""
          POMCP.jl reinvigorate! not implemented for reinvigorator $(typeof(r))\n
          argument types:
              b_old::$(typeof(b_old)),
              a::$(typeof(a)),
              o::$(typeof(o))
          Did you remember to explicitly import reinvigorate!?
          """)
end

"""
    handle_unseen_observation(r::ParticleReinvigorator, old_node::BeliefNode, a, o)

Create and return a new particle (state) collection of particles consistent with o.

This is called when o has not been previously seen during the state and action
"""
function handle_unseen_observation(r::ParticleReinvigorator, old_node::BeliefNode, a, o)
    error("""
          POMCP.jl: handle_unseen_observation not implemented for reinvigorator $(typeof(r))
          argument types:
              b_old::$(typeof(b_old)),
              a::$(typeof(a)),
              o::$(typeof(o))
          Did you remember to explicitly import handle_unseen_observation?
          """)
end

"""
Default reinvigorator - cannot do anything since there is no domain knowledge
"""
type DeadReinvigorator{S} <: ParticleReinvigorator{S} end

"""
Placeholder for a DeadReinvigorator so that the state type can be inferred during solve()
"""
type DefaultReinvigoratorStub end

function reinvigorate!(pc::ParticleCollection, r::DeadReinvigorator, old_node::BeliefNode, a, o)
    if length(pc.particles) == 0
        particle_depletion_error()
    end
    return pc # do nothing
end

handle_unseen_observation(r::DeadReinvigorator, old_node::BeliefNode, a, o) = particle_depletion_error()

particle_depletion_error() =
    error("""
          POMCP.jl: Particle Depletion! To fix this, you have three options:
                1) use more tree_queries (will only work for very small problems)
                2) implement a ParticleReinvigorator with reinvigorate!() and handle_unseen_observation()
                3) implement a more advanced updater for the agent (POMCP can use any
                   belief/state distribution that supports rand())
          """)

function update(::DeadReinvigorator, ::Any, ::Any, ::Any, b=nothing)
    error("update() is not implemented for DeadReinvigorator. update() should never be needed for a ParticleReinvigorator because it uses POMCP simulations to approximate belief updates.")
end


text:
include("../src/aliases.jl")
include("../src/entropy.jl")

P1 = Any[1,2,3,4]
P2 = Any[5,6,6,7]
P3 = Any[3,4,5,5]

m11 = mutual_information([P1,P1])
m12 = mutual_information([P1,P2])
m13 = mutual_information([P1,P3])
m23 = mutual_information([P2,P3])
m123 = mutual_information([P1,P2,P3])

println("m11: ",m11,"  m12: ",m12,"  m11: ",m11,"  m23: ",m23,"  m123: ",m123)

Pa = Any[1,1,2]   # For Dewar paper, p1 = 2/3
Pb = Any[1,2,2]   # For Dewar paper, p2 = 1/3
mPaPb =  mutual_information([P1,P2])
println("mPaPb: ",mPaPb)

function build_tbl( pvec::Vector{Float64} )
  m = length(pvec)
  result = zeros(Float64,(m,2))
  for i = 1:m
    result[i,:] = [pvec[i], 1-pvec[i]]
  end
  result/m
end

function build_tbl( parray::Array{Float64,2} )
  (m,n) = size(parray)
  result = zeros(Float64,(m,n+1))
  result[:,1:n] = parray
  diff_vec = [ 1.0-sum(parray[i,:]) for i = 1:m]
  @assert reduce( min, diff_vec ) >= 0.0
  result[:,n+1] = diff_vec
  result/m
end

# Generate a random m by n table whose rows sum to 1/m
function rand_tbl( m::Int64, n::Int64 )
  rtbl = rand(m,n)
  for i = 1:m
    rtbl[i,:] = rtbl[i,:]/sum(rtbl[i,:])/m
  end
  rtbl
end

pvec = [2/3,1/4,2/5,1/2]
ttbl = build_tbl(pvec)
std_mi = mutual_information(ttbl)
sh_mi = sherwin_mutual_information(ttbl)
@assert isapprox( std_mi, sh_mi )

parr = [[0.5 0.2]; [0.1 0.4]; [1/4 1/2]; [1/8 5/8]]
ttbl = build_tbl(parr)
std_mi = mutual_information(ttbl)
sh_mi = sherwin_mutual_information(ttbl)
@assert isapprox( std_mi, sh_mi )

m = 3
n = 5
rtbl = rand_tbl(m,n)
std_mi = mutual_information(rtbl)
sh_mi = sherwin_mutual_information(rtbl)
@assert isapprox( std_mi, sh_mi )

ext: jl

max_stars_repo_name: um-tech-evolution/nearly-neutral

max_issues_repo_name: um-tech-evolution/nearly-neutral

max_forks_repo_name: um-tech-evolution/nearly-neutral


text:
using MATLAB
using Base.Test

# test MMAT file I/O
fn = "$(tempname()).mat"

a32 = Int32[1 2 3; 4 5 6]
a64 = Int64[1 2 3; 4 5 6]
b = [1.2, 3.4, 5.6, 7.8]
c = [[0., 1.], [1., 2.], [1., 2., 3.]]
d = Dict("name"=>"MATLAB", "score"=>100.)

struct S
    x::Float64
    y::Bool
    z::Vector{Float64}
end

ss = S[S(1.0, true, [1., 2.]), S(2.0, false, [3., 4.])]

write_matfile(fn; a32=a32, a64=a64, b=b, c=c, d=d, ss=mxstructarray(ss))

r = read_matfile(fn)
@test isa(r, Dict{String, MxArray})
@test length(r) == 6

ra32 = jmatrix(r["a32"])
ra64 = jmatrix(r["a64"])
rb = jvector(r["b"])
rc = jvalue(r["c"])
rd = jdict(r["d"])
rss = r["ss"]

gc()  # make sure that ra, rb, rc, rd remain valid

@test ra32 == a32
@test ra64 == a64
@test rb == b
@test rc == c

@test rd["name"] == d["name"]
@test rd["score"] == d["score"]

@test is_struct(rss)
@test jscalar(get_field(rss, 1, "x")) == 1.0
@test jscalar(get_field(rss, 1, "y"))
@test jvector(get_field(rss, 1, "z")) == ss[1].z
@test jscalar(get_field(rss, 2, "x")) == 2.0
@test !jscalar(get_field(rss, 2, "y"))
@test jvector(get_field(rss, 2, "z")) == ss[2].z

# segfault on deleted references
s = MatFile(fn)
close(s)
@test_throws UndefRefError close(s)

rm(fn)


text:
module PrintUtils
  export
  print_iter_head,
  print_iter

  function print_iter_head()
    print_double_line()
    @printf("%10s  %20s\n", "Iter", "Cost")
    print_double_line()
  end

  function print_iter(iter :: Int, cost :: Real)
    @printf("%10d  %20.4e\n", iter, cost)
  end

  # some helper functions
  function print_double_line()
    @printf("==================================\n")
  end

  function print_blank_line()
    @printf("\n")
  end

end

max_stars_repo_name: upul/Senuk

max_issues_repo_name: upul/Senuk

max_forks_repo_name: upul/Senuk


text:
using JFVM
using PyPlot
include("CRT.jl")

#Calculating error as a function of timesteps
dt_mult = 0.1:0.1:2

# define the domain and constants
Lx = 1.0 # [m]
Nx = 100  # number of grids
ϕ  = 0.4 # porosity
u  = 1.0/(3600*24) # m/s, eq. 1 m/day
D  = [1e-9, 1e-9, 1e-9] # m^2/s diffusivity
K_eq = 10.0
k_f = 100.0
k_b = k_f/K_eq
nu_A = 1.0 # stoichiometric coefficient of component A
nu_B = 1.0 # stoichiometric coefficient of component B
nu_C = 1.0 # stoichiometric coefficient of component C
a = 1.0 #order of reaction
b = 1.0 #order of reaction
c = 1.0 #order of reaction
N_comp = 3 # number of components

# create the mesh
m = createMesh1D(Nx, Lx)  # mesh
# m = createMeshRadial2D(Nx, Nx, Lx, Lx)

# boundary conditions
C_bc = [1.0, 0.5, 0.0] # boundary concentration
C_0 = zeros(N_comp) # initial concentration

dt_ref = 0.1*minimum(dt_mult)*Lx/(u/ϕ)/Nx
T_ref = 1.0/(0.1*minimum(dt_mult))*Nx/2.0*dt_ref

err_array = zeros(5, length(dt_mult))

#Calculating reference "true solution"
tol = 1e-10 # convergence tolerance
C_ref = CRT.CRT_GIA_eps(m, Nx, ϕ, u, D, K_eq, nu_A, nu_B, nu_C, C_bc, C_0, N_comp, dt_ref, T_ref, tol)

for i=1:length(dt_mult)
  # Defining time-stepping parameters
  dt = dt_mult[i] * Lx/(u/ϕ)/Nx
  T = (1.0/dt_mult[i]) * Nx/2.0*dt

  #Calculating alternatives
  C_SNIA        = CRT.CRT_SNIA(m, Nx, ϕ, u, D, K_eq, nu_A, nu_B, nu_C, C_bc, C_0, N_comp, dt, T)
  C_SIA         = CRT.CRT_SIA(m, Nx, ϕ, u, D, K_eq, nu_A, nu_B, nu_C, C_bc, C_0, N_comp, dt, T, tol)
  C_GIA_eps     = CRT.CRT_GIA_eps(m, Nx, ϕ, u, D, K_eq, nu_A, nu_B, nu_C, C_bc, C_0, N_comp, dt, T, tol)
  C_GIA_kin_eq  = CRT.CRT_GIA_kin_eq(m, Nx, ϕ, u, D, K_eq, nu_A, nu_B, nu_C, a, b, c, C_bc, C_0, N_comp, dt, T, tol)
  C_GIA_kin_fb  = CRT.CRT_GIA_kin_fb(m, Nx, ϕ, u, D, k_f, k_b, nu_A, nu_B, nu_C, a, b, c, C_bc, C_0, N_comp, dt, T, tol)

  for j = 1:N_comp
    err_array[1, i] += sum(abs(internalCells(C_ref[j])-internalCells(C_SNIA[j])))/sum(abs(internalCells(C_ref[j])))
    err_array[2, i] += sum(abs(internalCells(C_ref[j])-internalCells(C_SIA[j])))/sum(abs(internalCells(C_ref[j])))
    err_array[3, i] += sum(abs(internalCells(C_ref[j])-internalCells(C_GIA_eps[j])))/sum(abs(internalCells(C_ref[j])))
    err_array[4, i] += sum(abs(internalCells(C_ref[j])-internalCells(C_GIA_kin_eq[j])))/sum(abs(internalCells(C_ref[j])))
    err_array[5, i] += sum(abs(internalCells(C_ref[j])-internalCells(C_GIA_kin_fb[j])))/sum(abs(internalCells(C_ref[j])))
  end

  println(i)
end

figure()
plot(dt_mult*Lx/(u/ϕ)/Nx, err_array[1,:])
plot(dt_mult*Lx/(u/ϕ)/Nx, err_array[2,:])
plot(dt_mult*Lx/(u/ϕ)/Nx, err_array[3,:])
plot(dt_mult*Lx/(u/ϕ)/Nx, err_array[4,:])
plot(dt_mult*Lx/(u/ϕ)/Nx, err_array[5,:])
plot([Lx/(u/ϕ)/Nx, Lx/(u/ϕ)/Nx],[0, maximum(err_array)], color=:black, linestyle="--")
ylim([0, maximum(err_array)])
grid("on")
legend(["SNIA", "SIA", "GIA ("L"$\epsilon$"")", "GIA ("L"$K_{eq}$"")", "GIA ("L"$k_fk_b$"")"])
xlabel(L"$\Delta t\ [s]$")
ylabel(L"||$C_{ref}-C||_1/||C_{ref}||_1$")
savefig("error_time_plot.png")

#CRT_SIA title:
#L"$\phi\frac{\partial C_{\alpha}}{\partial t}+\nabla\cdot(uC_{\alpha})-\nabla\cdot(D_{\alpha}\nabla C_{\alpha})+\nu_{\alpha}\epsilon=0,\quad K_{eq}=$""$K_eq"" (SIA)"

#CRT_SNIA title:
#L"$\phi\frac{\partial C_{\alpha}}{\partial t}+\nabla\cdot(uC_{\alpha})-\nabla\cdot(D_{\alpha}\nabla C_{\alpha})+\nu_{\alpha}\epsilon=0,\quad K_{eq}=$""$K_eq"" (SNIA)"

#CRT_GIA_eps title:
#L"$\phi\frac{\partial C_{\alpha}}{\partial t}+\nabla\cdot(uC_{\alpha})-\nabla\cdot(D_{\alpha}\nabla C_{\alpha})+\nu_{\alpha}\epsilon=0,\quad K_{eq}=$""$K_eq"" (FIM)"

#CRT_GIA_kin_eq title:
#L"$\phi\frac{\partial C_{\alpha}}{\partial t}+\nabla\cdot(uC_{\alpha})-\nabla\cdot(D_{\alpha}\nabla C_{\alpha})+\nu_{\alpha}(K_{eq}C_A^aC_B^b-C_C^c)=0,\quad K_{eq}=$""$K_eq"" (FIM)"

#CRT_GIA_kin_fb title:
#L"$\phi\frac{\partial C_{\alpha}}{\partial t}+\nabla\cdot(uC_{\alpha})-\nabla\cdot(D_{\alpha}\nabla C_{\alpha})+\nu_{\alpha}(k_{f}C_A^aC_B^b-k_{b}C_C^c)=0,\quad K_{eq}=$""$K_eq"" (FIM)"

max_stars_repo_name: FrederikLehn/coupled-solvers

max_issues_repo_name: FrederikLehn/coupled-solvers

max_forks_repo_name: FrederikLehn/coupled-solvers


text:
# @author:      0.382
# @description: 原子核相关的
#

# -------------------
# 同位素部分
# -------------------

"""
    struct Isotope
        Z::Int
        N::Int
    end
同位素。构造函数：
```julia
Isotope(iso::Isotope)
Isotope(Z::Integer, N::Integer)
Isotope(name::AbstractString)
```
"""
struct Isotope
    Z::Int
    N::Int
    Isotope(Z::Integer, N::Integer) = new(Z, N)
end

Isotope(iso::Isotope) = iso
Isotope(name::AbstractString) = parse(Isotope, name)

"""
    parse(::Type{Isotope}, str::AbstractString)
解析核素符号
"""
Base.parse(::Type{Isotope}, str::AbstractString)::Isotope = begin
    m = match(r"([A-Z][a-z]*)(\d{1,3})", str)
    m === nothing && error("match error in parsing '$str'")
    name = m.captures[1]
    _A = parse(Int, m.captures[2])
    ele = Element(name)
    _Z = getZ(ele)
    _N = _A - _Z
    @assert _N >= 0 "$str isotope has a negative neutron number $_N"
    Isotope(_Z, _N)
end

"获取同位素的质子数"
getZ(iso::Isotope) = iso.Z
"获取同位素的中子数"
getN(iso::Isotope) = iso.N
"获取同位素的核子数"
getA(iso::Isotope) = iso.Z + iso.N

"以文本形式显示同位素，比如`Ne20`"
Base.show(io::IO, ::MIME"text/plain", iso::Isotope) = begin
    ele = Element(iso.Z)
    is_none(ele) && error("no such element")
    print(io, symbol(ele), getA(iso))
end

"以latex形式显示同位素"
Base.show(io::IO, ::MIME"text/markdown", iso::Isotope) = begin
    ele = Element(iso.Z)
    is_none(ele) && error("no such element")
    show(io, "text/markdown", Markdown.parse("``^{$(getA(iso))}\\mathrm{$(symbol(ele))}``"))
end

# ------------------
# 单粒子轨道
# ------------------

const l_name_map = Dict(
     0 => "s",  1 => "p",  2 => "d",  3 => "f",
     4 => "g",  5 => "h",  6 => "i",  7 => "j",
     8 => "k",  9 => "l", 10 => "m", 11 => "n",
    12 => "o"
)

"""
    abstract type SingleParticleOrbit end
单粒子轨道
"""
abstract type SingleParticleOrbit end

"""
    struct JOrbit <: SingleParticleOrbit
        n::Int  # 主量子数
        l::Int  # 轨道角动量
        j::Int  # 总角动量
        tz::Int # 同位旋
    end
j-shceme 下的单粒子轨道。该结构体以及下面的`MOrbit`中，
总角动量和同位旋相关量都是使用物理上量子数的两倍来表示，以回避半整数
"""
struct JOrbit <: SingleParticleOrbit
    n::Int  # 主量子数
    l::Int  # 轨道角动量
    j::Int  # 总角动量
    tz::Int # 同位旋
    JOrbit(n, l, j, tz) = new(n, l, j, tz)
end

"""
    struct MOrbit <: SingleParticleOrbit
        n::Int  # 主量子数
        l::Int  # 轨道角动量
        j::Int  # 总角动量
        m::Int  # 总角动量的分量
        tz::Int # 同位旋
    end
m-scheme 下的单粒子轨道。
"""
struct MOrbit <: SingleParticleOrbit
    n::Int  # 主量子数
    l::Int  # 轨道角动量
    j::Int  # 总角动量
    m::Int  # 总角动量的分量
    tz::Int # 同位旋
    MOrbit(n, l, j, m, tz) = new(n, l, j, m, tz)
end

"""
    MOrbit(orbit::JOrbit, m::Integer)
通过 j-scheme 的单粒子轨道和``m``量子数来构造 m-scheme 轨道。
"""
MOrbit(orbit::JOrbit, m::Integer) = MOrbit(orbit.n, orbit.l, orbit.j, m, orbit.tz)

"""
    name(orbit::SingleParticleOrbit)
单粒子轨道的名字。例如`p(0d3/2)`，其中最外面的`p, n`表示质子或者中子，
`l > 12`以后的轨道角动量都用`"x"`表示。（我想也没人对这之后的轨道角动量感兴趣。）
"""
function name(orbit::SingleParticleOrbit)
    nucleon_type = orbit.tz == -1 ? "p" : "n"
    l_name = get(l_name_map, orbit.l, "x")
    "$(nucleon_type)($(orbit.n)$(l_name)$(orbit.j)/2)"
end

"显示j-scheme轨道"
Base.show(io::IO, orbit::JOrbit) = print(io, name(orbit))

"""
    struct NuclearShell
        orbits::Vector{JOrbit}
    end
原子核中壳的概念。实际上就是一系列的轨道组合而已，不保证真的具有物理意义。
按照壳模型计算的一般习惯，只需要j-scheme的轨道就够了。
"""
struct NuclearShell
    orbits::Vector{JOrbit}
    NuclearShell(orbits) = new(orbits)
end
NuclearShell(ns::NuclearShell) = ns

"""
    merge(ns1::NuclearShell, ns2::NuclearShell)
    merge(ns1::NuclearShell, ns2::NuclearShell, xns...)
合并两个壳。反正本代码中，壳只是一个计算上的概念，表示一系列轨道的集合。
把这些集合合并起来就是了。
"""
function Base.merge(ns1::NuclearShell, ns2::NuclearShell)
    NuclearShell(vcat(ns1.orbits, ns2.orbits))
end

function Base.merge(ns1::NuclearShell, ns2::NuclearShell, xns...)
    merge(NuclearShell(vcat(ns1.orbits, ns2.orbits)), xns...)
end

Base.merge(ns::NuclearShell) = ns

const s_shell = NuclearShell([
    JOrbit(0, 0, 1, -1), # 0s1/2
    JOrbit(0, 0, 1, 1)
])

const p_shell = NuclearShell([
    JOrbit(0, 1, 3, -1),  # 0p3/2
    JOrbit(0, 1, 3, 1),
    JOrbit(0, 1, 1, -1),  # 0p1/2
    JOrbit(0, 1, 1, 1)
])

const sd_shell = NuclearShell([
    JOrbit(0, 2, 5, -1),  # 0d5/2
    JOrbit(0, 2, 5, 1),
    JOrbit(1, 0, 1, -1),  # 1s1/2
    JOrbit(1, 0, 1, 1),
    JOrbit(0, 2, 3, -1),  # 0d3/2
    JOrbit(0, 2, 3, 1)
])

const pf_shell = NuclearShell([
    JOrbit(0, 3, 7, -1),  # 0f7/2
    JOrbit(0, 3, 7, 1),
    JOrbit(1, 1, 3, -1),  # 1p3/2
    JOrbit(1, 1, 3, 1),
    JOrbit(0, 3, 5, -1),  # 0f5/2
    JOrbit(0, 3, 5, 1),
    JOrbit(1, 1, 1, -1),  # 1p1/2
    JOrbit(1, 1, 1, 1)
])

"""
    j_orbits(ns::NuclearShell)
获得壳内所有的j-scheme轨道
"""
function j_orbits(ns::NuclearShell)
    ns.orbits
end

"""
    m_orbits(ns::NuclesrShell)
获得壳内所有的m-scheme轨道
"""
function m_orbits(ns::NuclearShell)
    orbits = Vector{MOrbit}()
    for orb in ns.orbits
        for m = -orb.j:2:orb.j
            push!(orbits, MOrbit(orb, m))
        end
    end
    return orbits
end

"""
    jsize(ns::NuclearShell)
壳内有几条j-scheme轨道。
"""
function jsize(ns::NuclearShell)
    length(ns.orbits)
end

"""
    msize(ns::NuclearShell)
计算壳内有几条m-scheme轨道，也就是能够容纳几个核子。
"""
function msize(ns::NuclearShell)
    sum = 0
    for orb in ns.orbits
        sum += orb.j + 1
    end
    return sum
end

"""
    psize(ns::NuclearShell)
计算壳内能够容纳几个质子。
"""
function psize(ns::NuclearShell)
    sum = 0
    for orb in ns.orbits
        (orb.tz == -1) && (sum += orb.j + 1)
    end
    return sum
end

"""
    nsize(ns::NuclearShell)
计算壳内能够容纳几个中子。
"""
function nsize(ns::NuclearShell)
    sum = 0
    for orb in ns.orbits
        (orb.tz == 1) && (sum += orb.j + 1)
    end
    return sum
end

"""
    max_pj(ns::NuclearShell)
计算壳内质子最大角动量。
"""
function max_pj(ns::NuclearShell)
    maxj = 1
    for orb in ns.orbits
        if orb.tz == -1
            maxj = max(maxj, orb.j)
        end
    end
    return maxj
end

"""
    max_nj(ns::NuclearShell)
计算壳内中子最大角动量
"""
function max_nj(ns::NuclearShell)
    maxj = 1
    for orb in ns.orbits
        if orb.tz == 1
            maxj = max(maxj, orb.j)
        end
    end
    return maxj
end

"""
    partition(n::Int, capacities::AbstractArray)
将 `n` 个小球放入 `length(capacities)` 个容器内，每个容器的容量由 `capacities` 数组指定
生成所有可能的放置方式
"""
function partition(n::Int, capacities::AbstractArray)
    result = Vector{Vector{Int}}()
    cap_size = length(capacities)
    next = (n::Int, pos::Int, line::Vector{Int}) -> begin
        if pos == cap_size
            if n <= capacities[end]
                xline = copy(line)
                xline[end] = n
                push!(result, xline)
            end
        else
            for i = 0:min(n, capacities[pos])
                xline = copy(line)
                xline[pos] = i
                next(n - i, pos + 1, xline)
            end
        end
        nothing
    end
    line = zeros(Int, length(capacities))
    next(n, 1, line)
    return result
end

"""
    binominal_safe(n::Integer, m::Integer)
判断这个二项式系数是否会溢出。
"""
function binominal_safe(n::Integer, m::Integer)
    n, m = promote(n, m)
    m = min(m, div(n, 2))
    try
        binomial(n, m)
    catch err
        if isa(err, OverflowError)
            return false
        else
            throw(err)
        end
    end
    return true
end

# 对于 julia 来说，`~m = -m - 1` 的小trick就不好用了
# 不过应该也无伤大雅
@inline m2index(m::Int) = m >= 0 ? m : (-m + 1)
@inline index2m(idx::Int) = isodd(idx) ? idx : (-idx + 1)

"""
    m_config_size(ns::NuclearShell, Z::Integer, N::Integer)
计算某一壳内，`Z+N`核子数的全组态大小。结果没有问题，效率也还行，但是大体系可能会压不住内存。
光靠GC确实还有不太行，可能用c++来写会好很一点。
"""
function m_config_size(ns::NuclearShell, Z::Integer, N::Integer)
    orbits = m_orbits(ns)
    orbit_number = msize(ns)

    p_m_map = zeros(Int64, max_pj(ns) + 1)
    n_m_map = zeros(Int64, max_nj(ns) + 1)

    for orb in orbits
        if orb.tz == -1
            p_m_map[m2index(orb.m)] += 1
        else
            n_m_map[m2index(orb.m)] += 1
        end
    end
    @assert binominal_safe(maximum(p_m_map), Z) "binomial of Int64 may overflow"
    @assert binominal_safe(maximum(n_m_map), N) "binomial of Int64 may overflow"

    hist_pM = Dict{Int, Int128}()
    p_ptn = partition(Z, p_m_map)
    for line in p_ptn
        num = Int128(1)
        m = 0
        for idx = 1:length(line)
            m += line[idx] * index2m(idx)
            num *= binomial(p_m_map[idx], line[idx])
        end
        hist_pM[m] = get(hist_pM, m, Int128(0)) + num
    end
    p_ptn = nothing
    GC.gc()

    hist_nM = Dict{Int, Int128}()
    n_ptn = partition(N, n_m_map)
    for line in n_ptn
        num = Int128(1)
        m = 0
        for idx = 1:length(line)
            m += line[idx] * index2m(idx)
            num *= binomial(n_m_map[idx], line[idx])
        end
        hist_nM[m] = get(hist_nM, m, Int128(0)) + num
    end
    n_ptn = nothing
    GC.gc()

    sum::BigInt = 0
    target = isodd(N + Z)
    for (pm, pnum) in hist_pM
        for (nm, nnum) in hist_nM
            if pm + nm == target
                sum += pnum * nnum
            end
        end
    end
    return sum
end

"""
    HO_shell(N::Integer)
构建一个`N = 2n + l`的谐振子壳层。
"""
function HO_shell(N::Integer)
    orbits = JOrbit[]
    for n = 0:div(N, 2)
        l = N - 2n
        for j in (l > 0 ? [2l+1, 2l-1] : [2l+1])
            push!(orbits, JOrbit(n, l, j, -1))
            push!(orbits, JOrbit(n, l, j, 1))
        end
    end
    NuclearShell(orbits)
end

"""
    HO_orbits(Nmax::Integer)
构建从`N = 0`到`N = Nmax`的所有谐振子轨道。
"""
function HO_orbits(Nmax::Integer)
    merge(ntuple(i -> HO_shell(i - 1), Nmax + 1)...)
end

"""
    struct ValenceSpace
        core_neutron::Int
        core_proton::Int
        ns::NuclearShell
    end
表示价空间的内。这个还是要有一定物理意义的，价轨道应该在核芯之上。
"""
struct ValenceSpace
    core_proton::Int  # 核芯的质子数
    core_neutron::Int # 核芯的中子数
    ns::NuclearShell  # 价轨道
    ValenceSpace(cp, cn, ns) = new(cp, cn, ns)
end

ValenceSpace(vs::ValenceSpace) = vs
"""
    ValenceSpace(core::Union{AbstractString, Isotope}, ns::Union{NuclearShell, Vector{JOrbit}})
`core`可以用一个同位素表示，比如``^{18}O``，价轨道可以用`NuclearShell`，也可以直接指定轨道。
"""
ValenceSpace(core::Union{AbstractString, Isotope}, ns::Union{NuclearShell, Vector{JOrbit}}) = begin
    iso = Isotope(core)
    ValenceSpace(iso.Z, iso.N, NuclearShell(ns))
end

"""
    Isotope(ns::ValenceSpace, Z::Integer, N::Integer)
给定一个价空间，给定价核子数目，得到核素。
"""
function Isotope(vs::ValenceSpace, Z::Integer, N::Integer)
    ns = vs.ns
    @assert Z < psize(ns) "Z = $Z is larger than valence size $(psize(ns))"
    @assert N < nsize(ns) "N = $N is larger than valence size $(nsize(ns))"
    Isotope(vs.core_proton + Z, vs.core_neutron + N)
end

"""
    Isotope(core::Union{AbstractString, Isotope}, Z::Integer, N::Integer)
给定一个核芯，以及在其之上的价核子数目，得到核素。这个函数没有价空间大小检查。
"""
function Isotope(core::Union{AbstractString, Isotope}, Z::Integer, N::Integer)
    iso = Isotope(core)
    Isotope(iso.Z + Z, iso.N + N)
end

"""
    valence(iso::Isotope, vs::ValenceSpace)
计算一个核素去掉核芯之后，价核子数目。
"""
function valence(iso::Isotope, vs::ValenceSpace)
    z = iso.Z - vs.core_proton
    n = iso.N - vs.core_neutron
    @assert z < psize(vs.ns) "valence proton number $z is larger than valence size $(psize(vs.ns))"
    @assert n < nsize(vs.ns) "valence neutron number $n is larger than valence size $(nsize(vs.ns))"
    return z, n
end

const p_space = ValenceSpace(2, 2, p_shell)
const sd_space = ValenceSpace(8, 8, sd_shell)
const pf_space = ValenceSpace(20, 20, pf_shell)


text:
@testset "Pressure" begin
    tolerances = Dict(Float16 => 1e-3, Float32 => 1e-5, Float64 => 1e-7, Real => 1e-7)
    @testset "Laplacian two dimensions" begin
        N, M = (4,4)
        typelist = [Float16 Float32 Float64 Real]
        for type in typelist
            testheight = zeros(type, (N,M))
            velocities = JuSwalbe.Twovector{Matrix{type}}(zeros(type, (N,M)),zeros(type, (N,M)))
            moment = JuSwalbe.Macroquant(ones(type, (N,M)), velocities, zeros(type, (N,M)), zeros(type, (N,M)))
            # Fill the array a little bit more demanding, otherwise the laplacian would be zero
            for i in 1:size(moment.height,1)
                for j in 1:size(moment.height,2)
                    moment.height[i,j] = i+j
                    testheight[i,j] = i+j
                end
            end
            # Make the laplace calculation
            test_lap1 = Δh(moment)
            @test isa(test_lap1, Array{type, 2})
            test_lap2 = Δh(testheight)
            @test isa(test_lap2, Array{type, 2})
            # Test the dimensions
            @test size(test_lap1) == (N,M)
            @test size(test_lap2) == (N,M)
            result = [8.0 4.0 4.0 0.0; 4.0 0.0 0.0 -4.0; 4.0 0.0 0.0 -4.0; 0.0 -4.0 -4.0 -8.0]
            resultcorrected = zeros(type, (N,M))
            for (ind, ele) in enumerate(result)
                resultcorrected[ind] = result[ind]
            end
            # Compare the results
            for j in 1:M, i in 1:N
                @test result[i,j] == test_lap1[i,j]
                @test result[i,j] == test_lap2[i,j]
            end
        end
    end

    @testset "Laplacian one dimension" begin
        N = 8
        typelist = [Float16 Float32 Float64 Real]
        for type in typelist
            testheight = ones(type, N)
            moment = JuSwalbe.Macroquant(ones(type, N), zeros(type, N), zeros(type, N), zeros(type, N))
            # Fill the array a little bit more demanding, otherwise the laplacian would be zero
            for i in 1:length(moment.height)
                moment.height[i] = i*i
                testheight[i] = i*i
            end
            # Make the laplace calculation
            test_lap1 = Δh(moment)
            @test isa(test_lap1, Array{type, 1})
            test_lap2 = Δh(testheight)
            @test isa(test_lap2, Array{type, 1})
            # Test the dimensions
            @test length(test_lap1) == N
            @test length(test_lap2) == N
            result = [66.0 2.0 2.0 2.0 2.0 2.0 2.0 -78.0]
            resultcorrected = zeros(type, N)
            for (ind, ele) in enumerate(result)
                resultcorrected[ind] = result[ind]
            end
            # Compare the results
            for i in 1:N
                @test result[i] == test_lap1[i]
                @test result[i] == test_lap2[i]
            end
        end
    end

    @testset "Disjoining pressure two dimensions" begin
        N, M = (4,4)

        typelist = [Float16 Float32 Float64]
        for type in typelist
            testheight = ones(type, (N,M)) .* reshape(collect(1:N*M),N,M)
            moment = JuSwalbe.Macroquant(testheight, JuSwalbe.Twovector(ones(type, (N,M)),ones(type, (N,M))), zeros(type, (N,M)), zeros(type, (N,M)))

            # Make the laplace calculation
            test_disj1 = Π(moment)
            @test isa(test_disj1, Array{type, 2})
            test_disj2 = Π(testheight)
            @test isa(test_disj1, Array{type, 2})
            # Test the dimensions
            @test size(test_disj1) == (N,M)
            @test size(test_disj2) == (N,M)
            result = [-1.6082e-5   -1.28656e-7  -2.20603e-8  -7.31997e-9;
                    -2.01025e-6  -7.44536e-8  -1.6082e-8   -5.86078e-9;
                    -5.95628e-7  -4.68862e-8  -1.20826e-8  -4.76503e-9;
                    -2.51281e-7  -3.14101e-8  -9.30669e-9  -3.92626e-9]
            # Compare the results
            resultcorrected = zeros(type, (N,M))
            for (ind, ele) in enumerate(result)
                resultcorrected[ind] = result[ind]
            end
            for i in 1:M, j in 1:N
                @test result[i,j] - test_disj1[i,j] ≈ 0 atol = tolerances[type]
                @test result[i,j] - test_disj2[i,j] ≈ 0 atol = tolerances[type]
            end
        end
    end
    @testset "Disjoining pressure one dimension" begin
        N = 4
        typelist = [Float16 Float32 Float64]
        for type in typelist
            testheight = ones(type, N) .* collect(1:N)
            moment = JuSwalbe.Macroquant(testheight, ones(type, N), zeros(type, N), zeros(type, N))

            # Make the laplace calculation
            test_disj1 = Π(moment)
            @test isa(test_disj1, Array{type, 1})
            test_disj2 = Π(testheight)
            @test isa(test_disj2, Array{type, 1})
            # Test the dimensions
            @test length(test_disj1) == N
            @test length(test_disj2) == N
            result = [-1.603e-5 -2.0e-6 -6.0e-7 -2.4e-7]

            # Compare the results
            resultcorrected = zeros(type, N)
            for (ind, ele) in enumerate(result)
                resultcorrected[ind] = result[ind]
            end
            for i in 1:N
                @test result[i] - test_disj1[i] ≈ 0 atol = tolerances[type]
                @test result[i] - test_disj2[i] ≈ 0 atol = tolerances[type]
            end
        end
    end

    @testset "Filmpressure nothing two dimensions" begin
        N, M = (4,4)
        typelist = [Float16 Float32 Float64]
        for type in typelist
            testheight = ones(type, (N,M))
            moment = JuSwalbe.Macroquant(testheight, JuSwalbe.Twovector(zeros(type ,(N,M)),zeros(type ,(N,M))), zeros(type ,(N,M)), zeros(type ,(N,M)))
            # Make the pressure calculation
            test_pressure1 = pressure(moment, γ=type(1.0), θ=zeros(type,(1,1)))
            test_pressure2 = pressure(moment, γ=type(1.0), θ=zeros(type,(1,1)))
            # Test the dimensions
            @test isa(test_pressure1, Array{type, 2})
            @test isa(test_pressure2, Array{type, 2})
            @test size(test_pressure1) == (N,M)
            @test size(test_pressure2) == (N,M)
            result = type(0.0)
            # Compare the results
            for j in 1:M, i in 1:N
                @test test_pressure1[i,j] == result
                @test test_pressure1[i,j] == result
            end
        end
    end

    @testset "Filmpressure nothing one dimension" begin
        N = 4
        typelist = [Float16 Float32 Float64]
        for type in typelist
            testheight = ones(type, N)
            moment = JuSwalbe.Macroquant(testheight, zeros(type ,N), zeros(type ,N), zeros(type ,N))
            # Make the pressure calculation
            test_pressure1 = pressure(moment, γ=type(1.0), θ=zeros(type, 1))
            test_pressure2 = pressure(moment, γ=type(1.0), θ=zeros(type, 1))
            # Test the dimensions
            @test isa(test_pressure1, Array{type, 1})
            @test isa(test_pressure2, Array{type, 1})
            @test length(test_pressure1) == N
            @test length(test_pressure2) == N
            result = type(0.0)
            # Compare the results
            for i in 1:N
                @test test_pressure1[i] == result
                @test test_pressure1[i] == result
            end
        end
    end

    @testset "Filmpressure nodisjoining two dimensional" begin
        N, M = (4,4)
        typelist = [Float16 Float32 Float64]
        for type in typelist
            testheight = ones(type, (N,M)) .* reshape(collect(1:N*M),N,M)
            moment = JuSwalbe.Macroquant(testheight, JuSwalbe.Twovector(zeros(type ,(N,M)),zeros(type ,(N,M))), zeros(type ,(N,M)), zeros(type ,(N,M)))

            test_pressure1 = pressure(moment, θ=zeros(type, (1,1)), γ=type(0.1))
            test_pressure2 = pressure(testheight, θ=zeros(type, (1,1)), γ=type(0.1))
            # Test the dimensions
            @test isa(test_pressure1, Array{type,2})
            @test isa(test_pressure2, Array{type,2})
            @test size(test_pressure1) == (N,M)
            @test size(test_pressure2) == (N,M)
            # Compare the results
            for i in 1:N, j in 1:M
                @test test_pressure1[i,j] - (type(-0.1)*Δh(moment))[i,j] ≈ type(0.0) atol = tolerances[type]
                @test test_pressure2[i,j] - (type(-0.1)*Δh(testheight))[i,j] ≈ type(0.0) atol = tolerances[type]
            end
        end
    end

    @testset "Filmpressure nodisjoining one dimensional" begin
        N = 4
        typelist = [Float16 Float32 Float64]
        for type in typelist
            testheight = ones(type, N) .* collect(1:N)
            moment = JuSwalbe.Macroquant(testheight, zeros(type ,N), zeros(type ,N), zeros(type ,N))

            test_pressure1 = pressure(moment, θ=zeros(type, 1), γ=type(0.1))
            test_pressure2 = pressure(testheight, θ=zeros(type, 1), γ=type(0.1))
            # Test the dimensions
            @test isa(test_pressure1, Array{type,1})
            @test isa(test_pressure2, Array{type,1})
            @test length(test_pressure1) == N
            @test length(test_pressure2) == N
            # Compare the results
            for i in 1:N
                @test test_pressure1[i] - (type(-0.1)*Δh(moment))[i] ≈ type(0.0) atol = tolerances[type]
                @test test_pressure2[i] - (type(-0.1)*Δh(testheight))[i] ≈ type(0.0) atol = tolerances[type]
            end
        end
    end

    @testset "Filmpressure NoLaplacian two dimensional" begin
        N, M = (4,4)
        typelist = [Float16 Float32 Float64]
        for type in typelist
            testheight = ones(type, (N,M))
            moment = JuSwalbe.Macroquant(testheight, JuSwalbe.Twovector(zeros(type ,(N,M)),zeros(type ,(N,M))), zeros(type ,(N,M)), zeros(type ,(N,M)))
            # Make the laplace calculation
            test_pressure1 = pressure(moment, γ=type(1.0))
            test_pressure2 = pressure(testheight, γ=type(1.0))
            # Test the dimensions
            @test isa(test_pressure1, Array{type,2})
            @test isa(test_pressure2, Array{type,2})
            @test size(test_pressure1) == (N,M)
            @test size(test_pressure2) == (N,M)
            # Compare the results
            for i in 1:N, j in 1:M
                @test test_pressure1[i,j] - Π(moment, γ=type(1.0))[i,j] ≈ type(0.0) atol = tolerances[type]
                @test test_pressure2[i,j] - Π(testheight, γ=type(1.0))[i,j] ≈ type(0.0) atol = tolerances[type]
            end
        end
    end

    @testset "Filmpressure NoLaplacian one dimensional" begin
        N = 4
        typelist = [Float16 Float32 Float64]
        for type in typelist
            testheight = ones(type, N)
            moment = JuSwalbe.Macroquant(testheight, zeros(type ,N), zeros(type ,N), zeros(type ,N))

            test_pressure1 = pressure(moment, γ=type(1.0))
            test_pressure2 = pressure(testheight, γ=type(1.0))
            # Test the dimensions
            @test isa(test_pressure1, Array{type,1})
            @test isa(test_pressure2, Array{type,1})
            @test length(test_pressure1) == N
            @test length(test_pressure2) == N
            # Compare the results
            for i in 1:N
                @test test_pressure1[i] - Π(moment, γ=type(1.0))[i] ≈ type(0.0) atol = tolerances[type]
                @test test_pressure2[i] - Π(testheight, γ=type(1.0))[i] ≈ type(0.0) atol = tolerances[type]
            end
        end
    end

    @testset "Filmpressure all two dimensional" begin
        N, M = (4,4)
        typelist = [Float16 Float32 Float64]
        for type in typelist
            testheight = ones(type, (N,M)) .* reshape(collect(1:N*M),N,M)
            moment = JuSwalbe.Macroquant(testheight, JuSwalbe.Twovector(zeros(type ,(N,M)),zeros(type ,(N,M))), zeros(type ,(N,M)), zeros(type ,(N,M)))

            # Make the laplace calculation
            test_pressure1 = pressure(moment, γ=type(1.0))
            test_pressure2 = pressure(testheight, γ=type(1.0))
            # Test the dimensions
            @test isa(test_pressure1, Array{type,2})
            @test isa(test_pressure2, Array{type,2})
            @test size(test_pressure1) == (N,M)
            @test size(test_pressure2) == (N,M)

            # Compare the results
            for j in 1:M, i in 1:N
                @test test_pressure1[i,j] - (type(-1.0)*Δh(moment)[i,j] + Π(moment, γ=type(1.0))[i,j]) ≈ type(0.0) atol = tolerances[type]
                @test test_pressure2[i,j] - (type(-1.0)*Δh(testheight)[i,j] + Π(testheight, γ=type(1.0))[i,j]) ≈ type(0.0) atol = tolerances[type]
            end
        end
    end

    @testset "Filmpressure all one dimensional" begin
        N = 4
        typelist = [Float16 Float32 Float64]
        for type in typelist
            testheight = ones(type, N) .* collect(1:N)
            moment = JuSwalbe.Macroquant(testheight, zeros(type ,N), zeros(type ,N), zeros(type ,N))

            # Make the laplace calculation
            test_pressure1 = pressure(moment, γ=type(1.0))
            test_pressure2 = pressure(testheight, γ=type(1.0))
            # Test the dimensions
            @test isa(test_pressure1, Array{type,1})
            @test isa(test_pressure2, Array{type,1})
            @test length(test_pressure1) == N
            @test length(test_pressure2) == N

            # Compare the results
            for i in 1:N
                @test test_pressure1[i] - (type(-1.0)*Δh(moment)[i] + Π(moment, γ=type(1.0))[i]) ≈ type(0.0) atol = tolerances[type]
                @test test_pressure2[i] - (type(-1.0)*Δh(testheight)[i] + Π(testheight, γ=type(1.0))[i]) ≈ type(0.0) atol = tolerances[type]
            end
        end
    end

    @testset "Pressure gradient h∇p two dimensions" begin
        N, M = (5,5)
        typelist = [Float16 Float32 Float64]
        for type in typelist
            testheight = ones(type, (N,M)) .* reshape(collect(1:N*M),N,M)
            testvel = JuSwalbe.Twovector(zeros(type, (N,M)), zeros(type, (N,M)))
            moment = JuSwalbe.Macroquant(height=testheight, velocity=testvel, pressure=testvel.x, energy=testvel.x)
            testforce = JuSwalbe.Forces(slip=testvel, thermal=testvel, h∇p=testvel, bathymetry=testvel)
            @test isa(testforce, JuSwalbe.Forces)
            @test isa(moment, JuSwalbe.Macroquant)
            testgrad = ∇p(moment, testforce)
            @test isa(testforce.h∇p, JuSwalbe.Twovector)
            @test size(testforce.h∇p.x) == (N,M)
            @test size(testforce.h∇p.y) == (N,M)
            for i in 1:N, j in 1:M
                @test testforce.h∇p.x[i,j] ≈ type(0.0) atol=tolerances[type]
                @test testforce.h∇p.y[i,j] ≈ type(0.0) atol=tolerances[type]
            end

            moment.pressure[3,3] = type(1.0)
            testgrad = ∇p(moment, testforce)
            # Upper row
            @test - testforce.h∇p.x[2,2] - moment.height[2,2]*type(1/12) ≈ type(0) atol=tolerances[type]
            @test - testforce.h∇p.x[2,3] - moment.height[2,3]*type(1/3) ≈ type(0) atol=tolerances[type]
            @test - testforce.h∇p.x[2,4] - moment.height[2,4]*type(1/12) ≈ type(0) atol=tolerances[type]
            # Lower row
            @test - testforce.h∇p.x[4,2] + moment.height[4,2]*type(1/12) ≈ type(0) atol=tolerances[type]
            @test - testforce.h∇p.x[4,3] + moment.height[4,3]*type(1/3) ≈ type(0) atol=tolerances[type]
            @test - testforce.h∇p.x[4,4] + moment.height[4,4]*type(1/12) ≈ type(0) atol=tolerances[type]
            # Left column
            @test - testforce.h∇p.y[2,2] - moment.height[2,2]*type(1/12) ≈ type(0) atol=tolerances[type]
            @test - testforce.h∇p.y[3,2] - moment.height[3,2]*type(1/3) ≈ type(0) atol=tolerances[type]
            @test - testforce.h∇p.y[4,2] - moment.height[4,2]*type(1/12) ≈ type(0) atol=tolerances[type]
            # Right column
            @test - testforce.h∇p.y[2,4] + moment.height[2,4]*type(1/12) ≈ type(0) atol=tolerances[type]
            @test - testforce.h∇p.y[3,4] + moment.height[3,4]*type(1/3) ≈ type(0) atol=tolerances[type]
            @test - testforce.h∇p.y[4,4] + moment.height[4,4]*type(1/12) ≈ type(0) atol=tolerances[type]
        end
    end
    @testset "Pressure gradient h∇p one dimension" begin
        N = 5
        typelist = [Float16 Float32 Float64]
        for type in typelist
            testheight = ones(type, N) .* collect(1:N)
            moment = JuSwalbe.Macroquant(height=testheight, velocity=testheight, pressure=testheight, energy=testheight)
            testforce = JuSwalbe.Forces(slip=testheight, thermal=testheight, h∇p=zeros(type,N), bathymetry=testheight)
            @test isa(testforce, JuSwalbe.Forces)
            @test isa(moment, JuSwalbe.Macroquant)
            testgrad = ∇p(moment, testforce)
            @test isa(testforce.h∇p, Array{type, 1})
            @test length(testforce.h∇p) == N

            moment.pressure = zeros(type, 5)
            moment.pressure[3] = type(1.0)
            testforce.h∇p = zeros(type, 5)
            testgrad = ∇p(moment, testforce)
            # Upper row
            @test - testforce.h∇p[1] + moment.height[1]*type(1/12) ≈ type(0) atol=tolerances[type]
            @test - testforce.h∇p[2] - moment.height[2]*type(2/3) ≈ type(0) atol=tolerances[type]
            @test - testforce.h∇p[3] + moment.height[3]*type(0) ≈ type(0) atol=tolerances[type]
            @test - testforce.h∇p[4] + moment.height[4]*type(2/3) ≈ type(0) atol=tolerances[type]
            @test - testforce.h∇p[5] - moment.height[5]*type(1/12) ≈ type(0) atol=tolerances[type]

        end
    end
end

max_stars_repo_name: Zitzeronion/JuSwalbe

max_issues_repo_name: Zitzeronion/JuSwalbe


text:
#== structure representing an Atom ==#
struct Atom
  atom_id::Int64
  symbol::String

  atom_center::SVector{3,Float64}
end
export Atom

#== structure representing a molecule ==#
struct Molecule
  atoms::Vector{Atom}
  mol_cxx::StdVector{JERI.Atom}
end
export Molecule

#== functions for Molecule ==#
function Base.getindex(mol::Molecule, index)
  return mol.atoms[index]
end

function Base.length(mol::Molecule)
  return length(mol.atoms)
end

function Base.iterate(mol::Molecule)
  return iterate(mol.atoms)
end

function Base.iterate(mol::Molecule, state)
  return iterate(mol.atoms, state)
end

function Base.push!(mol::Molecule, atom::Atom)
  return push!(mol.atoms, atom)
end


text:
export gather, gather!

"""
    gather!(dst, src, idx)

Reverse operation of [`scatter!`](@ref). Gathers data from source `src`
and writes it in destination `dst` according to the index array `idx`.
For each `k` in `CartesianIndices(idx)`, assign values to `dst` according to

    dst[:, ... , k] .= src[:, ... , idx[k]...]

Notice that if `idx` is a vector containing integers,
and both `dst` and `src` are matrices, previous expression simplifies to

    dst[:, k] .= src[:, idx[k]]

and `k` will run over `1:length(idx)`.

The elements of `idx` can be integers or integer tuples and may be repeated.
A single `src` column can end up being copied into zero, one,
or multiple `dst` columns.

See [`gather`](@ref) for an allocating version.
"""
function gather!(dst::AbstractArray{Tdst,Ndst},
                 src::AbstractArray{Tsrc,Nsrc},
                 idx::AbstractArray{Tidx, Nidx}) where
                    {Tdst, Tsrc, Ndst, Nsrc, Nidx, Tidx <: IntOrIntTuple}

    M = typelength(Tidx)
    d = Ndst - Nidx
    d == Nsrc - M || throw(ArgumentError("Incompatible input shapes."))
    size(dst)[1:d] == size(src)[1:d] || throw(ArgumentError("Incompatible input shapes."))
    size(dst)[d+1:end] == size(idx) || throw(ArgumentError("Incompatible input shapes."))

    colons = ntuple(i -> Colon(), d)
    for k in CartesianIndices(idx)
        view(dst, colons..., k) .= view(src, colons..., idx[k]...)
    end
    return dst
end

"""
    gather(src, idx) -> dst

Reverse operation of [`scatter`](@ref). Gathers data from source `src`
and writes it in a destination `dst` according to the index
array `idx`.
For each `k` in `CartesianIndices(idx)`, assign values to `dst`
according to

    dst[:, ... , k] .= src[:, ... , idx[k]...]

Notice that if `idx` is a vector containing integers
and `src` is a matrix, previous expression simplifies to

    dst[:, k] .= src[:, idx[k]]

and `k` will run over `1:length(idx)`.

The elements of `idx` can be integers or integer tuples and may be repeated.
A single `src` column can end up being copied into zero, one,
or multiple `dst` columns.

See [`gather!`](@ref) for an in-place version.
"""
function gather(src::AbstractArray{Tsrc, Nsrc},
                idx::AbstractArray{Tidx, Nidx}) where
                    {Tsrc, Nsrc, Nidx, Tidx<:IntOrIntTuple}

    M = typelength(Tidx)
    dstsize = (size(src)[1:Nsrc-M]..., size(idx)...)
    dst = similar(src, Tsrc, dstsize)
    return gather!(dst, src, idx)
end

# Simple implementation with getindex for integer array.
# Perf equivalent to the one above (which can also handle the integer case)
# leave it here to show the simple connection with getindex.
function gather(src::AbstractArray{Tsrc, Nsrc},
                idx::AbstractArray{<:Integer}) where {Tsrc, Nsrc}
    colons = ntuple(i -> Colon(), Nsrc-1)
    return src[colons..., idx]
end


text:
using Test
using NeXLUncertainties

@testset "UncertainValues" begin
    @testset "Label" begin
        l1a = label(123)
        l1b = label(123)
        l2a = nl"X"
        l2b = nl"X"
        l3a = label(123.0)
        l3b = label(123.0)
        l3c = label(123.1)
        @test isequal(l1a, l1b)
        @test isequal(l2a, l2b)
        @test isequal(l3a, l3b)
        @test !isequal(l1a, l2a)
        @test !isequal(l1a, l3a)
        @test !isequal(l2a, l3a)
        @test !isequal(l3a, l3c)
    end

    @testset "Covariances" begin

        cov = [
            0.1^2 -0.3*0.1*0.2 0.0
            -0.3*0.1*0.2 0.2^2 0.2*0.2*0.3
            0.0 0.2*0.2*0.3 0.3^2
        ]
        @test checkcovariance!(cov)

        bad = [
            0.1^2 -0.29*0.1*0.2 0.0
            -0.3*0.1*0.2 0.2^2 0.2*0.2*0.3
            0.0 0.2*0.2*0.3 0.3^2
        ]

        @test_throws ErrorException checkcovariance!(bad)

        lbls = [label("X0"), label("X1"), label("Z")]
        vals = [2.0, 8.0, 12.0]

        uvs1 = uvs(lbls, vals, cov)

        @test σ(uvs1, label("X0")) == 0.100
        @test σ(uvs1, label("X1")) == 0.200
        @test σ(uvs1, label("Z")) == 0.3
        @test σ(uvs1, label("Z")) == 0.3

        @test filter(uvs1, [label("Z"), label("X1")]) ==
              [cov[3, 3] cov[3, 2]; cov[2, 3] cov[2, 2]]

        @test uvs1[label("Z")] == uv(12.0, 0.3)
        @test uvs1[label("X1")] == uv(8.0, 0.2)
        @test get(uvs1, label("X1"), uv(8.1, 0.0)) == uv(8.0, 0.2)
        @test get(uvs1, label("Z"), uv(8.1, 0.0)) ≠ uv(8.0, 0.2)
        @test get(uvs1, label("Zz"), uv(8.1, 0.0)) == uv(8.1, 0.0)

        @test length(uvs1) == 3

        @test value(uvs1, label("Z")) == 12.0
        @test value(uvs1, label("X1")) == 8.0

        @test covariance(uvs1, label("X1"), label("Z")) == 0.2 * 0.2 * 0.3
        @test covariance(uvs1, label("Z"), label("X1")) == 0.2 * 0.2 * 0.3

        @test variance(uvs1, label("Z")) == 0.3^2
        @test variance(uvs1, label("X1")) == 0.2^2

        @test uncertainty(uvs1, label("Z")) == 0.3
        @test uncertainty(uvs1, label("X1")) == 0.2
        @test uncertainty(uvs1, label("Z"), 2.0) == 2.0 * 0.3
        @test uncertainty(uvs1, label("X1"), 3.0) == 3.0 * 0.2
    end

    @testset "cat" begin
        uvs1 = uvs(
            [nl"A", nl"B", nl"C"],
            [1.0, 2.0, 3.0],
            [0.1 0.05 -0.03; 0.05 0.2 0.06; -0.03 0.06 0.3],
        )
        uvs2 = uvs(
            [nl"D", nl"E", nl"F"],
            [4.0, 5.0, 6.0],
            [0.2 0.1 -0.13; 0.1 0.3 -0.12; -0.13 -0.12 0.1],
        )
        uvs3 = cat([uvs1, uvs2])
        @test value(uvs3, nl"A") == 1.0
        @test value(uvs3, nl"C") == 3.0
        @test value(uvs3, nl"A") == value(uvs1, nl"A")
        @test value(uvs3, nl"C") == value(uvs1, nl"C")
        @test value(uvs3, nl"E") == value(uvs2, nl"E")
        @test value(uvs3, nl"F") == value(uvs2, nl"F")
        @test covariance(uvs3, nl"A", nl"B") == 0.05
        @test covariance(uvs3, nl"E", nl"E") == variance(uvs3, nl"E")
        @test covariance(uvs3, nl"E", nl"E") == 0.3
        @test covariance(uvs3, nl"E", nl"F") == -0.12
        @test covariance(uvs3, nl"B", nl"D") == 0.0
        @test covariance(uvs3, nl"E", nl"E") == variance(uvs3, nl"E")
        @test σ(uvs3, nl"E") == sqrt(0.3)
        @test σ(uvs2, nl"E") == σ(uvs3, nl"E")

        ex = filter(uvs3, [nl"A", nl"C", nl"E"])
        @test ex[1, 1] == covariance(uvs3, nl"A", nl"A")
        @test ex[1, 2] == covariance(uvs3, nl"A", nl"C")
        @test ex[1, 2] == ex[2, 1]
        @test ex[3, 1] == covariance(uvs3, nl"E", nl"A")
        @test ex[1, 3] == ex[3, 1]

        @test uvs3[nl"E"] == UncertainValue(5.0, sqrt(0.3))
        @test uvs3[nl"C"] == UncertainValue(3.0, sqrt(0.3))
        @test get(uvs3, nl"C", missing) == UncertainValue(3.0, sqrt(0.3))
        @test ismissing(get(uvs3, nl"XX", missing))
        @test get(uvs3, nl"C", missing) == uvs3[nl"C"]
        @test get(uvs3, nl"C", UncertainValue(3.2, sqrt(0.99))) == uvs3[nl"C"]
        @test get(uvs3, nl"zZ", UncertainValue(3.2, sqrt(0.99))) ==
              UncertainValue(3.2, sqrt(0.99))

        @test length(uvs3) == 6
        @test length(uvs1) == 3

        vals = values(uvs3)
        lbls = labels(uvs3)
        @test value(uvs3[lbls[1]]) == vals[1]
        @test value(uvs3[lbls[3]]) == vals[3]

        @test uncertainty(uvs3, nl"A", 2.0) == 2.0 * sqrt(0.1)
    end
end


text:
export ProductFunction

struct ProductFunction{T, N} <: AbstractMathFun{T, N}
    terms::MFCollection

    ProductFunction{T, N}() where {T, N} = new{T, N}(MFCollection(undef, 0))
    ProductFunction{T, N}(terms::MFCollection) where {T, N} = new{T, N}(terms)

    Base.insert!(p::ProductFunction, i::Integer, term::AbstractMathFun) = insert!(p.terms, i, term)
    Base.push!(p::ProductFunction, term::AbstractMathFun) = push!(p.terms, term)
    Base.deleteat!(p::ProductFunction, range) = deleteat!(p.terms, range)
    Base.length(p::ProductFunction) = length(p.terms)
    Base.getindex(p::ProductFunction, i::Integer) = getindex(p.terms, i)
    Base.setindex!(p::ProductFunction, val::AbstractMathFun, i::Integer) = setindex!(p.terms, val, i)
    Base.iterate(p::ProductFunction) = Base.iterate(p.terms)
    Base.iterate(p::ProductFunction, state) = Base.iterate(p.terms, state)

    function Functions.string(p::ProductFunction, emitter::FunctionEmitter)
        str = "∏("
        first_print = true
        for term in p.terms
            term_str = string(term, emitter)
            length(term_str) != 0 && (str *= (first_print ? term_str : ", " * term_str); first_print = false)
        end
        str * ")"
    end

    function (p::ProductFunction)(call::MFCall)
        product = 1.0
        for term in p.terms
            product *= term(call::MFCall)
        end
        call(product)
    end

    function Functions.emit(f::ProductFunction{T, N}, emitter::FunctionEmitter{T}) where {T, N}
        prod_expr = Expr(:call, :∏)
        for term in f.terms
            push!(prod_expr.args, emit(term, emitter))
        end
        push!(prod_expr.args, Expr(:kw, :T, Symbol(string(T))))
        prod_expr
    end

    function Functions.simplify(f::ProductFunction{T, N}, emit_constants::Bool) where {T, N}
        constprod = 1
        termcpy = MFCollection(undef, 0)

        for item in f.terms
            item = simplify(item, emit_constants)
            if is_constant(item, emit_constants)
                constprod *= item.value
            else
                if item isa ProductFunction
                    append!(termcpy, item.terms)
                else
                    push!(termcpy, item)
                end
            end
        end

        (constprod != 1 && constprod != 0) && push!(termcpy,  ConstFunction{T}(constprod))
        (length(termcpy) == 0) && return ConstFunction{T}(constprod)
        length(termcpy) == 1 && return termcpy[1]
        return ProductFunction{T, N}(termcpy)
    end
end


text:
# ------------------------------------------------------------------------------
# F R E Q U E N T   P A T T E R N   A N A L Y S I S
#
# Using FP Growth
#
# Watch https://www.youtube.com/watch?v=VB8KWm8MXss for a walk-through of the
# algorithm and https://www.youtube.com/watch?v=ToswH_dA7KU for a walk-through
# of how to mine the rules.
#
# TODO Research which packages exist and still have active support
# TODO Implement minimal memory management where needed by assigning objects to
#      nothing and running the garbage collector – gc().
# TODO Make it possible to select restrictions when building the FP tree - e.g.
#      provide a minimum support level required. This will knock out items
#      early in the process and save subsequent calculations, but also put a
#      limit on the rules that can be mined afterwards.
# ------------------------------------------------------------------------------

# Load packages ----------------------------------------------------------------

using CSV
using DataFrames
using Statistics

# Define structures ------------------------------------------------------------
# Building on the code from:
# https://stackoverflow.com/questions/36593490/tree-data-structure-in-julia

mutable struct TreeNode
    parent::Int
    children::Vector{Int}

    # item_id allows two nodes to represent the same item at different locations
    item_id::Int
    item_label::String
    count::Int
end

struct Tree
    nodes::Vector{TreeNode}
end

Tree() = Tree([TreeNode(0, Vector{Int}())])

function addchild(tree::Tree, id::Int)
    1 <= id <= length(tree.nodes) || throw(BoundsError(tree, id))
    push!(tree.nodes, TreeNode(id, Vector{}()))
    child = length(tree.nodes)
    push!(tree.nodes[id].children, child)
    child
end

children(tree, id) = tree.nodes[id].children

parent(tree,id) = tree.nodes[id].parent

# ------------------------------------------------------------------------------
# T e s t   r u n s
# Using data from the Groceries dataset in the arules package for R
# ------------------------------------------------------------------------------

# Read sample data: groceries
filename = "groceries.tsv"
filepath = joinpath(@__DIR__, filename)
#println(filepath)

groceries = CSV.read(filepath, delim = ",")

# ------------------------------------------------------------------------------
# Step 1  F o r m a t   t h e   d a t a
#
# Organise the data into a tidy format with transaction ID and Item as the two
# variables.
# ------------------------------------------------------------------------------

function tidy_data(df)
    tidy_df = df
    return tidy_df
end

# ------------------------------------------------------------------------------
# Step 2  C a l c u l a t e   s u p p o r t
# ------------------------------------------------------------------------------

function calculate_support(df)
    # Takes as input a dataframe where each row is a transaction and each
    # column an item. Values are assumed boolean.
    map(eachcol(df)) do col
        mean(col)
    end
end

# This code is more concise than the function above :D
mean.(eachcol(groceries))

# ------------------------------------------------------------------------------
# Step 3  R e o r d e r   i t e m s e t s
# ------------------------------------------------------------------------------

function reorder_items(df)
    # Get support for all columns: support
    support = calculate_support(df)
    n_items = length(support)
    item_names = names(groceries)

    # Build dataframe of support and id: df_ordered
    df_ordered = DataFrame()
    df_ordered.support = support
    df_ordered.id = 1:n_items
    df_ordered.label = item_names
    sort!(df_ordered, rev = true)
end

# Test function
# sorted = reorder_items(groceries)
# sorted = nothing

# ------------------------------------------------------------------------------
# Step 4  R e o r d e r   t r a n s a c t i o n s
# ------------------------------------------------------------------------------

#
function reorder_transactions(df)
    # Get dimensions of dataframe
    n_transactions = nrow(df)
    n_items = ncol(df)

    # Get support for all items
    sorted = reorder_items(df)

    # Add transaction id (:tid) to df
    insertcols!(df, 1, tid = 1:n_transactions)

    # Stack the columns of the dataframe to make it tidy
    df_long = stack(df, 2:(n_items + 1))

    filter!(row -> row[:value] == true, df_long)
    rename!(df_long, :variable => :label)
    tidy = join(df_long, sorted, on = :label)
    sort!(tidy, [:tid, order(:support, rev = true)])

    return tidy
end

# Test function
# tidy_df = reorder_transactions(groceries)

# ------------------------------------------------------------------------------
# Step 5  B u i l d   t i d y   F P   t r e e
# ------------------------------------------------------------------------------

connections = DataFrame(parent = Int[],
                        child = Int[])
nodes = DataFrame(id = Int[],
                  item_id = Int[],
                  item_label = String[],
                  count = Int[])

# Function to add sub trees to a node. Takes a vector of items to add which
# allows the same function to be used for adding single nodes.
function add_sub_tree!(tree::Tree, id::Int, itemset::DataFrame)
  # Set active node to id
  current_node = tree.nodes[id]

  # Cycle through all items in itemset
  for item in itemset
      # Append id.children with item
      append!(current_node.children, item)

      # Append new node with item as id
      push!(tree.nodes,
            TreeNode(current_node.item_id, Vector{Int}(), item, "", 1))

      # Make new node the current node
      current_node = tree.nodes[item]
  end
end

# Function to build the FP tree
function build_fp_tree(df)
    # df should be in a tidy format with columns for transaction id, item id,
    # item label and support.

    # Reorder the transactions using reorder_transactions()

    # Initialise the tree (build the root node)
    fp_tree = Tree([TreeNode(0, Vector{Int}(), 0, "root", 0)])

    # Cycle through all transactions
    for transaction in transactions

        # Cycle through all items in transaction
        for item in itemset
            # Start at root node
            # Lookup current item in children
            # If success increase count and move to relevant child node
            # If fail add current and remaining transactions as sub tree
        end
    end # for

    # Return FP tree
end

# ------------------------------------------------------------------------------
# Step 5  B u i l d   F P   t r e e
# ------------------------------------------------------------------------------

# Function to add sub trees to a node. Takes a vector of items to add which
# allows the same function to be used for adding single nodes.
function add_sub_tree!(tree::Tree, id::Int, itemset::Vector{Int})
    # Set active node to id
    current_node = tree.nodes[id]

    # Cycle through all items in itemset
    for item in itemset
        # Append id.children with item
        append!(current_node.children, item)

        # Append new node with item as id
        push!(tree.nodes,
              TreeNode(current_node.item_id, Vector{Int}(), item, "", 1))

        # Make new node the current node
        current_node = tree.nodes[item]
    end
end

# fp_tree.nodes.item_id
#
# current_node = fp_tree.nodes[1]
# append!(current_node.children, 16)
# current_node.children
#
# push!(fp_tree.nodes, TreeNode(1, Vector{Int}(), 16, "inner", 1))
#
# [fp_tree.nodes[i].item_id .== 16 for i = 1:2]
#
# (fp->fp_tree.nodes).(1:2)
#
# getfield.(fp_tree.nodes, item_id)
#
# dump(fp_tree)

# Function to build the FP tree
function build_fp_tree(df)
    # df should be in a tidy format with columns for transaction id, item id,
    # item label and support.

    # Reorder the transactions using reorder_transactions()

    # Initialise the tree (build the root node)
    fp_tree = Tree([TreeNode(0, Vector{Int}(), 0, "root", 0)])

    # Cycle through all transactions
    for transaction in transactions

        # Cycle through all items in transaction
        for item in itemset
            # Start at root node
            # Lookup current item in children
            # If success increase count and move to relevant child node
            # If fail add current and remaining transactions as sub tree
        end
    end # for

    # Return FP tree
end

# ------------------------------------------------------------------------------
# Step 6a  B u i l d   c o n d i t i o n a l   p a t t e r n   b a s e
# ------------------------------------------------------------------------------

function build_conditional_pattern_base(args)
    body
end

# ------------------------------------------------------------------------------
# Step 6b  B u i l d   c o n d i t i o n a l   F P   t r e e
# ------------------------------------------------------------------------------

function build_conditional_fp_tree(args)
    body
end

# ------------------------------------------------------------------------------
# Step 6c  G e n e r a t e   f r e q u e n t   p a t t e r n s
# ------------------------------------------------------------------------------

function generate_frequent_patterns(args)
    body
end

# ------------------------------------------------------------------------------
# Q u e r y   f u n c t i o n s
# ------------------------------------------------------------------------------

function query_rules(df, support, confidence, lift)
    # Filter the dataframe to meet the criterias for support, confidence and
    # lift: rules
    return rules
end

max_stars_repo_name: DopedCow/FPGrowth

max_issues_repo_name: DopedCow/FPGrowth

max_forks_repo_name: DopedCow/FPGrowth


text:
using BlackBoxOptim: fitness, index, tag, rank_by_fitness!, bbsetup
using BlackBoxOptim: hassolution, setsolution!, solution

# The ask tell interface aims to allow flexibility in cases where
# the evaluator might not be easy/possible to define up-front
# or when the optimization needs to be "driven" from the outside
# rather than from BBO itself.
@testset "Top-level ask/tell interface" begin
    @testset "Single float fitness, with intermediate solutions" begin

        # A typical use case for this is that we need to map parameter sets to
        # intermediate "solutions" before we can evaluate their fitness. We want
        # to avoid the repeated generation of solutions since that might be costly.
        # In this example it is not but in general it might be, and then we want
        # to "cache" the intermediate solutions with the parameters.
        params2solution(params) = string(sum(params))
        fitn(sol) = length(sol)

        # Since the user will step through the ask/tell cycle themselves we should
        # not require knowledge of the optimization function / evaluator. For cases
        # where the fitness is not a single Float64 one should state the fitness
        # scheme explicitly as for the multi-objective case?
        oc = bbsetup(; SearchRange = (-1.0, 1.0), NumDimensions = 5)

        # When we ask for solutions we get either solutions for params
        # we previously supplied or empty Nullables. The alternative would
        # be to embed the solutions in the Candidate struct and let user
        # have getters/setters for them. Seems unnecessary to add this
        # level of complexity. Easier to have a general top-level interface
        # where we can get and set solutions/phenotypes associated with each
        # set of parameters. OTOH we will probably save the solutions in the
        # candidates anyway so users that require a more control can just ask
        # to get the candidates themselves and then return them (sorted).
        # Here we show the simpler, first version. It will unpack the params
        # and the solutions from the (internal) candidates and thus does not
        # require the user to know anything about the Candidate type or internals
        # of BBO.
        params, solutions = ask(oc; withSolutions = true)

        @test length(params) > 0
        @test length(params) == length(solutions)

        # In this case there are no solutions yet but in general there might be.
        # Our task is now to sort the params and solutions based on their fitness.
        fs = Array{Float64}()
        for i in eachindex(params)
            p = params[i]
            s = solutions[i]
            @test typeof(p) <: Vector{Float64}
            @test typeof(s) <: Nullable
            @test fitness(c) == NaN       # Since not yet evaluated (and since we did not specify a non-standard fitness scheme).
            @test isnull(s)               # Since not yet evaluated. If one of the params is one for which we have previously supplied the solution is should be returned (in a Nullable) instead.

            ns = isnull(s) ? params2solution(p) : s # Don't calc new solution is already given before.
            push!(newsolutions, ns)
            push!(fs, fitn(ns))
        end

        # A usability risk here is that the user only sorts one of the arrays
        # returned? Another risk is that the fitnesses might not be comparable on
        # a global scale, i.e. between calls to ask/tell. This is likely to mess up
        # tracing and archives as well as the saving of the "best" individuals. OTOH,
        # since optimization is "driven" from the outside maybe the resposibility for
        # keeping history, tracing etc should be on the user rather than on BBO?
        p = sortperm(fs)
        tell!(oc, params[p], fs[p], newsolutions[p])
    end
end


text:
# Group Interface

Base.one(G::PermGroup) = Permutation(Perm(degree(G)), G)
GroupsCore.order(::Type{T}, G::AbstractPermutationGroup) where {T} =
    order(T, StabilizerChain(G))
GroupsCore.gens(G::PermGroup) = Permutation.(gens_raw(G), Ref(G))

function Base.rand(
    rng::Random.AbstractRNG,
    rs::Random.SamplerTrivial{Gr}
) where {Gr <: PermGroup}
    G = rs[]
    tr = transversals(G)
    # Ref(rng) is needed for julia-1.3
    img =rand.(Ref(rng), tr)
    return perm_by_baseimages(G, img)
end
### iteration protocol for PermGroups

Base.eltype(::Type{GT}) where {I,GT<:PermGroup{I}} = Permutation{I,GT}
Base.IteratorSize(::Type{<:AbstractPermutationGroup}) = Base.HasLength()

struct PermGroupIter{V, I, S, W}
    base_images::V
    itr::I
    state::S
    tmp_word::W
end

function Base.iterate(G::PermGroup)
    itr = base_images(G)
    tmp_word = Word(eltype(eltype(G))[])
    img, st = iterate(itr)
    permiter = PermGroupIter(img, itr, st, tmp_word)
    g = perm_by_baseimages(G, img, false, permiter.tmp_word)
    return g, (permiter, 1)
end

function Base.iterate(G::PermGroup, state)
    permiter, count = state
    count >= length(G) && return nothing
    img, st = iterate(permiter.itr, permiter.state)
    g = perm_by_baseimages(G, img, false, permiter.tmp_word)
    return g, (permiter, count+1)
end

# GroupElement Interface
Base.parent(g::Permutation) = g.parent

Base.:(==)(g::Permutation, h::Permutation) = parent(g) === parent(h) && g.perm == h.perm

Base.deepcopy_internal(g::Permutation, stackdict::IdDict) =
    Permutation(Base.deepcopy_internal(g.perm, stackdict), parent(g))

Base.inv(g::Permutation) = Permutation(inv(g.perm), parent(g))

function Base.:(*)(g::Permutation, h::Permutation)
    parent(g) === parent(h) ||
        error("Cannot multiply elements from different permutation groups")
    return Permutation(g.perm * h.perm, parent(g))
end

# IO

function Base.show(io::IO, G::PermGroup)
    init = isdefined(G, :stabchain) ? " of order $(order(StabilizerChain(G)))" : ""
    ngen = length(gens(G))

    print(io, "Permutation group on ", ngen, " generator", ngen > 1 ? "s" : "", init)
end

function Base.show(io::IO, ::MIME"text/plain", G::PermGroup)
    init = isdefined(G, :stabchain) ? " of order $(order(StabilizerChain(G)))" : ""
    ngen = length(gens(G))

    println(
        io,
        "Permutation group on ",
        ngen,
        " generator",
        ngen > 1 ? "s" : "",
        init,
        " generated by",
    )
    Base.print_array(io, gens(G))
end

Base.show(io::IO, g::Permutation) = show(io, perm(g))

### Performance optimizations

Base.similar(p::Permutation{T}, ::Type{S} = T) where {T,S} =
    Permutation(similar(perm(p), T), parent(p))
Base.isone(p::AbstractPerm) = nfixedpoints(p) == degree(p)
Base.:(^)(p::Permutation, n::Integer) = parent(p)(p.perm^n)
GroupsCore.order(::Type{T}, p::Permutation) where {T} = order(T, p.perm)
Base.hash(p::Permutation, h::UInt) = hash(typeof(p), hash(p.perm, hash(parent(p), h)))

# Mutable API

function GroupsCore.one!(g::AbstractPerm)
    for i in eachindex(g)
        g[i] = i
    end
    return g
end

Base.@propagate_inbounds function GroupsCore.mul!(
    out::AbstractPerm,
    g::AbstractPerm,
    h::AbstractPerm,
)
    out = (out === h ? similar(out) : out)
    @boundscheck @assert degree(out) >= max(degree(g), degree(h))

    @inbounds for i in eachindex(out)
        out[i] = (i^g)^h #h[g[i]]
    end

    return out
end

#### end of Group interface

# accessors
gens_raw(G::PermGroup) = G.gens

# AbstractPerm Interface??
"""
	perm(g::AbstractPerm)::Perm
Return an instance of `Perm`, a parent-less permutation backed by simple vector storage.
"""
perm(g::Permutation) = g.perm
perm(g::AbstractPerm) = Perm([i^g for i in eachindex(g)], false)

Base.setindex!(g::Permutation, v::Integer, n::Integer) = g.perm[n] = v

@doc doc"""
	degree(G::PermGroup)
Return the degree of `G`, i.e. the length of the storage of permutations in `G`.
"""
degree(G::PermGroup) = return G.deg
degree(G::SymmetricGroup) = G.n
degree(p::Permutation) = degree(parent(p))
degree(p::Perm{I}) where {I} = I(length(p.d))

permtype(p::Permutation) = permtype(perm(p))
Base.sign(p::Permutation) = sign(perm(p))


text:
function cancel_notification(params::CancelParams, server::LanguageServerInstance, conn)
end

function setTrace_notification(params::SetTraceParams, server::LanguageServerInstance, conn)
end

# TODO Provide type for params
function setTraceNotification_notification(params, server::LanguageServerInstance, conn)
end

function julia_getCurrentBlockRange_request(tdpp::VersionedTextDocumentPositionParams, server::LanguageServerInstance, conn)
    fallback = (Position(0, 0), Position(0, 0), tdpp.position)
    uri = URI2(tdpp.textDocument.uri)

    hasdocument(server, uri) || return nodocuemnt_error(uri)

    doc = getdocument(server, uri)

    if doc._version !== tdpp.version
        return mismatched_version_error(uri, doc, tdpp, "getCurrentBlockRange")
    end

    offset = get_offset(doc, tdpp.position)
    x = getcst(doc)
    loc = 0

    if headof(x) === :file
        for (i, a) in enumerate(x)
            if loc <= offset <= loc + a.fullspan
                # we'll try to use the next expression instead if the current expr is a
                # NOTHING or we're in the whitespace after an expr _and_ that whitespace wraps
                # to a new line
                if !(loc <= offset <= loc + a.span) || headof(a) === :NOTHING
                    if length(x) > i
                        thisline, _ = get_position_at(doc, loc + a.span)
                        if tdpp.position.line > thisline || headof(a) === :NOTHING
                            loc += a.fullspan
                            a = x[i + 1]
                        end
                    end
                end

                if a.head === :macrocall && a.args[1].head === :globalrefdoc && length(a.args) == 4 && CSTParser.defines_module(a.args[4])
                    for i in 1:3
                        loc += a.args[i].fullspan
                    end
                    a = a.args[4]
                end

                if CSTParser.defines_module(a) # Within module at the top-level, lets see if we can select on of the block arguments
                    if loc <= offset <= loc + a.trivia[1].span # Within `module` keyword, so return entire expression
                        return Position(get_position_at(doc, loc)...), Position(get_position_at(doc, loc + a.span)...), Position(get_position_at(doc, loc + a.fullspan)...)
                    end
                    if loc + a.trivia[1].fullspan <= offset <= loc + a.trivia[1].fullspan + a.args[2].span # Within name of the module, so return entire expression
                        return Position(get_position_at(doc, loc)...), Position(get_position_at(doc, loc + a.span)...), Position(get_position_at(doc, loc + a.fullspan)...)
                    end

                    if loc + a.trivia[1].fullspan + a.args[2].fullspan <= offset <= loc + a.trivia[1].fullspan + a.args[2].fullspan + a.args[3].span # Within the body of a module
                        loc += a.trivia[1].fullspan + a.args[2].fullspan
                        for b in a.args[3].args
                            if loc <= offset <= loc + b.span
                                return Position(get_position_at(doc, loc)...), Position(get_position_at(doc, loc + b.span)...), Position(get_position_at(doc, loc + b.fullspan)...)
                            end
                            loc += b.fullspan
                        end
                    elseif loc + a.trivia[1].fullspan + a.args[2].fullspan + a.args[3].fullspan < offset <= loc + a.trivia[1].fullspan + a.args[2].fullspan + a.args[3].fullspan + a.trivia[2].span # Within `end` of the module, so return entire expression
                        return Position(get_position_at(doc, loc)...), Position(get_position_at(doc, loc + a.span)...), Position(get_position_at(doc, loc + a.fullspan)...)
                    end
                else
                    return Position(get_position_at(doc, loc)...), Position(get_position_at(doc, loc + a.span)...), Position(get_position_at(doc, loc + a.fullspan)...)
                end
            end
            loc += a.fullspan
        end
    end
    return fallback
end

function julia_activateenvironment_notification(params::NamedTuple{(:envPath,),Tuple{String}}, server::LanguageServerInstance, conn)
    server.env_path = params.envPath

    trigger_symbolstore_reload(server)
end

julia_refreshLanguageServer_notification(_, server::LanguageServerInstance, conn) =
    trigger_symbolstore_reload(server)


text:
abstract type Enumerable end

Base.IteratorSize(::Type{T}) where {T <: Enumerable} = Base.SizeUnknown()
IteratorInterfaceExtensions.isiterable(x::Enumerable) = true

haslength(S) = Base.IteratorSize(S) isa Union{Base.HasLength, Base.HasShape} ? Base.HasLength() : Base.IteratorSize(S)


text:
import RNG: AbstractRNG

"""
```julia
WrappedRNG{R, T1, T2} <: AbstractRNG{T2}
WrappedRNG(base_rng, T2)
WrappedRNG(R, T2, args...)
```

Wrap a RNG which originally provides output in T1 into a RNG that provides output in T2.

# Examples
```jldoctest
julia> r = Xorshifts.Xorshift128Star(123);

julia> RNG.output_type(r)
UInt64

julia> r1 = WrappedRNG(r, UInt32);

julia> RNG.output_type(r1)
UInt32

julia> r2 = WrappedRNG(Xorshifts.Xorshift128Star, UInt32, 123);

julia> RNG.output_type(r2)
UInt32

julia> @Test.test rand(r1, UInt32, 3) == rand(r2, UInt32, 3)
Test Passed
  Expression: rand(r1,UInt32,3) == rand(r2,UInt32,3)
   Evaluated: UInt32[0x18a21796,0x20241598,0x63c65407] == UInt32[0x18a21796,0x20241598,0x63c65407]
```
"""
type WrappedRNG{R<:AbstractRNG, T1<:BitTypes, T2<:BitTypes} <: AbstractRNG{T2}
    base_rng::R
    x::T1
    p::Int
    function WrappedRNG(r, x, p)
        @assert T1 ≠ T2
        wr = new(r, x, p)
        if sizeof(T1) > sizeof(T2)
            wr.x = rand(wr.base_rng, T1)
        end
        wr
    end
end
WrappedRNG{T1<:BitTypes, T2<:BitTypes}(
    base_rng::AbstractRNG{T1}, ::Type{T2}) = WrappedRNG{typeof(base_rng), T1, T2}(base_rng, 0, 0)
function WrappedRNG{R<:AbstractRNG, T2<:BitTypes}(::Type{R}, ::Type{T2}, args...)
    base_rng = R(args...)
    WrappedRNG{R, output_type(base_rng), T2}(base_rng, 0, 0)
end

@inline rand{R<:AbstractRNG, T1<:BitTypes}(rng::WrappedRNG{R, T1}, ::Type{T1}) = rand(rng.base_rng, T1)

@inline function rand{R<:AbstractRNG, T1<:BitTypes, T2<:BitTypes}(rng::WrappedRNG{R, T1, T2}, ::Type{T2})
    s1 = sizeof(T1)
    s2 = sizeof(T2)
    if s2 >= s1
        t = rand(rng.base_rng, T1) % T2
        for i in 2:(s2 ÷ s1)
            t |= (rand(rng.base_rng, T1) % T2) << ((s1 << 3) * (i - 1))
        end
    else
        t = rng.x % T2
        rng.p += 1
        if rng.p == s1 ÷ s2
            rng.p = 0
            rng.x = rand(rng.base_rng, T1)
        else
            rng.x >>= s2 << 3
        end
    end
    return t
end


text:
using Random
using Parameters
#using Debugger

Random.seed!(1)

abstract type TreeSearch end

@with_kw mutable struct Mcts <: TreeSearch
	depth::Int64 = 12
	iters::Int64 = 100
	tMax::Int64 = 200

	c::Float64 = 1.0

	mdp = nothing
	nnet = nothing
	reuseTree::Bool = false
	pi0 = nothing

	#Tree = Set()
	#Nsa = Dict()
	#Ns = Dict()
	#Q = Dict()

	Tree = Set{Array{Float64,1}}()
	Nsa = Dict{Tuple{Array{Float64,1}, Float64}, Float64}()
	Ns = Dict{Array{Float64,1}, Float64}()
	Q = Dict{Tuple{Array{Float64,1}, Float64}, Float64}()
end

function init!(m::Mcts;
               mdp = nothing,
			   piRollout = nothing,
			   depth::Int64 = 12,
			   iters::Int64 = 100,
			   explorConst::Float64 = 1.0,
			   tMaxRollouts::Int64 = 200,
			   reuseTree::Bool = true,
			   nnet = nothing)
	m.depth = depth
	m.iters = iters
	m.c = explorConst
	m.tMax = tMaxRollouts
	m.reuseTree = reuseTree
	m.mdp = mdp
	m.nnet = nnet
	m.pi0 = piRollout

	resetTree!(m)
end

function resetTree!(m::Mcts)
	m.Tree = Set{Array{Float64,1}}()
	m.Nsa = Dict{Tuple{Array{Float64,1}, Float64}, Float64}()
	m.Ns = Dict{Array{Float64,1}, Float64}()
	m.Q = Dict{Tuple{Array{Float64,1}, Float64}, Float64}()

	sizehint!(m.Tree, m.iters)
	sizehint!(m.Nsa, m.iters)
	sizehint!(m.Ns, m.iters)
	sizehint!(m.Q, m.iters)
end

function act(m::Mcts, state::Array{Float64,1}, obstacles::Array{Tuple{Float64, Float64}, 1})::Float64
	if m.reuseTree
		resetTree!(m)
	end
	action = selectAction(m, state, m.depth)
	return action
end

function selectAction(m::Mcts, s::Array{Float64,1}, d::Int64)::Float64
	for i in 1:m.iters
		simulate(m, s, d, m.pi0)
	end
	q, action = maximum([(m.Q[(s,a)], a) for a in actions(m.mdp, s)])
	#println("q=$q")
	return action
end

function simulate(m::Mcts, s::Array{Float64,1}, d::Int64, pi0)::Float64
	if isEnd(m.mdp, s)>0 #[1]
		return 0
	end
	if d == 0 # we stop exploring the tree, just estimate Qval here
		return rollout(m, s, m.tMax, pi0)
	end
	if !in(s, m.Tree)
		for a in actions(m.mdp, s)
			m.Nsa[(s,a)], m.Ns[s], m.Q[(s,a)] =  0, 1, 0.
		end
		push!(m.Tree, s)
		return rollout(m, s, m.tMax, pi0)
	end
	a::Float64 = maximum([(m.Q[(s,a)]+m.c*sqrt(log(m.Ns[s])/(1e-5 + m.Nsa[(s,a)])), a) for a in actions(m.mdp, s)])[2]
	sp::Vector{Float64}, r::Float64, _ = sampleSuccReward(m.mdp, s, a)
	q::Float64 = r + discount(m.mdp) * simulate(m, sp, d-1, pi0)
	m.Nsa[(s,a)] += 1
	#println(typeof(s))
	m.Ns[s] += 1
	m.Q[(s,a)] += (q-m.Q[(s,a)])/m.Nsa[(s,a)]
	return q
end

function rollout(m::Mcts, s::Array{Float64,1}, d::Int64, pi0)::Float64
	if d==0 || isEnd(m.mdp, s)>0#[1]
		return 0
	else
		a::Float64 = pi0(m.mdp, s)
		sp::Vector{Float64}, r::Float64, _ = sampleSuccReward(m.mdp, s, a)
		return r + discount(m.mdp) * rollout(m, sp, d-1, pi0)
	end
end


text:
@system Density begin
    weather ~ hold

    PD0: initial_planting_density => 55 ~ preserve(u"m^-2", parameter)

    CDSF: cold_damage_shape_factor => 0.9 ~ preserve(u"K^-1", parameter)
    CDCT: cold_damage_critical_temperature => -15 ~ preserve(u"°C", parameter)

    enable_cold_damage => true ~ flag(parameter)

    "mortality due to cold damage (2019-08-09: KHM, KDY)"
    CDM(s=CDSF, Tc=CDCT, T=weather.T_air): cold_damage_mortality => begin
        x = exp(-s * (T - Tc))
        x / (1 + x)
    end ~ track(u"d^-1", when=enable_cold_damage)

    CDS(CDS, CDM): cold_damage_survival_rate => (CDS * -CDM) ~ accumulate(init=1, min=0)

    PD(PD0, CDS): planting_density => (PD0 * CDS) ~ track(u"m^-2")
end


text:
using Zygote: @adjoint, gradient, @adjoint!

@adjoint setdiff(args...) = setdiff(args...), _ -> nothing
@adjoint intersect(args...) = intersect(args...), _ -> nothing

nograd(f, args...) = f(args...)
@adjoint nograd(f, args...) = f(args...), _ -> nothing
@adjoint optimaltree(args...) = optimaltree(args...), _ -> nothing

@adjoint get_size_dict(args...) = get_size_dict(args...), _ -> nothing
@adjoint leg_analysis(args...) = leg_analysis(args...), _ -> nothing

@adjoint! function copyto!(xs::AbstractArray, ys::AbstractArray)
  xs_ = copy(xs)
  copyto!(xs, ys), function (dxs)
    copyto!(xs_, xs)
    (nothing, dxs)
  end
end


text:
using Base.Test
using Tensors

# Single mode QR

t = rand([Mode(k,2) for k in 1:4])
q,r = qr(t, 3)
@test_approx_eq_eps(norm(t - q*r), 0, 1e-12)

# Multiple modes QR

t = rand([Mode(k,2) for k in 1:4])
q,r = qr(t, Any[3,1], 5)
@test_approx_eq_eps(norm(t - q*r), 0, 1e-12)

# Single mode SVD

t = rand([Mode(k,2) for k in 1:4])
u,s,v = svd(t, 3)
@test_approx_eq_eps(norm(t - u*scale(s,v)), 0, 1e-12)

# Multiple modes SVD

t = rand([Mode(k,2) for k in 1:4])
u,s,v = svd(t, Any[4,2], 5)
@test_approx_eq_eps(norm(t - scale(u,s)*v), 0, 1e-12)


text:
type TrainState{n,m}
    ss::NTuple{n, Array{Float64}}   # holds the layers' output
    as::NTuple{n, Array{Float64}}   # holds the activation
    ds::NTuple{n, Array{Float64}}   # temporary storage for gradient
    dws::NTuple{n, Matrix{Float64}} # holds the gradient of weight matrices
    dbs::NTuple{n, Vector{Float64}} # holds the gradient of bias vector
    ls::NTuple{Int}
end

# creates a TrainState adapted to a given Neural net and # of samples
function TrainState(nn::NN, nsamples::Int) # nsamples = 3
    siz = size(nn)
    ss  = tuple([Array(Float64, s, nsamples) for s in siz]...)
    # dws = tuple( Array(Float64, siz[1], inputsize(nn)),
    #             [Array(Float64, siz[i], siz[i-1]) for i in 2:length(siz)]...)
    dws = tuple( map(similar, nn.ws)... )
    dbs = tuple( map(similar, nn.bs)... )

    TrainState{depth(nn),nsamples}(ss,
                                   map(similar,ss),
                                   map(similar,ss),
                                   dws,
                                   dbs,
                                   siz)
end

compatible(nn::NN, ss::TrainState) = size(nn) == size(ss)

size(ss::TrainState) = ss.ls
nsamples{n,m}(ss::TrainState{n,m}) = m

function cycle{T<:FloatingPoint}(nn::NN,
    input::StridedArray{T},
    output::StridedArray{T},
    metric::PreMetric,
    weights::StridedVector{T})

    cycle!(TrainState(nn, size(input,2)), nn, input, output, metric, weights)
end

function cycle!{n,m,T<:FloatingPoint}(s::TrainState{n,m},
    net::NN{n},
    input::StridedArray{T},
    output::StridedArray{T},
    metric::PreMetric,
    weights::StridedVector{T})

    @assert compatible(net,s) "incompatible network and state"
    @assert size(input,1) == inputsize(net) "incompatible size of input and network"
    @assert size(output,1) == size(net.ws[end],1) "incompatible size of output and network"
    @assert size(output,2) == size(input,2) "incompatible # samples between input and output"
    if ndims(input) == 1
        @assert m == 1 "incompatible # of samples in input"
    else
        @assert m == size(input,2) "incompatible # of samples in input"
    end

    # forward pass
    prev = input
    for i in 1:depth(net)
        A_mul_B!(s.as[i], net.ws[i], prev)
        # bias term
        for j in 1:size(s.as[i],1)
            bias = net.bs[i][j]
            for k in 1:size(s.as[i],2)
                s.as[i][j,k] += bias
            end
        end
        nfunc!(s.ss[i], s.as[i], net.ns[i])
        prev = s.ss[i]
    end

    # backward pass
    dmetric!(s.ds[end], s.ss[end], output, metric, weights)
    for i in depth(net):-1:2
        dnfunc!(s.ds[i], s.as[i], s.ss[i], net.ns[i])
        sum!(s.dbs[i], s.ds[i])                    # gradient on bias term
        A_mul_Bt!( s.dws[i],   s.ds[i], s.ss[i-1]) # gradient on weight matrix
        At_mul_B!(s.ds[i-1], net.ws[i], s.ds[i])
    end
    dnfunc!(s.ds[1], s.as[1], s.ss[1], net.ns[1])
    sum!(s.dbs[1], s.ds[1])                        # gradient on bias term
    A_mul_Bt!( s.dws[1],   s.ds[1], input)

    s
end

function calc{T<:FloatingPoint}(net::NN, input::StridedArray{T})
    @assert size(input,1) == inputsize(net) "incompatible size of input and network"

    # forward pass
    prev = input
    for i in 1:depth(net)
        prev = net.ws[i] * prev
        # bias term
        for j in 1:size(prev,1)
            bias = net.bs[i][j]
            for k in 1:size(prev,2)
                prev[j,k] += bias
            end
        end
        nfunc!(prev, prev, net.ns[i])
    end

    prev
end

max_stars_repo_name: fredo-dedup/SimpleNeuralNets

max_issues_repo_name: fredo-dedup/SimpleNeuralNets

max_forks_repo_name: fredo-dedup/SimpleNeuralNets


text:
using ThreeBodyDecay
using Test

Ps = ThreeBodyParities('+','-','+'; P0='+')

@testset "Three Body Parities structure" begin
    #
    @test Ps.P1 == Ps[1] == '+'
    @test Ps.P2 == Ps[2] == '-'
    @test Ps.P3 == Ps[3] == '+'
    @test Ps.P0 == Ps[0] == Ps[4] == '+'
    #
    @test_throws ErrorException Ps[5]
    @test_throws ErrorException ThreeBodyParities(0,1,0; two_js0=1)
    @test_throws ErrorException ThreeBodyParities('+','+','+')
end

@testset "operations and interate" begin
    @test sum(i == '+' for i in Ps) == 3
end

max_stars_repo_name: mmikhasenko/ThreeBodyDecay

max_issues_repo_name: mmikhasenko/ThreeBodyDecay

max_forks_repo_name: mmikhasenko/ThreeBodyDecay


text:
"""
    master_query(master_df::DataFrame, key_type::Symbol, key_value::Any)

Wrapper for a query, takes in an already loaded DataFrame `master_df`, a `key_type` to
search over (e.g. `obsid`), and a `key_value` to find (e.g. `0123456789`)

Returns the full row for any observations matching the search criteria
"""
function master_query(master_df::DataFrame, key_type::Symbol, key_value::Any)
    observations = filter(row -> row[key_type] == key_value, master_df)

    if size(observations, 1) == 0
        @warn "master_query returned no results for $key_type with $key_value search"
    end

    if size(observations, 1) == 1
        # Return DataFrameRow for single observation result
        observations = observations[1, :]
    end

    return observations
end

function master_query(mission::Mission, key_type::Symbol, key_value::Any)
    return master_query(master(mission), key_type, key_value)
end

function master_query_public(master_df::DataFrame, key_type::Symbol, key_value::Any)
    query_result = master_query(master_df, key_type, key_value)

    if typeof(query_result) == DataFrames.DataFrameRow
        # Single row returned by master query
        if query_result[:time] >= now() # time is in the future
            query_result = DataFrame() # empty dataframe returned
        end
    else
        query_result = filter(x->x[:time]<=now(), query_result)
    end

    return query_result
end

function master_query_public(mission::Mission, key_type::Symbol, key_value::Any)
    return master_query_public(master(mission), key_type, key_value)
end

function master_query_public(mission::Mission)
    return filter(x->x[:time]<=now(), master(mission))
end


text:
struct SynchronousExecution{T <: AbstractFloat,
                            A <: AbstractVector,
                            S <: LQSolver} <: AbstractExecution
    subworkers::Vector{SubWorker{T,A,S}}

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

function init_subproblems!(ph::AbstractProgressiveHedgingSolver, subsolver::QPSolver, execution::SynchronousExecution)
    return init_subproblems!(ph, subsolver, execution.subworkers)
end

function iterate!(ph::AbstractProgressiveHedgingSolver, ::SynchronousExecution)
    return iterate_nominal!(ph)
end

function init_workers!(::AbstractProgressiveHedgingSolver, ::SynchronousExecution)
    return nothing
end

function close_workers!(::AbstractProgressiveHedgingSolver, ::SynchronousExecution)
    return nothing
end

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

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

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

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

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

function fill_submodels!(ph::AbstractProgressiveHedgingSolver, scenarioproblems, execution::SynchronousExecution)
    return fill_submodels!(ph, scenarioproblems, execution.subworkers)
end
# API
# ------------------------------------------------------------
struct Synchronous <: Execution end

function (execution::Synchronous)(::Type{T}, ::Type{A}, ::Type{S}) where {T <: AbstractFloat, A <: AbstractVector, S <: LQSolver}
    return SynchronousExecution(T,A,S)
end

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