text:
### TEST FOR THE GROUND STATE OF ELLIPTICAL WAVEGUIDE FOR DIFFERENT ECCENTRICITY
## PLOTTING RESULTS IN NORMALIZED COORDINATES AS Salasnich - SciPost Physics (2022)

using SolitonBEC
using Printf
using Plots
using Elliptic
hbar = 6.62607015e-34 / (2 * pi)

# --------- Numerics ---------
L = 100e-6

num = Numerics(
  50e-3, #T
  5e-5, #dt
  L, #S
  L * 1e-3, #ds
)

# --------- Simulation ---------
khaykovich_gpe = Simulation(
  "NPSE",
  "ellipse",
)

#--------- Potential ---------
function elliptical(ecc::Float64)
  r = Potential(
    "ellipse", # type
    5e-6, #width
    15e-6, #position
    5e-6, #energy
    ecc, #ϵ
    L / E(2 * pi, ecc) #a
  )
  return r
end

# --------- Apparata ---------
std_apparatus = Apparatus(
  6.941 * 1.660539e-27, #m (conversion AMU -> kg)
  0, # as
  2 * pi * 710, # ω_perp
  4e3, #N
  1.77e-11, #γ
  2 * pi * 4,# ω_z
)

# --------- InitialStates ---------
InitialState1 = InitialState(
  "sech", #type
  1.7e-6, # width
  0e6 #v0
)

## Configurations
configs = []
for ecc in [0, 0.5, 0.75, 0.9]
  potential = elliptical(ecc)
  push!(configs, (num, khaykovich_gpe, potential, std_apparatus, InitialState1))
end

display(configs)
mem_limit = 15000000000 #byte
cnt = 1

pyplot()
plt_width = 800
plt_height = 600
## Curvature potential--------------------------------------------

fig1 = plot(title="",
  xlabel="s/L",
  ylabel="U_Q/E_a",
  size=(plt_width, plt_height))

for (num, sim, pot, app, state) in configs
  characteristic_energy = hbar^2 / (app.m * pot.a^2)

  space_steps = Int(floor(num.S / num.ds))
  space = LinRange(-num.S / 2, num.S / 2, space_steps)
  plot!(fig1, space / L, -abs.(map(s -> potential(sim, app, pot, s), space)) / characteristic_energy, reuse=false,
    label="ϵ = $(pot.ϵ)")
end

## Non Interacting ground state--------------------------------------------

fig2 = plot(title="",
  xlabel="s/L",
  ylabel="ρ(s) * a",
  reuse=false,
  size=(plt_width, plt_height))

for (num, sim, pot, app, state) in configs
  global cnt
  @printf("\n------Running simulation # %3i \n", cnt)
  estimate = mem_estimate(num)

  if estimate < mem_limit

    @printf("Estimated memory consumption: %4.1f MiB\n", estimate / (1024^2))

    # Natural orientation choice for plot
    coeffs = get_coefficients(sim, app, pot, state)

    time_steps = Int(floor(num.T / num.dt))
    space_steps = Int(floor(num.S / num.ds))
    time = LinRange(0, num.T, time_steps)
    space = LinRange(-num.S / 2, num.S / 2, space_steps)

    time, space, ψ = @time ground_state_solve(num, coeffs)
    plot!(fig2, space / L,
      abs.(ψ) .^ 2 * pot.a,
      label="ϵ = $(pot.ϵ)")

    display(fig1)
    display(fig2)
  else
    @printf("Estimated memory consumption (%4.1f MiB) exceed maximum!\n", estimate / (1024^2))
  end

  cnt += 1

end

## Spontaneous Symmetry Breaking--------------------------------------------
configs_gamma = []

fixed_potential = elliptical(0.9)
mass = 6.941 * 1.660539e-27
ω_perp = 2 * pi * 710
Num = 4e3

as_apparatus(as) = Apparatus(
  mass, #m (conversion AMU -> kg)
  as, # as
  ω_perp, # ω_perp
  Num, #N
  1.77e-11, #γ
  2 * pi * 4,# ω_z
)
l_perp = sqrt(hbar / (mass * ω_perp))

for gamma in [-1, -1.5, -2, -3, 2, 3]
  as = gamma * l_perp^2 / (2 * fixed_potential.a * (Num - 1))
  display(as_apparatus(as))
  push!(configs_gamma, (num, khaykovich_gpe, fixed_potential, as_apparatus(as), InitialState1))
end
p = Plots.palette(:rainbow_bgyr_35_85_c72_n256, length(configs_gamma) + 3)

fig3 = plot(title="",
  xlabel="s/L",
  ylabel="ρ(s) * a",
  reuse=false,
  size=(plt_width, plt_height),
  palette=p)

for (num, sim, pot, app, state) in configs_gamma
  global cnt
  @printf("\n------Running simulation # %3i \n", cnt)
  estimate = mem_estimate(num)

  if estimate < mem_limit

    @printf("Estimated memory consumption: %4.1f MiB\n", estimate / (1024^2))

    # Natural orientation choice for plot
    coeffs = get_coefficients(sim, app, pot, state)

    time_steps = Int(floor(num.T / num.dt))
    space_steps = Int(floor(num.S / num.ds))
    time = LinRange(0, num.T, time_steps)
    space = LinRange(-num.S / 2, num.S / 2, space_steps)

    time, space, ψ = @time ground_state_solve(num, coeffs)
    γ = (2 * pot.a * app.as * (app.N - 1)) / (l_perp^2)
    plot!(fig3, space / L,
      abs.(ψ) .^ 2 * pot.a,
      label="γ = $(γ)")

    display(fig3)
  else
    @printf("Estimated memory consumption (%4.1f MiB) exceed maximum!\n", estimate / (1024^2))
  end

  cnt += 1

end

max_stars_repo_name: lorenzifrancesco/soliton-BEC

max_issues_repo_name: lorenzifrancesco/soliton-BEC

max_forks_repo_name: lorenzifrancesco/soliton-BEC


text:
using DimensionalData, Test, Unitful
using DimensionalData.LookupArrays, DimensionalData.Dimensions
using DimensionalData.LookupArrays: _slicespan, isrev, _bounds
using DimensionalData.Dimensions: _slicedims

@testset "isrev" begin
    @test isrev(ForwardOrdered()) == false
    @test isrev(ForwardOrdered()) == false
end

@testset "reverse" begin
    @test reverse(ForwardOrdered()) == ReverseOrdered()
    @test reverse(ReverseOrdered()) == ForwardOrdered()
    @test reverse(Unordered()) == Unordered()
    lu = Sampled(order=ForwardOrdered(), span=Regular(1))
    @test order(reverse(lu)) == ReverseOrdered()
    lu = Categorical(order=ReverseOrdered())
    @test order(reverse(lu)) == ForwardOrdered()
end

@testset "getindex" begin
    ind = [10.0, 20.0, 30.0, 40.0, 50.0]

    @testset "Irregular forwards" begin
        m = Sampled(ind, order=ForwardOrdered(), span=Irregular((10.0, 60.0)), sampling=Intervals(Start()))
        @test bounds(getindex(m, 3:3)) == (30.0, 40.0)
        @test bounds(getindex(m, 1:5)) == (10.0, 60.0)
        @test bounds(getindex(m, 2:3)) == (20.0, 40.0)
        m = Sampled(ind, ForwardOrdered(), Irregular((0.0, 50.0)), Intervals(End()), NoMetadata())
        @test bounds(getindex(m, 3:3)) == (20.0, 30.0)
        @test bounds(getindex(m, 1:5)) == (0.0, 50.0)
        @test bounds(getindex(m, 2:3)) == (10.0, 30.0)
        m = Sampled(ind, ForwardOrdered(), Irregular((5.0, 55.0)), Intervals(Center()), NoMetadata())
        @test bounds(getindex(m, 3:3)) == (25.0, 35.0)
        @test bounds(getindex(m, 1:5)) == (5.0, 55.0)
        @test bounds(getindex(m, 2:3)) == (15.0, 35.0)
    end

    @testset "Irregular reverse" begin
        revind = [50.0, 40.0, 30.0, 20.0, 10.0]
        m = Sampled(revind; order=ReverseOrdered(), span=Irregular(0.0, 50.0), sampling=Intervals(Start()))
        @test bounds(getindex(m, 1:5)) == (0.0, 50.0)
        @test bounds(getindex(m, 1:2)) == (30.0, 50.0)
        @test bounds(getindex(m, 2:3)) == (20.0, 40.0)
        m = Sampled(revind, ReverseOrdered(), Irregular(10.0, 60.0), Intervals(End()), NoMetadata())
        @test bounds(getindex(m, 1:5)) == (10.0, 60.0)
        @test bounds(getindex(m, 1:2)) == (40.0, 60.0)
        @test bounds(getindex(m, 2:3)) == (30.0, 50.0)
        m = Sampled(revind, ReverseOrdered(), Irregular(0.5, 55.0), Intervals(Center()), NoMetadata())
        @test bounds(getindex(m, 1:5)) == (0.5, 55.0)
        @test bounds(getindex(m, 1:2)) == (35.0, 55.0)
        @test bounds(getindex(m, 2:3)) == (25.0, 45.0)
    end

    @testset "Irregular with no bounds" begin
        m = Sampled(ind, ForwardOrdered(), Irregular(nothing, nothing), Intervals(Start()), NoMetadata())
        @test bounds(getindex(m, 3:3)) == (30.0, 40.0)
        @test bounds(getindex(m, 2:4)) == (20.0, 50.0)
        # TODO should this be built into `identify` to at least get one bound?
        @test bounds(getindex(m, 1:5)) == (10.0, nothing)
        m = Sampled(ind, ForwardOrdered(), Irregular(nothing, nothing), Intervals(End()), NoMetadata())
        @test bounds(getindex(m, 3:3)) == (20.0, 30.0)
        @test bounds(getindex(m, 2:4)) == (10.0, 40.0)
        @test bounds(getindex(m, 1:5)) == (nothing, 50.0)
        m = Sampled(ind, ForwardOrdered(), Irregular(nothing, nothing), Intervals(Center()), NoMetadata())
        @test bounds(getindex(m, 3:3)) == (25.0, 35.0)
        @test bounds(getindex(m, 2:4)) == (15.0, 45.0)
        @test bounds(getindex(m, 1:5)) == (nothing, nothing)
    end

end

@testset "bounds" begin

    @testset "Intervals" begin
        @testset "Regular bounds are calculated from interval type and span value" begin
            @testset "forward ind" begin
                ind = 10.0:10.0:50.0
                dim = X(Sampled(ind, order=ForwardOrdered(), sampling=Intervals(Start()), span=Regular(10.0)))
                @test bounds(dim) == (10.0, 60.0)
                dim = X(Sampled(ind, order=ForwardOrdered(), sampling=Intervals(End()), span=Regular(10.0)))
                @test bounds(dim) == (0.0, 50.0)
                dim = X(Sampled(ind, ForwardOrdered(), Regular(10.0), Intervals(Start()), NoMetadata()))
                @test bounds(dim) == (10.0, 60.0)
                dim = X(Sampled(ind, ForwardOrdered(), Regular(10.0), Intervals(End()), NoMetadata()))
                @test bounds(dim) == (0.0, 50.0)
                dim = X(Sampled(ind, ForwardOrdered(), Regular(10.0), Intervals(Center()), NoMetadata()))
                @test bounds(dim) == (5.0, 55.0)
            end
            @testset "reverse ind" begin
                revind = [10.0, 9.0, 8.0, 7.0, 6.0]
                dim = X(Sampled(revind, ; order=ReverseOrdered(), sampling=Intervals(Start()), span=Regular(-1.0)))
                dim = X(Sampled(revind, ReverseOrdered(), Regular(-1.0), Intervals(Start()), NoMetadata()))
                @test bounds(dim) == (6.0, 11.0)
                dim = X(Sampled(revind, ; order=ReverseOrdered(), sampling=Intervals(End()), span=Regular(-1.0)))
                dim = X(Sampled(revind, ReverseOrdered(), Regular(-1.0), Intervals(End()), NoMetadata()))
                @test bounds(dim) == (5.0, 10.0)
                dim = X(Sampled(revind, ; order=ReverseOrdered(), sampling=Intervals(Center()), span=Regular(-1.0)))
                dim = X(Sampled(revind, ReverseOrdered(), Regular(-1.0), Intervals(Center()), NoMetadata()))
                @test bounds(dim) == (5.5, 10.5)
            end
        end
        @testset "Irregular bounds are whatever is stored in span" begin
            ind = 10.0:10.0:50.0
            dim = X(Sampled(ind, ForwardOrdered(), Irregular(0.0, 50000.0), Intervals(Start()), NoMetadata()))
            @test bounds(dim) == (0.0, 50000.0)
            @test bounds(getindex(dim, 2:3)) == (20.0, 40.0)
        end
        @testset "Explicit bounds are is stored in span matrix" begin
            ind = 10.0:10.0:50.0
            bnds = vcat(ind', (ind .+ 10)')
            dim = X(Sampled(ind, ForwardOrdered(), Explicit(bnds), Intervals(Start()), NoMetadata()))
            @test bounds(dim) == (10.0, 60.0)
            @test bounds(_slicedims(getindex, dim, 2:3)[1][1]) == (20.0, 40.0)
        end
    end

    @testset "Points" begin
        ind = 10:15
        dim = X(Sampled(ind; order=ForwardOrdered(), sampling=Points()))
        @test bounds(dim) == (10, 15)
        ind = 15:-1:10
        dim = X(Sampled(ind; order=ReverseOrdered(), sampling=Points()))
        last(dim), first(dim)
        @test bounds(dim) == (10, 15)
        dim = X(Sampled(ind; order=Unordered(), sampling=Points()))
        @test bounds(dim) == (nothing, nothing)
    end

    @testset "Categorical" begin
        ind = [:a, :b, :c, :d]
        dim = X(Categorical(ind; order=ForwardOrdered()))
        @test order(dim) == ForwardOrdered()
        @test_throws ErrorException step(dim)
        @test span(dim) == NoSpan()
        @test sampling(dim) == NoSampling()
        @test dims(lookup(dim)) === nothing
        @test locus(dim) == Center()
        @test bounds(dim) == (:a, :d)
        dim = X(Categorical(ind; order=ReverseOrdered()))
        @test bounds(dim) == (:d, :a)
        @test order(dim) == ReverseOrdered()
        dim = X(Categorical(ind; order=Unordered()))
        @test bounds(dim) == (nothing, nothing)
        @test order(dim) == Unordered()
    end

end

@testset "dims2indices with Transformed" begin
    tdimz = Dim{:trans1}(Transformed(identity, X())),
            Dim{:trans2}(Transformed(identity, Y())),
            Z(NoLookup(1:1))
    @test dims2indices(tdimz, (X(1), Y(2), Z())) == (1, 2, Colon())
    @test dims2indices(tdimz, (Dim{:trans1}(1), Dim{:trans2}(2), Z())) == (1, 2, Colon())
end

ext: jl


text:
__precompile__()
module Firelight

using Sockets
using Electron
using Electron: URI_file
import HTTP
import JSON
# import FunctionalCollections

"""
    @LOCAL(filespec)

Construct an absolute URI to `filespec` relative to the source file containing the macro call.
"""
macro LOCAL(filespec)
    # v0.7: base = String(__source__.file)
    #       filespec isa String && return URI_file(base, filespec) # can construct eagerly
    #       return :(URI_file($base, $(esc(filespec))))
    return :(URI_file(@__DIR__, $(esc(filespec))))
end

application() = Electron.default_application()
#let app
#    global function application()
#        if !@isdefined(app) || !app.exists
#            app = Application()
#        end
#        return app::Application
#    end
#end

function escapehtml(io::IO, bytes::AbstractVector{UInt8}, attribute::Bool=true)
    for b in bytes
        if b < 0x20
            write(io, "&#$b;")
        elseif b == UInt8('&')
            write(io, "&")
        elseif b == UInt8('<')
            write(io, "<")
        elseif b == UInt8('>')
            write(io, ">")
        elseif attribute && b == UInt8('\'')
            write(io, "'") # for html4, use #39
        elseif attribute && b == UInt8('"')
            write(io, """)
        else
            write(io, b)
        end
    end
end

Base.position(io::IOContext) = position(io.io)

include("dom.jl")

mutable struct ServerState
    key::UInt64
    content::String
    shadowdom::Vector{Node}
    reason::String
    exists::Bool
    function ServerState(key::UInt64, reason::String)
        state = new(key, "", Vector{Node}(), reason, true)
        sessions[key] = state
        return state
    end
end

global sessions = Dict{UInt64, ServerState}()

function start_server()
    process_id = getpid()
    local server_name
    id = UInt(1)
    while true
        server_name = Electron.generate_pipe_name("firelight-$process_id-$id")
        ispath(server_name) || break
        id += 1
    end
    let server = listen(server_name)
        @schedule HTTP.listen(server_name, tcpref=Ref(server)) do request::HTTP.Request
            key = tryparse(UInt64, HTTP.header(request, "Session-id", ""), 16)
            key === nothing && return HTTP.Response(400, "Session-id header missing")
            session = get(sessions, key, nothing)
            session === nothing && return HTTP.Response(410, "Session non-existant")
            return invokelatest(render, request, session)
        end
    end
    return server_name
end

function invokelatest(f, args...)
    return eval(Expr(:body, Expr(:return, Expr(:call, Core._apply, QuoteNode(f), QuoteNode(args)))))
end

function start_session()
    key = rand(UInt64)
    while haskey(sessions, key)
        key += 1
    end
    session = ServerState(key, "ans")
    return session
end

function start_session(server_name::String)
    app = application()
    session = start_session()
    keystr = hex(session.key, 16)
    w = Window(app, @LOCAL("views/main.html?session-id=$keystr&server=$(HTTP.escapeuri(server_name))"))
    return session
end

function end_session(key::UInt64)
    session = pop!(sessions, key)
    session.exists = false
    nothing
end

function render(req::HTTP.Request, session::ServerState)
    rep = req.response
    HTTP.setheader(rep, "Content-Type" => "text/plain;charset=utf-8")
    target = req.target
    if target == "/"
        rep.body = session.content
    elseif target == "/fork"
        session = start_session()
        keystr = hex(session.key, 16)
        rep.body = keystr
    elseif target == "/reason"
        rep.body = session.reason
    elseif ismatch(r"/select/(\d+)", target)
        id = tryparse(Int, match(r"/select/(\d+)", target)[1])
        if id === nothing || id < 1 || id > length(session.shadowdom)
            return HTTP.Response(404, "Id not found in shadowdom")
        end
        o = session.shadowdom[id].object
        Core.eval(Main, Expr(:block, Expr(:(=), :sel, QuoteNode(o)), Expr(:return, nothing)))
    elseif ismatch(r"/dump/(\d+)", target)
        id = tryparse(Int, match(r"/dump/(\d+)", target)[1])
        if id === nothing || id < 1 || id > length(session.shadowdom)
            return HTTP.Response(404, "Id not found in shadowdom")
        end
        o = session.shadowdom[id].object
        rep.body = sprint(o) do io, r
            dump(io, r, maxdepth=1)
        end
    else
        return HTTP.Response(404)
    end
    return rep
end

end # module

max_stars_repo_name: vtjnash/Firelight


text:
module RayCastWorlds

import ReinforcementLearningBase as RLBase

abstract type AbstractGame end

function reset! end
function act! end
function cast_rays! end
function get_action_keys end
function get_action_names end
function play! end
function update_top_view! end
function update_camera_view! end

include("utils.jl")
include("collision_detection.jl")
include("rlbase.jl")
include("single_room.jl")

end


text:
function add_node!(ls::Union{LinkStream,DirectedLinkStream}, n::AbstractString)
    push!(ls.V,n)
end

function add_node!(s::Union{StreamGraph,DirectedStreamGraph}, n::Node)
    if n.name ∉ s
        push!(s.V,n.name)
        s.W[n.name] = n
    else
        merge!(s.W[n.name],n)
    end
end

function add_link!(ls::DirectedLinkStream, l::Link)
    if l.from ∉ ls
        add_node!(ls,l.from)
    end
    if l.to ∉ ls
        add_node!(ls,l.to)
    end
    if !haskey(ls.E,l.from)
        ls.E[l.from] = Dict()
    end
    if !haskey(ls.E[l.from],l.to)
        ls.E[l.from][l.to] = l
    else
        merge!(ls.E[l.from][l.to],l)
    end
end

function add_link!(ls::LinkStream, l::Link)
    if l.from ∉ ls
        add_node!(ls,l.from)
    end
    if l.to ∉ ls
        add_node!(ls,l.to)
    end
    if l.to < l.from
        l.from,l.to=l.to,l.from
    end
    if !haskey(ls.E,l.from)
        ls.E[l.from] = Dict()
    end
    if !haskey(ls.E[l.from],l.to)
        ls.E[l.from][l.to] = l
    else
        merge!(ls.E[l.from][l.to],l)
    end
end

function add_link!(s::DirectedStreamGraph, l::Link)
    if l.from ∉ s
        new_from = Node(l.from, l.presence)
        add_node!(s,new_from)
    elseif l.presence ⊈ s.W[l.from].presence
        new_from = Node(l.from, l.presence)
        merge!(s.W[l.from],new_from)
    end
    if l.to ∉ s
        new_to = Node(l.to, l.presence)
        add_node!(s,new_to)
    elseif l.presence ⊈ s.W[l.to].presence
        new_to = Node(l.to, l.presence)
        merge!(s.W[l.to],new_to)
    end
    if !haskey(s.E,l.from)
        s.E[l.from] = Dict()
    end
    if !haskey(s.E[l.from],l.to)
        s.E[l.from][l.to] = l
    else
        merge!(s.E[l.from][l.to],l)
    end
end

function add_link!(s::StreamGraph, l::Link)
    if l.from ∉ s
        new_from = Node(l.from, l.presence)
        add_node!(s,new_from)
    elseif l.presence ⊈ s.W[l.from].presence
        new_from = Node(l.from, l.presence)
        merge!(s.W[l.from],new_from)
    end
    if l.to ∉ s
        new_to = Node(l.to, l.presence)
        add_node!(s,new_to)
    elseif l.presence ⊈ s.W[l.to].presence
        new_to = Node(l.to, l.presence)
        merge!(s.W[l.to],new_to)
    end
    if l.to < l.from
        l.from,l.to=l.to,l.from
    end
    if !haskey(s.E,l.from)
        s.E[l.from] = Dict()
    end
    if !haskey(s.E[l.from],l.to)
        s.E[l.from][l.to] = l
    else
        merge!(s.E[l.from][l.to],l)
    end
end

function record!(s::AbstractStream, t0::Float64, t1::Float64, from::AbstractString, to::AbstractString)
    if (t0,t1) ⊈ s
        throw("Stream $(s.name) is defined over $(s.T). Invalid link between t0=$t0 and t1=$t1.")
    end
    new = Link("$from to $to", Intervals([(t0,t1)]), from, to, 1)
    add_link!(s,new)
end

function record!(s::Union{StreamGraph,DirectedStreamGraph}, t0::Float64, t1::Float64, n::AbstractString)
    if (t0,t1) ⊈ s
        throw("Stream $(s.name) is defined over $(s.T). Invalid link between t0=$t0 and t1=$t1.")
    end
    new = Node(n,Intervals([(t0,t1)]))
    add_node!(s,new)
end


text:
export Simulate

function Initialize_MD!(sys::System,comm::MPI.Comm,root::Int64)
    CalculateForce!(sys,comm)
    UpdateThermo!(sys,comm,root)
    DumpInformation(sys,comm,root)
    LoggerInformation!(sys,comm,root)
    sys.first_step[1]=false
end

function Initialize_LocalSystem!(sys,
                                local_types,
                                local_masses,
                                local_energies,
                                local_coords,
                                local_forces,
                                local_vels,
                                local_bonds)
        for i in 1:size(local_types)[1]
            push!(sys.atoms.types,local_types[i])
        end

        for i in 1:size(local_masses)[1]
            push!(sys.atoms.masses,local_masses[i])
        end

        for i in 1:size(local_energies)[1]
            push!(sys.atoms.energies,local_energies[i])
        end

        for i in 1:size(local_coords)[1]
            push!(sys.atoms.coords,local_coords[i,1:3])
        end

        for i in 1:size(local_forces)[1]
            push!(sys.atoms.forces,local_forces[i,1:3])
        end

        for i in 1:size(local_vels)[1]
            push!(sys.atoms.vels,local_vels[i,1:3])
        end

        for i in 1:size(local_bonds)[1]
            push!(sys.atoms.bonds,local_bonds[i,1:3])
        end
        if sys.field!=NoField()
            for i in 1:size(local_types)[1]
                push!(sys.field.particle2cell,zeros(8))
                push!(sys.field.cellindex,zeros(Int64,8,3))
            end
        end

end

function Simulate(inputs::Dict,config::Dict)

    MPI.Init()
    comm=MPI.COMM_WORLD
    root=0
    my_rank = MPI.Comm_rank(comm)

    # Data Structure

    if my_rank==root
        sys=System(inputs,config,comm,root)
    else
        sys=System(Dict(),Dict(),comm,root)
    end
    if my_rank==root
        p = Progress(sys.steps; showspeed=true)
    end

    sys = MPI.bcast(sys, root, comm)

    local_types = MPI.Scatterv(sys.atoms.types_all, sys.mpi.commsizes, root, comm)
    local_masses = MPI.Scatterv(sys.atoms.masses_all, sys.mpi.commsizes, root, comm)
    local_energies = MPI.Scatterv(sys.atoms.energies_all, sys.mpi.commsizes, root, comm)

    local_coords = DeSerializeArray(MPI.Scatterv(SerializeArray(sys.atoms.coords_all), sys.mpi.commsizes*3, root, comm))
    local_forces = DeSerializeArray(MPI.Scatterv(SerializeArray(sys.atoms.forces_all), sys.mpi.commsizes*3, root, comm))
    local_vels = DeSerializeArray(MPI.Scatterv(SerializeArray(sys.atoms.vels_all), sys.mpi.commsizes*3, root, comm))

    if sys.bdinter!=BdInter[]
        local_bonds = DeSerializeBondArray(MPI.Scatterv(SerializeArray(sys.atoms.bonds_all), sys.mpi.bond_commsizes*3, root, comm))
    else
        local_bonds=zeros(0,3)
    end

    Initialize_LocalSystem!(sys,
                            local_types,
                            local_masses,
                            local_energies,
                            local_coords,
                            local_forces,
                            local_vels,
                            local_bonds)

    # Main Algorithm
    Initialize_MD!(sys,comm,root)

    for stepi in sys.current_step[1]+1:sys.current_step[1]+sys.steps-1
        sys.current_step[1]=stepi
        UpdatePosition!(sys,sys.integrator)
        UpdateVelocity!(sys,sys.integrator,comm)
        MPI.Barrier(comm)
        ClearForce!(sys)
        CalculateForce!(sys,comm)
        MPI.Barrier(comm)
        UpdateThermo!(sys,comm,root)
        MPI.Barrier(comm)
        UpdateVelocity!(sys,sys.integrator,comm)
        MPI.Barrier(comm)
        Thermostat!(sys,comm,root)
        MPI.Barrier(comm)
        if my_rank==root
            ProgressMeter.next!(p)
        end
        UpdateThermo!(sys,comm,root)
        LoggerInformation!(sys,comm,root)
        DumpInformation(sys,comm,root)
        RestartInformation(sys,comm,root)
    end
    MPI.Finalize()
end


text:
@testset "Test space state simulation" begin
    # Declare functions
    function f(x_next, x, k)
        x_next[1] = x[2]
        x_next[2] = x[1]
    end
    function g(y, x, k)
        y[1] = x[1]
    end
    # Initialize vector
    time_span = (1, 5)
    x0 = [1, 0]
    y = Vector{typeof(x0)}(time_span[2] - time_span[1]+1)
    for i = 1:length(y)
        y[i] = [0]
    end
    buffer = deepcopy(x0)
    # Simulate
    x = ms.simulate_space_state!(y, f, g, x0, time_span;
                                 x0_2=buffer)
    # Test
    @test y == [[1], [0], [1], [0], [1]]
    @test x == [0, 1]
    @test x0 == x
end

@testset "Test space state simulation extra args" begin
    # Declare functions
    function f(x_next, x, k, a, b, c)
        x_next[1] = a[1]*b[1]*c[2]*x[2]
        x_next[2] = a[1]*b[1]*c[2]*x[1]
    end
    function g(y, x, k, a, b, c)
        y[1] = a[2]*b[2]*c[1]*x[1]
    end
    # Initialize vector
    time_span = (1, 5)
    x0 = [1, 0]
    y = Vector{typeof(x0)}(time_span[2] - time_span[1]+1)
    for i = 1:length(y)
        y[i] = [0]
    end
    buffer = deepcopy(x0)
    # Simulate
    a = Float64[1, 1]
    b = Float64[1, 1]
    c = Float64[1, 1]
    args = (a, b, c)
    x = ms.simulate_space_state!(y, f, g, x0, time_span, args;
                                 x0_2=buffer)
    # Test
    @test y == [[1], [0], [1], [0], [1]]
    @test x == [0, 1]
    @test x0 == x
end

@testset "Test space state simulation tuple" begin
    # Declare functions
    function f(x_next, x, k, a, b, c)
        x_next1, x_next2 = x_next
        x1, x2 = x
        x_next1[1] = a[1]*b[1]*c[2]*x1[2]
        x_next1[2] = a[1]*b[1]*c[2]*x1[1]
        x_next2[1] = a[1]*b[1]*c[2]*x2[2]
        x_next2[2] = a[1]*b[1]*c[2]*x2[1]
    end
    function g(y, x, k, a, b, c)
        y1, y2 = y
        x1, x2 = x
        y1[1] = a[2]*b[2]*c[1]*x1[1]
        y2[1] = a[2]*b[2]*c[1]*x2[1]
    end
    # Initialize vector
    time_span = (1, 5)
    x0 = ([1, 0], [0, 1])
    y = Vector{typeof(x0)}(time_span[2] - time_span[1]+1)
    for i = 1:length(y)
        y[i] = ([0], [0])
    end
    buffer = deepcopy(x0)
    # Simulate
    a = Float64[1, 1]
    b = Float64[1, 1]
    c = Float64[1, 1]
    args = (a, b, c)
    x = ms.simulate_space_state!(y, f, g, x0, time_span, args;
                                 x0_2=buffer)
    # Test
    @test y == [([1], [0]),
                ([0], [1]),
                ([1], [0]),
                ([0], [1]),
                ([1], [0])]
    @test x == ([0, 1], [1, 0])
    @test x0 == x
end


text:
using DFControl, Test
import DFControl: Exec, data, execs

prevdefault = getdefault_server()
setdefault_server("localhost")
setdefault_pseudodir(:test, joinpath(testdir, "testassets", "pseudos"))
configuredefault_pseudos(pseudo_dirs=Dict(:test => getdefault_pseudodirs()[:test]))

@test DFControl.getdefault_pseudo(:Pt, :test, specifier="UPF") == Pseudo("Pt.UPF",joinpath(testdir, "testassets", "pseudos"))
@test DFControl.getdefault_server() == "localhost"


text:
export @sum

function sum_expr(ex)
    @match ex begin
        alpha_Number * expr_ => :($alpha * $(sum_expr(expr)))
        expr_ * alpha_Number => :($alpha * $(sum_expr(expr)))

        expr1_ + expr2_ => :($(sum_expr(expr1)) + $(sum_expr(expr2)))
        expr1_ - expr2_ => :($(sum_expr(expr1)) - $(sum_expr(expr2)))
        _ => :(sum($(make_subham(SubHamExpr(ex)))))
    end
end

"""
    @sum local_hamiltonian

Sum over a region defined in `local_hamiltonian`.
"""
macro sum(expr)
    sum_expr(expr)
end

abstract type Operator end

Base.eltype(x::Operator) = eltype(x.expr)

"""
    Sum{H} <: Operator

Symbolic summation.

```math
\\sum H
```
"""
struct Sum{H} <: Operator
    expr::H
end

# NOTE: prod only operates on surround/face?
"""
    Prod{H} <: Operator

Symbolic product.

```math
\\prod H
```
"""
struct Prod{H} <: Operator
    expr::H
end

struct ScalarProd{T <: Number, H} <: Operator
    val::T
    expr::H
end

ScalarProd(val::Number, ex::ScalarProd) = ScalarProd(val * ex.val, ex.expr)

abstract type AddOrSub{LHS, RHS} <: Operator end

struct Add{LHS, RHS} <: AddOrSub{LHS, RHS}
    lhs::LHS
    rhs::RHS
end

struct Sub{LHS, RHS} <: AddOrSub{LHS, RHS}
    lhs::LHS
    rhs::RHS
end

Base.eltype(ex::AddOrSub) = promote_type(eltype(ex.lhs), eltype(ex.rhs))

Base.sum(x::SubHam) = Sum(x)
# move scalar out of sum
Base.sum(x::ScalarProd{<:Number}) = ScalarProd(x.val, Sum(x.expr))
# move add/sub out of sum
Base.sum(x::Add) = Add(Sum(x.lhs), Sum(x.rhs))
Base.sum(x::Sub) = Sub(Sum(x.lhs), Sum(x.rhs))

import Base: *, +, -
*(lhs::Number, rhs::Sum) = ScalarProd(lhs, rhs)
*(lhs::Number, rhs::SubHam) = ScalarProd(lhs, rhs)
*(lhs::Number, rhs::ScalarProd) = ScalarProd(lhs * rhs.val, rhs.expr)
*(lhs::Number, rhs::Add) = Add(lhs * rhs.lhs, lhs * rhs.rhs)
*(lhs::Number, rhs::Sub) = Sub(lhs * rhs.lhs, lhs * rhs.rhs)
# scalar multiplication is exchangable
*(lhs::Sum, rhs::Number) = *(rhs, lhs)
*(lhs::SubHam, rhs::Number) = *(rhs, lhs)
*(lhs::AddOrSub, rhs::Number) = *(rhs, lhs)

# minus
-(x::SubHam) = ScalarProd(-1, x)
-(x::Sum) = ScalarProd(-1, x)
-(x::ScalarProd) = ScalarProd(-x.val, x.expr)
-(x::Add) = Sub(-x.lhs, x.rhs)
-(x::Sub) = Add(-x.lhs, x.rhs)

+(lhs::T1, rhs::T2) where {T1 <: Union{Operator, SubHam}, T2 <: Union{Operator, SubHam}} = Add(lhs, rhs)
-(lhs::T1, rhs::T2) where {T1 <: Union{Operator, SubHam}, T2 <: Union{Operator, SubHam}} = Sub(lhs, rhs)
# merge same op
+(lhs::T, rhs::T) where {T <: Union{Operator, SubHam}} = lhs === rhs ? ScalarProd(2, lhs) : Add(lhs, rhs)
-(lhs::T, rhs::T) where {T <: Union{Operator, SubHam}} = lhs === rhs ? ScalarProd(2, lhs) : Sub(lhs, rhs)

# add/sub merge rules

## ScalarProd
+(lhs::ScalarProd{<:Number, T}, rhs::T) where {T <: Union{Operator, SubHam}} = lhs.expr === rhs ? ScalarProd(lhs.val+1, lhs.expr) : Add(lhs, rhs)
-(lhs::ScalarProd{<:Number, T}, rhs::T) where {T <: Union{Operator, SubHam}} = lhs.expr === rhs ? ScalarProd(lhs.val-1, lhs.expr) : Sub(lhs, rhs)
+(lhs::ScalarProd{<:Number, T}, rhs::ScalarProd{<:Number, T}) where T =
    lhs.expr === rhs.expr ? ScalarProd(lhs.val + rhs.val, lhs.expr) : Add(lhs, rhs)
-(lhs::ScalarProd{<:Number, T}, rhs::ScalarProd{<:Number, T}) where T =
    lhs.expr === rhs.expr ? ScalarProd(lhs.val - rhs.val, lhs.expr) : Sub(lhs, rhs)

+(lhs::ScalarProd{<:Number, T}, rhs::Sum{T}) where T =
    lhs.expr === rhs.expr ? ScalarProd(lhs.val + 1, rhs) : Add(lhs, rhs)
-(lhs::ScalarProd{<:Number, T}, rhs::Sum{T}) where T =
    lhs.expr === rhs.expr ? ScalarProd(lhs.val - 1, rhs) : Sub(lhs, rhs)

+(lhs::ScalarProd{<:Number, Sum{T}}, rhs::Add{Sum{T}, <:Operator}) where T =
    lhs.expr === rhs.lhs.expr ? Add(ScalarProd(lhs.val+1, lhs.expr), rhs.rhs) : Add(lhs, rhs)
+(lhs::ScalarProd{<:Number, Sum{T}}, rhs::Add{<:Operator, Sum{T}}) where T =
    lhs.expr === rhs.rhs.expr ? Add(ScalarProd(lhs.val+1, lhs.expr), rhs.lhs) : Add(lhs, rhs)
+(lhs::ScalarProd{<:Number, Sum{T}}, rhs::Sub{Sum{T}, <:Operator}) where T =
    lhs.expr === rhs.lhs.expr ? Sub(ScalarProd(lhs.val+1, lhs.expr), rhs.rhs) : Add(lhs, rhs)
+(lhs::ScalarProd{<:Number, Sum{T}}, rhs::Sub{<:Operator, Sum{T}}) where T =
    lhs.expr === rhs.rhs.expr ? Add(ScalarProd(lhs.val-1, lhs.expr), rhs.lhs) : Add(lhs, rhs)

-(lhs::ScalarProd{<:Number, Sum{T}}, rhs::Add{Sum{T}, <:Operator}) where T =
    lhs.expr === rhs.lhs.expr ? Sub(ScalarProd(lhs.val-1, lhs.expr), rhs.rhs) : Sub(lhs, rhs)
-(lhs::ScalarProd{<:Number, Sum{T}}, rhs::Add{<:Operator, Sum{T}}) where T =
    lhs.expr === rhs.rhs.expr ? Sub(ScalarProd(lhs.val-1, lhs.expr), rhs.lhs) : Sub(lhs, rhs)
-(lhs::ScalarProd{<:Number, Sum{T}}, rhs::Sub{Sum{T}, <:Operator}) where T =
    lhs.expr === rhs.lhs.expr ? Add(ScalarProd(lhs.val-1, lhs.expr), rhs.rhs) : Sub(lhs, rhs)
-(lhs::ScalarProd{<:Number, Sum{T}}, rhs::Sub{<:Operator, Sum{T}}) where T =
    lhs.expr === rhs.rhs.expr ? Sub(ScalarProd(lhs.val+1, lhs.expr), rhs.lhs) : Sub(lhs, rhs)

+(lhs::T, rhs::ScalarProd{<:Number, T}) where T = +(rhs, lhs)
-(lhs::T, rhs::ScalarProd{<:Number, T}) where T = -(rhs, lhs)
+(lhs::Sum, rhs::ScalarProd) = +(rhs, lhs)
+(lhs::Add, rhs::ScalarProd) = +(rhs, lhs)
+(lhs::Sub, rhs::ScalarProd) = +(rhs, lhs)
-(lhs::Add, rhs::ScalarProd) = +(rhs, lhs)
-(lhs::Sub, rhs::ScalarProd) = +(rhs, lhs)

################################# Geometry #####################################

struct OnRegion{E, R} <: Operator
    expr::E
    region::R

    OnRegion(e::E, r::R) where {E, R} = issum(e) ? new{E, R}(e, r) : error("Expect a sum expression, got $e")
end

Base.eltype(x::OnRegion) = eltype(x.expr)

issum(ex::Add) = issum(ex.lhs) && issum(ex.rhs)
issum(ex::Sub) = issum(ex.lhs) && issum(ex.rhs)
issum(ex::Sum) = true
issum(ex::ScalarProd) = issum(ex.expr)
issum(ex::Prod) = false
issum(ex::SubHam) = false

# NOTE: this is a workaround
(ex::Add)(region) = OnRegion(ex, region)
(ex::Sub)(region) = OnRegion(ex, region)
(ex::Sum)(region) = OnRegion(ex, region)
(ex::ScalarProd)(region) = OnRegion(ex, region)

################################# Printing #####################################
Base.show(io::IO, x::Sum{<:SubHam}) = print(io, "sum(", x.expr, ")")
Base.show(io::IO, x::Prod{<:SubHam}) = print(io, "prod(", x.expr, ")")
Base.show(io::IO, x::ScalarProd) = print(io, x.val, "*", x.expr)
Base.show(io::IO, x::Add) = print(io, x.lhs, " + ", x.rhs)
Base.show(io::IO, x::Sub) = print(io, x.lhs, " - ", x.rhs)

function Base.show(io::IO, x::OnRegion)
    printstyled(io, "Expr:\n"; bold=true)
    println(io, x.expr)
    printstyled(io, "Region:\n"; bold=true)
    print(io, x.region)
end

using LaTeXStrings
import LaTeXStrings: latexstring

function _latexstring(x::Sum{<:SubHam})
    index_str = join(x.expr.placeholders, ", ")
    "\\sum_{$(index_str)} $(_latexstring(x.expr))"
end

function _latexstring(x::Prod{<:SubHam})
    index_str = join(x.expr.placeholders, ", ")
    "\\prod_{$(index_str)} $(_latexstring(x.expr))"
end

_latexstring(x::ScalarProd) = "$(x.val) $(_latexstring(x.expr))"
_latexstring(x::Add) = "$(_latexstring(x.lhs) + _latexstring(x.rhs))"
_latexstring(x::Sub) = "$(_latexstring(x.lhs) - _latexstring(x.rhs))"
latexstring(x::Operator) = latexstring("\$\$", _latexstring(x), "\$\$")

Base.show(io::IO, ::MIME"text/latex", x::Operator) = print(io, latexstring(x))

function Base.show(io::IO, ::MIME"text/latex", x::OnRegion)
    printstyled(io, "Expr:\n"; bold=true)
    println(io, latexstring(x.expr))
    printstyled(io, "Region:\n"; bold=true)
    print(io, x.region)
end

max_stars_repo_name: Roger-luo/QuHamiltonian

max_issues_repo_name: Roger-luo/QuHamiltonian

max_forks_repo_name: Roger-luo/QuHamiltonian


text:
@testset "Reddit_full" begin
    data  = Reddit(full=true)
    @test length(data) == 1
    g = data[1]
    @test g.num_nodes == 232965
    @test g.num_edges == 114615892
    @test size(g.node_data.features) == (602, g.num_nodes)
    @test size(g.node_data.labels) == (g.num_nodes,)
    @test size(data.metadata["split"]["train"]) == (153431,)
    @test size(data.metadata["split"]["val"]) == (23831,)
    @test size(data.metadata["split"]["test"]) == (55703,)
    s, t = g.edge_index
    @test length(s) == length(t) == g.num_edges
    @test minimum(s) == minimum(t) == 1
    @test maximum(s) == maximum(t) == g.num_nodes
    @test sum(length.(values(data.metadata["split"]))) == g.num_nodes
end

@testset "Reddit_subset" begin
    data  = Reddit(full=false)
    @test length(data) == 1
    g = data[1]
    @test g.num_nodes == 231443
    @test g.num_edges == 23213838
    @test size(g.node_data.features) == (602, g.num_nodes)
    @test size(g.node_data.labels) == (g.num_nodes,)
    @test size(data.metadata["split"]["train"]) == (152410,)
    @test size(data.metadata["split"]["val"]) == (23699,)
    @test size(data.metadata["split"]["test"]) == (55334,)
    s, t = g.edge_index
    @test length(s) == length(t) == g.num_edges
    @test minimum(s) == minimum(t) == 1
    @test maximum(s) == maximum(t) == g.num_nodes
    @test sum(length.(values(data.metadata["split"]))) == g.num_nodes
end

@testset "TUDataset - PROTEINS" begin
    data  = TUDataset("PROTEINS")

    @test data.num_nodes == 43471
    @test data.num_edges == 162088
    @test data.num_graphs == 1113

    @test data.num_nodes == sum(g->g.num_nodes, data.graphs)
    @test data.num_edges == sum(g->g.num_edges, data.graphs)
    @test data.num_edges == sum(g->length(g.edge_index[1]), data.graphs)
    @test data.num_edges == sum(g->length(g.edge_index[2]), data.graphs)
    @test data.num_graphs == length(data) == length(data.graphs)

    i = rand(1:length(data))
    di = data[i]
    @test di isa NamedTuple
    g, targets = di.graphs, di.targets
    @test targets == 1
    @test g isa Graph
    @test all(1 .<= g.edge_index[1] .<= g.num_nodes)
    @test all(1 .<= g.edge_index[2] .<= g.num_nodes)

    # graph data
    @test size(data.graph_data.targets) == (data.num_graphs,)

    # node data
    @test size(g.node_data.features) == (1, g.num_nodes)
    @test size(g.node_data.targets) == (g.num_nodes,)

    # edge data
    @test g.edge_data === nothing
end

@testset "TUDataset - QM9" begin
    data  = TUDataset("QM9")

    @test data.num_nodes == 2333625
    @test data.num_edges == 4823498
    @test data.num_graphs == 129433

    @test data.num_nodes == sum(g->g.num_nodes, data.graphs)
    @test data.num_edges == sum(g->g.num_edges, data.graphs)
    @test data.num_edges == sum(g->length(g.edge_index[1]), data.graphs)
    @test data.num_edges == sum(g->length(g.edge_index[2]), data.graphs)
    @test data.num_graphs == length(data) == length(data.graphs)

    i = rand(1:length(data))
    g, features = data[i]
    @test g isa Graph
    @test all(1 .<= g.edge_index[1] .<= g.num_nodes)
    @test all(1 .<= g.edge_index[2] .<= g.num_nodes)

    # graph data
    @test size(data.graph_data.features) == (19, data.num_graphs)

    # node data
    @test size(g.node_data.features) == (16, g.num_nodes)

    # edge data
    @test size(g.edge_data.features) == (4, g.num_edges)
end


text:
# This is part of Calendars.jl. See the copyright note there.
# ========================= Hebrew dates ====================

# Day number of the start of Hebrew calendar.
# The Jewish calendar's epoch, 1 Tishri AM 1,
# is equivalent to Monday, 7 October 3761 BCE,
# in the proleptic Julian calendar.
# AM-0001-01-01  = JD--3761-10-07 ~ DN#-1373427 = EpochHebrew
# AM-3761-10-16  = JD-0001-01-01  ~ DN#-1
# AM-4872-09-09  = JD-1111-11-11  ~ DN#0405740
# AM-13760-11-13 = JD-9999-12-31  ~ DN#33652132
const EpochHebrew = -1373427
const ValidYearsHebrew = (3761, 13759) # by convention

# True if Heshvan is long in Hebrew year.
function isLongHeshvan(year::DPart)
    rem(DaysInYearHebrew(year), 10) == 5
end

# True if Kislev is short in Hebrew year.
function isShortKislev(year::DPart)
    rem(DaysInYearHebrew(year), 10) == 3
end

# True if year is an Hebrew leap year
function isLeapYearHebrew(year::DPart)
    rem(7 * year + 1, 19) < 7
end

# True if month is long (30 days).
function isLongMonth(year::DPart, month::DPart)
    (  month == 2   || month == 4 || month == 6
    || month == 10  || month == 13
    || (month == 8  && ! isLongHeshvan(year))
    || (month == 9  && isShortKislev(year))
    || (month == 12 && ! isLeapYearHebrew(year))
    )
end

# Last day of month in Hebrew year.
function LastDayOfMonthHebrew(year::DPart, month::DPart)
    isLongMonth(year, month) ? 29 : 30
end

# Last month of Hebrew year.
function LastMonthOfYearHebrew(year::DPart)
    isLeapYearHebrew(year) ? 13 : 12
end

# Return the number of days in the given Hebrew year.
function DaysInYearHebrew(year::DPart)
    DelayHebrew(year + 1) - DelayHebrew(year)
end

# Is the date a valid Hebrew date?
function isValidDateHebrew(cd::CDate, warn=true)
    cal, year, month, day = cd
    val = ( cal == AM
        && (ValidYearsHebrew[1] <= year && year <= ValidYearsHebrew[2])
        && (1 <= month && month <= LastMonthOfYearHebrew(year))
        && (1 <= day && day <= LastDayOfMonthHebrew(year, month)))
    ! val && warn && @warn(Warning(cd))
    return val
end

# Represents the day number a valid Hebrew date?
function isValidDateHebrew(dn::DPart, warn=true)
    val = EpochHebrew <= dn
    ! val && warn && @warn(Warning(dn))
    return val
end

# Number of days elapsed from the Sunday prior to the start of the
# Hebrew calendar to the mean DelayHebrew of Tishri of Hebrew year.
function DelayHebrew(year::DPart)

    MonthsElapsed = ( 235 * div(year - 1, 19)  # Months in complete cycles so far.
                     + 12 * rem(year - 1, 19)  # Regular months in this cycle.
                     + div(7 * rem(year - 1, 19) + 1, 19)) # Leap months this cycle

    PartsElapsed = 204 + 793 * rem(MonthsElapsed, 1080)
    HoursElapsed = (5 + 12 * MonthsElapsed
                 + 793 * div(MonthsElapsed, 1080) + div(PartsElapsed, 1080))
    DelayHebrewDay = 1 + 29 * MonthsElapsed + div(HoursElapsed, 24)
    DelayHebrewParts = 1080 * rem(HoursElapsed, 24) + rem(PartsElapsed, 1080)
    DayOfWeek = rem(DelayHebrewDay, 7)

    if (
        DelayHebrewParts >= 19440           # If new moon is at or after midday,
        || (DayOfWeek == 2                  # ...or is on a Tuesday...
            && DelayHebrewParts >= 9924     # at 9 hours, 204 parts or later...
            && ! isLeapYearHebrew(year))    # ...of a common year,
        || (DayOfWeek == 1                  # ...or is on a Monday at...
            && DelayHebrewParts >= 16789    # 15 hours, 589 parts or later...
            && isLeapYearHebrew(year - 1))  # at the end of a leap year
    )
        # Then postpone Rosh HaShanah one day
        DelayHebrewDay += 1
        DayOfWeek = rem(DelayHebrewDay, 7)
    end

    if ( DayOfWeek == 0  # If Rosh HaShanah would occur on Sunday,
      || DayOfWeek == 3  # or Wednesday,
      || DayOfWeek == 5  # or Friday
    )
      # Then postpone it one (more) day
      DelayHebrewDay += 1
    end

    return DelayHebrewDay
end

# Return the days this year so far.
function DayOfYearHebrew(year::DPart, month::DPart, day::DPart)
    DaysInYear = day

    if month < 7   # Before Tishri, so add days in prior months
        # this year before and after Nisan.
        for m in 7:LastMonthOfYearHebrew(year)
            DaysInYear += LastDayOfMonthHebrew(year, m)
        end
        for m in 1:month - 1
            DaysInYear += LastDayOfMonthHebrew(year, m)
        end
    else # Add days in prior months this year
        for m in 7:month - 1
            DaysInYear += LastDayOfMonthHebrew(year, m)
        end
    end
    return DaysInYear - 1
end

# Computes the day number from a valid Hebrew date.
function DayNumberHebrew(year::DPart, month::DPart, day::DPart)
    ( EpochHebrew - 1
    + DayOfYearHebrew(year, month, day)
    + DelayHebrew(year)
    )
end

# Computes the day number of a Hebrew date.
DayNumberHebrew(d::CDate) = DayNumberHebrew(d[2], d[3], d[4])

# Computes the Hebrew date from an DayNumber.
function DateHebrew(dn::DPart)
    ! isValidDateHebrew(dn) && return InvalidDate

    year = div(dn + EpochHebrew, 366) # Approximation from below.
    # Search forward for year from the approximation.
    while dn >= DayNumberHebrew(year + 1, 7, 1)
        year += 1
    end

    # Search forward for month from either Tishri (7) or Nisan (1).
    month = dn < DayNumberHebrew(year, 1, 1) ? 7 : 1

    while dn > DayNumberHebrew(year, month, (LastDayOfMonthHebrew(year, month)))
        month += 1
    end

    # Calculate the day by subtraction.
    day = dn - DayNumberHebrew(year, month, 1) + 1

    return (AM, year, month, day)::CDate
end

# Return the day number of the first day in the given Hebrew year.
function YearStartHebrew(year::DPart)
    return DayNumberHebrew(year, 7, 1)
end

# Return the day number of the last day in the given Hebrew year.
function YearEndHebrew(year::DPart)
    return YearStartHebrew(year + 1) - 1
end

# Return the number of months in the given Hebrew year.
function MonthsInYearHebrew(year::DPart)
    return LastMonthOfYearHebrew(year)
end


text:
# ---
# title: 566. Reshape the Matrix
# id: problem566
# author: Qling
# date: 2021-02-17
# difficulty: Easy
# categories: Array
# link: <https://leetcode.com/problems/reshape-the-matrix/description/>
# hidden: true
# ---
#
# In MATLAB, there is a very useful function called 'reshape', which can reshape
# a matrix into a new one with different size but keep its original data.
#
# You're given a matrix represented by a two-dimensional array, and two
# **positive** integers **r** and **c** representing the **row** number and
# **column** number of the wanted reshaped matrix, respectively.
#
# The reshaped matrix need to be filled with all the elements of the original
# matrix in the same **row-traversing** order as they were.
#
# If the 'reshape' operation with given parameters is possible and legal, output
# the new reshaped matrix; Otherwise, output the original matrix.
#
# **Example 1:**
#
#
#
#     Input:
#     nums =
#     [[1,2],
#      [3,4]]
#     r = 1, c = 4
#     Output:
#     [[1,2,3,4]]
#     Explanation:
#      The **row-traversing** of nums is [1,2,3,4]. The new reshaped matrix is a 1 * 4 matrix, fill it row by row by using the previous list.
#
#
# **Example 2:**
#
#
#
#     Input:
#     nums =
#     [[1,2],
#      [3,4]]
#     r = 2, c = 4
#     Output:
#     [[1,2],
#      [3,4]]
#     Explanation:
#      There is no way to reshape a 2 * 2 matrix to a 2 * 4 matrix. So output the original matrix.
#
#
# **Note:**
#
#   1. The height and width of the given matrix is in range [1, 100].
#   2. The given r and c are all positive.
#
#
## @lc code=start
using LeetCode

function matrixReshape(nums::Vector{Vector{Int}}, r::Int, c::Int)
    m, n = length(nums), length(nums[1])

    m * n != r * c && return nums

    ans = [[0 for _ in 1:c] for _ in 1:r]

    for x in 0:(m * n - 1)
        ans[x ÷ c + 1][(x % c) + 1] = nums[x ÷ n + 1][(x % n) + 1]
    end

    return ans
end

## @lc code=end


text:
using ProtoBuf
import ProtoBuf.meta

mutable struct V1 <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function V1(; kwargs...)
        obj = new(meta(V1), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end
const __meta_V1 = Ref{ProtoMeta}()
function meta(::Type{V1})
    if !isassigned(__meta_V1)
        __meta_V1[] = target = ProtoMeta(V1)
        allflds = Pair{Symbol,Union{Type,String}}[:f1 => Int32, :f2 => Bool]
        meta(target, V1, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
    end
    __meta_V1[]
end
function Base.getproperty(obj::V1, name::Symbol)
    if name === :f1
        return (obj.__protobuf_jl_internal_values[name])::Int32
    elseif name === :f2
        return (obj.__protobuf_jl_internal_values[name])::Bool
    else
        getfield(obj, name)
    end
end

mutable struct V2 <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function V2(; kwargs...)
        obj = new(meta(V2), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end
const __meta_V2 = Ref{ProtoMeta}()
function meta(::Type{V2})
    if !isassigned(__meta_V2)
        __meta_V2[] = target = ProtoMeta(V2)
        allflds = Pair{Symbol,Union{Type,String}}[:f1 => Int64, :f2 => Bool, :f3 => Int64]
        meta(target, V2, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
    end
    __meta_V2[]
end
function Base.getproperty(obj::V2, name::Symbol)
    if name === :f1
        return (obj.__protobuf_jl_internal_values[name])::Int64
    elseif name === :f2
        return (obj.__protobuf_jl_internal_values[name])::Bool
    elseif name === :f3
        return (obj.__protobuf_jl_internal_values[name])::Int64
    else
        getfield(obj, name)
    end
end

function check_samestruct()
    iob = PipeBuffer();
    writeval = V1(; f1=1, f2=true);
    writeproto(iob, writeval);
    readval = readproto(iob, V1());
    @test readval == writeval

    iob = PipeBuffer();
    writeval = V2(; f1=1, f2=true, f3=20);
    writeproto(iob, writeval);
    readval = readproto(iob, V2());
    @test readval == writeval
end

# write V1, read V2
# write V2, read V1
function check_V1_v2()
    iob = PipeBuffer();
    writeval = V1(; f1=1, f2=true);
    writeproto(iob, writeval);
    readval = readproto(iob, V2());
    checkval = V2(; f1=1, f2=true, f3=0)
    @test readval == checkval

    iob = PipeBuffer();
    writeval = V2(; f1=1, f2=true, f3=20);
    writeproto(iob, writeval);
    readval = readproto(iob, V1());
    checkval = V1(; f1=1, f2=true)
    @test readval == checkval

    iob = PipeBuffer();
    writeval = V2(; f1=typemax(Int64), f2=true, f3=20);
    writeproto(iob, writeval);
    readval = readproto(iob, V1());
    checkval = V1(; f1=-1, f2=true)  # overflow can't be detected as negative int32 is sign extended for serialization
    @test readval == checkval
end

@testset "Serialization across type versions" begin
    check_samestruct()
    check_V1_v2()
end


text:
module ThAD

using MacroTools: @forward
using ..ThArrays: Tensor, Scalar, TorchNumber
using ..ThC

@static if Sys.islinux()
    using LibTorchCAPI_jll
elseif Sys.isapple()
    const libtorch_capi = :libtorch_capi
end

import ..ThC: grad, requires_grad!

function has_grad(a::Tensor)
    ret = ccall((:tensor_method_has_grad, libtorch_capi),
                Cint, (Ptr{Cvoid},), a.pointer)
    return ret != 0
end

function get_grad(a::Tensor, default=nothing)
    if has_grad(a)
        return grad(a)
    else
        default == nothing ? ThC.zeros_like(a) : default
    end
end

function backward(a::Tensor, d::Union{Ptr{Nothing}, Tensor}=C_NULL;
                  keep_graph::Bool=false, create_graph::Bool=false)
    if d isa Tensor
        d = d.pointer
    end
    ccall((:tensor_method_backward, libtorch_capi),
          Ptr{Cvoid}, (Ptr{Cvoid}, Ptr{Cvoid}, Cint, Cint),
          a.pointer, d, keep_graph, create_graph)
    nothing
end

reset_grad!(t::Tensor) = ThC.zero!(ThC.grad(t))
requires_grad!(t::Tensor, r::Bool) = requires_grad!(t, r ? 1 : 0)

function gradient(f, data...; d::Union{Ptr{Nothing}, Tensor}=C_NULL)
    tensors = map(d -> Tensor(d, requires_grad=true), data)
    return gradient(f, tensors...; d=d)
end

function gradient(f, tensors::Vararg{Tensor}; d::Union{Ptr{Nothing}, Tensor}=C_NULL)
    result = f(tensors...)
    backward(result, d)
    return ThC.grad.(tensors)
end

# tracker compatible API

struct Params
    order::Vector{Any}
    params::IdDict{Any, Bool}
    Params() = new([], IdDict())
end

@forward Params.order Base.iterate, Base.length

function Base.push!(ps::Params, x)
  if !haskey(ps.params, x)
      push!(ps.order, x)
      ps.params[x] = true
  end
  return ps
end

Base.push!(ps::Params, x...) = (foreach(x -> push!(ps, x), x); ps)
Params(xs) = push!(Params(), xs...)

data(t::Tensor) = convert(Array, t)
data(t::Tensor{T, 0}) where T = t[]
param(x) = x
param(x::Number) = Tensor(float(x); requires_grad=true)
param(xs::AbstractArray) = Tensor(float.(xs); requires_grad=true)

function forward(f, ps::Params)
    y = f()
    back = (d) -> begin
        g = IdDict()
        # reset grad!
        foreach((t) -> has_grad(t) && reset_grad!(t), ps)
        backward(y, param(d))
        foreach((t) -> g[t] = ThC.grad(t), ps)
        return g
    end
    return data(y), back
end

function forward(f, args...)
    args = param.(args)
    y, back = forward(()->f(args...), Params(args))
    y, (d) -> getindex.(Ref(back(d)), args)
end

end


text:
function main()
    arr = collect.(readlines())

    sy, sx = length(arr), length(arr[1])
    is_valid(y, x) = 1 <= y <= sy && 1 <= x <= sx

    flashes = 0
    function flash(y, x)
        if arr[y][x] > '9'
            flashes += 1
            arr[y][x] = '0'
            for (y2, x2) in [(y+a, x+b) for a=-1:1, b=-1:1]
                if is_valid(y2, x2) && arr[y2][x2] != '0'
                    arr[y2][x2] += 1
                    flash(y2, x2)
                end
            end
        end
    end

    total_flashes = 0
    for step in 1:1000
        flashes = 0
        (a -> a .+= 1).(arr)
        for y in 1:sy, x in 1:sx
            flash(y, x)
        end
        total_flashes += flashes
        step == 100 && println(total_flashes)
        flashes == sy * sx && return println(step)
    end
end

main()

max_stars_repo_name: LiquidFun/adventofcode

max_issues_repo_name: LiquidFun/adventofcode

max_forks_repo_name: LiquidFun/adventofcode


text:
fun(x) =  (1.0 - x[1])^2 + 100.0 * (x[2] - x[1]^2)^2

function fun_grad!(g, x)
g[1] = -2.0 * (1.0 - x[1]) - 400.0
* (x[2] - x[1]^2) * x[1]
g[2] = 200.0 * (x[2] - x[1]^2)
end

function fun_hess!(h, x)
h[1, 1] = 2.0 - 400.0 * x[2] + 1200.0 * x[1]^2
h[1, 2] = -400.0 * x[1]
h[2, 1] = -400.0 * x[1]
h[2, 2] = 200.0
end;

x0 = [0.0, 0.0]
df = TwiceDifferentiable(fun, fun_grad!, fun_hess!, x0)

lx = [-0.5, -0.5]; ux = [0.5, 0.5]
dfc = TwiceDifferentiableConstraints(lx, ux)

res = optimize(df, dfc, x0, IPNewton())

ux = fill(Inf, 2)
dfc = TwiceDifferentiableConstraints(lx, ux)

clear!(df)
res = optimize(df, dfc, x0, IPNewton())

lx = fill(-Inf, 2); ux = fill(Inf, 2)
dfc = TwiceDifferentiableConstraints(lx, ux)

clear!(df)
res = optimize(df, dfc, x0, IPNewton())

lx = Float64[]; ux = Float64[]
dfc = TwiceDifferentiableConstraints(lx, ux)

clear!(df)
res = optimize(df, dfc, x0, IPNewton())

con_c!(c, x) = (c[1] = x[1]^2 + x[2]^2; c)
function con_jacobian!(J, x)
    J[1,1] = 2*x[1]
    J[1,2] = 2*x[2]
    J
end
function con_h!(h, x, λ)
    h[1,1] += λ[1]*2
    h[2,2] += λ[1]*2
end;

lx = Float64[]; ux = Float64[]
lc = [-Inf]; uc = [0.5^2]
dfc = TwiceDifferentiableConstraints(con_c!, con_jacobian!, con_h!,
                                     lx, ux, lc, uc)
res = optimize(df, dfc, x0, IPNewton())

lc = [0.1^2]
dfc = TwiceDifferentiableConstraints(con_c!, con_jacobian!, con_h!,
                                     lx, ux, lc, uc)
res = optimize(df, dfc, x0, IPNewton())

function con2_c!(c, x)
    c[1] = x[1]^2 + x[2]^2     ## First constraint
    c[2] = x[2]*sin(x[1])-x[1] ## Second constraint
    c
end
function con2_jacobian!(J, x)
    # First constraint
    J[1,1] = 2*x[1]
    J[1,2] = 2*x[2]
    # Second constraint
    J[2,1] = x[2]*cos(x[1])-1.0
    J[2,2] = sin(x[1])
    J
end
function con2_h!(h, x, λ)
    # First constraint
    h[1,1] += λ[1]*2
    h[2,2] += λ[1]*2
    # Second constraint
    h[1,1] += λ[2]*x[2]*-sin(x[1])
    h[1,2] += λ[2]*cos(x[1])
    # Symmetrize h
    h[2,1]  = h[1,2]
    h
end;

x0 = [0.25, 0.25]
lc = [-Inf, 0.0]; uc = [0.5^2, 0.0]
dfc = TwiceDifferentiableConstraints(con2_c!, con2_jacobian!, con2_h!,
                                     lx, ux, lc, uc)
res = optimize(df, dfc, x0, IPNewton())

# This file was generated using Literate.jl, https://github.com/fredrikekre/Literate.jl


text:
export Counter
export add!

Counter = Dict{String, Int}

function add!(c::Counter, s::String)
    if s ∈ keys(c)
        c[s] += 1
    else
        c[s]  = 1
    end
end

max_stars_repo_name: neherlab/pangraph

max_issues_repo_name: neherlab/pangraph

max_forks_repo_name: neherlab/pangraph


text:
using
  FourierFlows.KuramotoSivashinsky,
  PyPlot

nx = 512
Lx = 32π
dt = 0.1
nt = 2000
prob = Problem(nx=nx, Lx=Lx, dt=dt, stepper="ETDRK4")

x = prob.grid.x
u0 = @. cos(x/16) * (1 + sin(x/16))
set_u!(prob, u0)

T = [ dt*j for i=1:nx, j=1:nt ]
X = [ x[i] for i=1:nx, j=1:nt ]
u = zeros(nx, nt)

for i = 1:nt
  stepforward!(prob)
  updatevars!(prob)
  u[:, i] .= prob.vars.u
end

close("all")
fig, ax = subplots(figsize=(6, 8))
pcolormesh(X/2π, T, u)
ax[:tick_params](bottom=false, left=false, labelbottom=false, labelleft=false)
show()


text:
struct LocalHamiltonian{O<:SumOfLocalTerms, S<:Real}
    h::O
    domain::Tuple{S,S}
end

∫(o::SumOfLocalTerms, domain::Tuple{<:Real,<:Real}) =
    LocalHamiltonian(o, promote(domain...))

∫(o::LocalOperator, domain::Tuple{<:Real,<:Real}) = ∫(1*o, domain)

domain(H::LocalHamiltonian) = H.domain

density(H::LocalHamiltonian) = H.h

Base.:*(α::Number, H::LocalHamiltonian) = LocalHamiltonian(α * H.h, H.domain)
Base.:*(H::LocalHamiltonian, α::Number) = LocalHamiltonian(H.h * α, H.domain)
Base.:\(α::Number, H::LocalHamiltonian) = LocalHamiltonian(α \ H.h, H.domain)
Base.:/(H::LocalHamiltonian, α::Number) = LocalHamiltonian(H.h / α, H.domain)
Base.:+(H::LocalHamiltonian) = LocalHamiltonian(+H.h, H.domain)
Base.:-(H::LocalHamiltonian) = LocalHamiltonian(-H.h, H.domain)

function Base.:+(H1::LocalHamiltonian, H2::LocalHamiltonian)
    H1.domain == H2.domain || error("non-matching domains")
    return LocalHamiltonian(H1.h + H2.h, H1.domain)
end
function Base.:-(H1::LocalHamiltonian, H2::LocalHamiltonian)
    H1.domain == H2.domain || error("non-matching domains")
    return LocalHamiltonian(H1.h - H2.h, H1.domain)
end


text:
using Base.Test
using EGR

println("TestSAG")

createOracleOutputLevel = 1
numEquivalentPasses = 5
algOutputLevel = 0
maxOutputNum=20
numChunks=3

sLikeGd(k,numTP) = k==0 ? 0 : numTP
uLikeGd(k,numTP) = k==0 ? numTP : 0

myREfunction(problem, opts, sd, expIndices,n) = false
myWriteFunction(problem, sd, opts, k, gnum, origWant, fromOutputsFunction) = false

for (gradientOracle, numVars, numTP, csDataType, LossFunctionString, myOutputter, L2reg, thisDataName, thisProblem) in Oracles

	println(" $thisDataName $LossFunctionString L2reg = $L2reg")

	println(" numTP = $numTP")

	maxG  = Int(round(numEquivalentPasses*numTP))

	expIndices=unique(round(Int, logspace(0,log10(maxG+1),maxOutputNum)))-1

	myOutputOpts =  OutputOpts(myOutputter,expIndices)
	myOptss = Opts(zeros(numVars,1); stepSizePower=0, maxG=maxG, outputLevel=algOutputLevel, stepOutputLevel=0)

	problem = thisProblem(Task(() -> getSequentialFinite(numTP, gradientOracle)))
	y=csDataType[]
	for i in 1:numChunks
		println("  Creating chunk $i")
		for j in (i-1)*round(Int,floor(numTP/numChunks))+1:i*round(Int,floor(numTP/numChunks))
			(f,sampleG) = (problem.getSampleFunctionAt(j))(myOptss.init)
			push!(y,sampleG)
		end
	end

	alg(problem,  myOptss, SAG(numVars, copy(y), numTP,numChunks,1), myOutputOpts, myWriteFunction)

	alg(problem,  myOptss, SAG(numVars, copy(y), numTP,numChunks,2), myOutputOpts, myWriteFunction)

end

println("TestSAG successful!")
println()


text:
const use_return  = haskey(ENV, "USE_RETURN" ) ? parse(Bool, ENV["USE_RETURN"] ) : false
const USE_GPU     = haskey(ENV, "USE_GPU"    ) ? parse(Bool, ENV["USE_GPU"]    ) : false
const do_viz      = haskey(ENV, "DO_VIZ"     ) ? parse(Bool, ENV["DO_VIZ"]     ) : true
const do_save     = haskey(ENV, "DO_SAVE"    ) ? parse(Bool, ENV["DO_SAVE"]    ) : false
const do_save_viz = haskey(ENV, "DO_SAVE_VIZ") ? parse(Bool, ENV["DO_SAVE_VIZ"]) : false
const nx          = haskey(ENV, "NX"         ) ? parse(Int , ENV["NX"]         ) : 512
const ny          = haskey(ENV, "NY"         ) ? parse(Int , ENV["NY"]         ) : 512
###
using ParallelStencil
using ParallelStencil.FiniteDifferences2D
@static if USE_GPU
    @init_parallel_stencil(CUDA, Float64, 2)
else
    @init_parallel_stencil(Threads, Float64, 2)
end
using Plots, Printf, LinearAlgebra, MAT

@parallel function compute_flux!(qHx, qHy, qHx2, qHy2, H, D, θr_dτ, dx, dy)
    @all(qHx)  = (@all(qHx) * θr_dτ - D * @d_xi(H) / dx) / (1.0 + θr_dτ)
    @all(qHy)  = (@all(qHy) * θr_dτ - D * @d_yi(H) / dy) / (1.0 + θr_dτ)
    @all(qHx2) = -D * @d_xi(H) / dx
    @all(qHy2) = -D * @d_yi(H) / dy
    return
end

@parallel function compute_update!(H, Hold, qHx, qHy, dτ_ρ, dt, dx, dy)
    @inn(H) = (@inn(H) +  dτ_ρ * (@inn(Hold) / dt - (@d_xa(qHx) / dx + @d_ya(qHy) / dy))) / (1.0 + dτ_ρ / dt)
    return
end

@parallel function check_res!(ResH, H, Hold, qHx2, qHy2, dt, dx, dy)
    @all(ResH) = -(@inn(H) - @inn(Hold)) / dt - (@d_xa(qHx2) / dx + @d_ya(qHy2) / dy)
    return
end

@views function diffusion_2D_()
    # Physics
    lx, ly  = 10.0, 10.0    # domain size
    D       = 1.0           # diffusion coefficient
    ttot    = 1.0           # total simulation time
    dt      = 0.2           # physical time step
    # Numerics
    # nx, ny  = 2*256, 2*256  # numerical grid resolution
    tol     = 1e-8          # tolerance
    itMax   = 1e5           # max number of iterations
    nout    = 10            # tol check
    CFL     = 1/sqrt(2)     # CFL number
    # Derived numerics
    dx, dy  = lx / nx, ly / ny  # grid size
    Vpdτ    = CFL * min(dx, dy)
    Re      = π + sqrt(π^2 + (max(lx, ly)^2 / D / dt)) # Numerical Reynolds number
    θr_dτ   = max(lx, ly) / Vpdτ / Re
    dτ_ρ    = Vpdτ * max(lx, ly) / D / Re
    xc, yc  = LinRange(dx/2, lx - dx/2, nx), LinRange(dy/2, ly - dy/2, ny)
    # Array allocation
    qHx     = @zeros(nx-1, ny-2)
    qHy     = @zeros(nx-2, ny-1)
    qHx2    = @zeros(nx-1, ny-2)
    qHy2    = @zeros(nx-2, ny-1)
    ResH    = @zeros(nx-2, ny-2)
    # Initial condition
    H0      = Data.Array(exp.(-(xc .- lx/2).^2 .- ((yc .- ly/2)').^2))
    Hold    = @ones(nx,ny) .* H0
    H       = @ones(nx,ny) .* H0
    t = 0.0; it = 0; ittot = 0; nt = Int(ceil(ttot/dt))
    # Physical time loop
    while it < nt
        iter = 0; err = 2 * tol
        # Pseudo-transient iteration
        while err > tol && iter < itMax
            @parallel compute_flux!(qHx, qHy, qHx2, qHy2, H, D, θr_dτ, dx, dy)
            @parallel compute_update!(H, Hold, qHx, qHy, dτ_ρ, dt, dx, dy)
            iter += 1
            if iter % nout == 0
                @parallel check_res!(ResH, H, Hold, qHx2, qHy2, dt, dx, dy)
                err = norm(ResH) / sqrt(length(ResH))
            end
        end
        ittot += iter; it += 1; t += dt
        Hold .= H
        if isnan(err) error("NaN") end
    end
    @printf("Total time = %1.2f, time steps = %d, nx = %d, iterations tot = %d \n", round(ttot, sigdigits=2), it, nx, ittot)
    # Visualise
    if do_viz display(heatmap(xc, yc, Array(H'), aspect_ratio=1, framestyle=:box, xlims=(xc[1], xc[end]), ylims=(yc[1], yc[end]), xlabel="lx", ylabel="ly", c=:viridis, clims=(0, 1), title="linear diffusion (nt=$it, iters=$ittot)")) end
    if do_save
        !ispath("../output") && mkdir("../output")
        open("../output/out_diff_2D_lin.txt","a") do io
            println(io, "$(nx) $(ny) $(ittot) $(nt)")
        end
    end
    if do_save_viz
        !ispath("../out_visu") && mkdir("../out_visu")
        matwrite("../out_visu/diff_2D_lin.mat", Dict("H_2D"=> Array(H), "xc_2D"=> Array(xc), "yc_2D"=> Array(yc)); compress = true)
    end
    return xc, yc, H
end

if use_return
    xc, yc, H = diffusion_2D_();
else
    diffusion_2D = begin diffusion_2D_(); return; end
end


text:
include("UpdateExtrema.jl")
function tank_triggers(wn::PyCall.PyObject, wn_dict::Dict{Any,Any}, controls_dict::Dict{Any,Any}, control_names::Base.KeyIterator{Dict{Any,Any}}, Qmax_dict::Dict{String, Float64}, Qmin_dict::Dict{String,Float64}, Hmax_dict::Dict{String, Float64}, Hmin_dict::Dict{String,Float64}, tight_coef::Float64, Q_lb::Float64)
    num_controls = length(control_names)
    for name in control_names
        if contains(controls_dict[name][:_condition][:name], "SimTime") == false && contains(controls_dict[name][:_condition][:name], "ClockTime") == false
            source_obj = controls_dict[name][:_condition][:_source_obj]
            source_obj_name = source_obj[:name]
            tank = wn_dict["tanks"][source_obj_name]
            interval = tank["max_level"] - tank["min_level"]
            cond1 = controls.ValueCondition(source_obj, "level", ">", tank["max_level"] - interval*0.05)
            cond2 = controls.ValueCondition(source_obj, "level", "<", tank["min_level"] + interval*0.05)

            target_obj = controls_dict[name][:_then_actions][1][:_target_obj]
            target_obj_name = target_obj[:name]
            act1 = controls.ControlAction(target_obj, "status", 0)
            act2 = controls.ControlAction(target_obj, "status", 1)

            control1 = controls.Control(cond1, act1)
            control2 = controls.Control(cond2, act2)

            wn[:add_control]("$name-1", control1)
            wn[:add_control]("$name-2", control2)
        end
    end
    UpdateExtrema(wn, wn_dict, Qmax_dict, Qmin_dict, Hmax_dict, Hmin_dict, 2, 2, tight_coef, Q_lb)
end

function tank_level(wn::PyCall.PyObject, wn_dict::Dict{Any,Any}, Qmax_dict::Dict{String, Float64}, Qmin_dict::Dict{String,Float64}, Hmax_dict::Dict{String, Float64}, Hmin_dict::Dict{String,Float64}, tight_coef::Float64, Q_lb::Float64)
    num_tanks = length(wn_dict["tanks"])
    num_sim = 2*num_tanks

    for sim = 1:num_sim
        switch = bin(sim, num_tanks)
        for (i, name) in enumerate(keys(wn_dict["tanks"]))
            switch[i] == '0' ? wn[:get_node](name)[:init_level] = wn_dict["tanks"][name]["max_level"] : wn[:get_node](name)[:init_level] = wn_dict["tanks"][name]["min_level"]
        end
        UpdateExtrema(wn, wn_dict, Qmax_dict, Qmin_dict, Hmax_dict, Hmin_dict, num_sim, sim, tight_coef, Q_lb)
    end
end

function tank_triggers_rand(wn::PyCall.PyObject, wn_dict::Dict{Any,Any}, controls_dict::Dict{Any,Any}, control_names::Base.KeyIterator{Dict{Any,Any}}, Qmax_dict::Dict{String, Float64}, Qmin_dict::Dict{String,Float64}, Hmax_dict::Dict{String, Float64}, Hmin_dict::Dict{String,Float64}, tight_coef::Float64, Q_lb::Float64)
    controls_dict = wn_dict["controls"]
    control_names = keys(controls_dict)
    for name in control_names
        if contains(controls_dict[name][:_condition][:name], "SimTime") == false && contains(controls_dict[name][:_condition][:name], "ClockTime") == false
            source_obj = controls_dict[name][:_condition][:_source_obj]
            source_obj_name = source_obj[:name]
            tank = wn_dict["tanks"][source_obj_name]
            interval = tank["max_level"] - tank["min_level"]

            rand_coefs = rand(2)
            sort!(rand_coefs)
            cond1 = controls.ValueCondition(source_obj, "level", ">", tank["min_level"]+ interval*rand_coefs[2])
            cond2 = controls.ValueCondition(source_obj, "level", "<", tank["min_level"]+ interval*rand_coefs[1])

            target_obj = controls_dict[name][:_then_actions][1][:_target_obj]
            target_obj_name = target_obj[:name]
            act1 = controls.ControlAction(target_obj, "status", 0)
            act2 = controls.ControlAction(target_obj, "status", 1)

            control1 = controls.Control(cond1, act1)
            control2 = controls.Control(cond2, act2)

            wn[:add_control]("$name-1", control1)
            wn[:add_control]("$name-2", control2)
        end
    end
    UpdateExtrema(wn, wn_dict, Qmax_dict, Qmin_dict, Hmax_dict, Hmin_dict, 2, 2, tight_coef, Q_lb)
end

function tank_level_rand(wn::PyCall.PyObject, wn_dict::Dict{Any,Any}, Qmax_dict::Dict{String, Float64}, Qmin_dict::Dict{String,Float64}, Hmax_dict::Dict{String, Float64}, Hmin_dict::Dict{String,Float64}, tight_coef::Float64, Q_lb::Float64)

    num_tanks = length(wn_dict["tanks"])
    num_sim = 2*num_tanks

    for sim = 1:num_sim
        for (i, name) in enumerate(keys(wn_dict["tanks"]))
            rand_coef = rand(1)
            wn[:get_node](name)[:init_level] = wn_dict["tanks"][name]["min_level"] + rand_coef[1]*(wn_dict["tanks"][name]["max_level"]- wn_dict["tanks"][name]["min_level"])
        end
        UpdateExtrema(wn, wn_dict, Qmax_dict, Qmin_dict, Hmax_dict, Hmin_dict, num_sim, sim, tight_coef, Q_lb)
    end
end

# function remove_controls(wn::PyCall.PyObject)
#     for name in wn[:control_name_list]
#         wn[:remove_control](name)
#     end
# end
#
# function run_sims(wn::PyCall.PyObject, wn_dict::Dict{Any,Any}, Qmax_dict::Dict{String, Float64}, Qmin_dict::Dict{String,Float64}, num_sim::Int64, sim::Int64, tight_coef::Float64, Q_lb::Float64)
#     wns = sims.EpanetSimulator(wn)
#     results = wns[:run_sim]()
#     link_results = results[:link]
#     pipes = wn_dict["pipes"]
#     valves = wn_dict["valves"]
#     links = merge(pipes, valves)
#     for (key, link) in links
#         name = link["name"]
#         flows = link_results["flowrate"][name][:values]
#         convert(Array{Float64}, flows)
#         if sim == 1
#             Qmin_dict[name] = convert(Float64,minimum(flows))
#             Qmax_dict[name] = convert(Float64, maximum(flows))
#         elseif sim <= num_sim
#             Qmin = Qmin_dict[name]
#             Qmax = Qmax_dict[name]
#
#             Qmin_dict[name] = minimum([Qmin, minimum(flows)])
#             Qmax_dict[name] = maximum([Qmax, maximum(flows)])
#         else
#             Qmax_val = Qmax_dict[name]
#             Qmin_val = Qmin_dict[name]
#             gap = Qmax_val - Qmin_val
#             if Qmin_val <= 0
#                 Qmin_dict[name] = Qmin_val - tight_coef*gap
#             else
#                 Qmin_dict[name] = minimum([Q_lb, Qmin_val])
#             end
#             if Qmax_val >= 0
#                 Qmax_dict[name] = Qmax_val + tight_coef*gap
#             else
#                 Qmax_dict[name] = maximum([Qmax_val, -Q_lb])
#             end
#         end
#     end
#
#     for (key, pump) in wn_dict["pumps"]
#         name = pump["name"]
#         flow = link_results["flowrate"][name][:values]
#         convert(Array{Float64}, flow)
#         if sim != 1
#             flow = flow[flow .> 0]
#         end
#
#         if length(flow) !=0 || sim == num_sim
#             if sim == 1
#                 Qmin_dict[name] = minimum(flow)
#                 Qmax_dict[name]= maximum(flow)
#             elseif sim <= num_sim
#                 Qmax = Qmax_dict[name]
#                 Qmin = Qmin_dict[name]
#                 Qmin_dict[name] = minimum([Qmin, minimum(flow)])
#                 Qmax_dict[name] = maximum([Qmax, maximum(flow)])
#             else
#                 Qmin_val = Qmin_dict[name]
#                 Qmax_val = Qmax_dict[name]
#
#                 gap = Qmax_val - Qmin_val
#                 if pump["pump_type"] == "HEAD"
#                     pump_curve_name = pump["pump_curve_name"]
#                     pump_curve = wn_dict["curves"][pump_curve_name]["points"]
#                 else
#                     pump_curve = [(pump["power"],0.0)] #power pump types gives fixed power value,
#                 end
#
#                 Qmin_dict[name] = maximum([pump_curve[1][1], Qmin_val - tight_coef*gap])
#                 Qmax_dict[name] = minimum([pump_curve[length(pump_curve)][1], Qmax_val + tight_coef*gap])
#             end
#         end
#     end
# end
# else
#     Qmin_val = Qmin_dict[name]
#     Qmax_val = Qmax_dict[name]
#
#     gap = Qmax_val - Qmin_val
#
#     pump_curve_name = pump["pump_curve_name"]
#     pump_curve = wn["curves"][pump_curve_name]["points"]
#     Qmin[name] = maximum(pump_curve[1][1], Qmin_val - tight_coef*gap)
#     Qmax[name] = minimum(pump_curve[length(pump_curve)][1], Qmax_val +tight_coef*gap)
# end

#     else
#         Qmax_val = Qmax_dict[name]
#         Qmin_val = Qmin_dict[name]
#         gap = Qmax_val - Qmin_val
#         if Qmin_val <= 0
#             Qmin_dict[name] = Qmin_val - tight_coef*gap
#         else
#             Qmin_dict[name] = minimum(Q_lb, Qmin_val)
#         end
#         if Qmax_val >= 0
#             Qmax_dict[name] = Qmax_val + tight_coef*gap
#         else
#             Qmax_dict[name] = maximum(Qmax_val, -Q_lb)
#         end
#     end
# end
# elseif sim == num_sim
#     for (key, tank) in wn_dict["tanks"]
#         wn[:get_node](key)[:init_level] = tank["init_level"]
#     end
#     control_names = wn[:control_name_list]
#     for name in control_names
#         println(name)
#         wn[:remove_control](name)
#     end
#     for (key, control) in wn_dict["controls"]
#         println("key $key")
#         wn[:add_control](key, control)
#     end
# end
# else
#     controls_dict = wn_dict["controls"]
#     control_names = keys(controls_dict)
#     for name in control_names
#         println(name)
#         wn[:remove_control](name)
#     end
#     for name in control_names
#         source_obj = controls_dict[name][:_condition][:_source_obj]
#         source_obj_name = source_obj[:name]
#         tank = wn_dict["tanks"][source_obj_name]
#         interval = tank["max_level"] - tank["min_level"]
#
#         rand_coefs = rand(2)
#         sort!(rand_coefs)
#         cond1 = controls.ValueCondition(source_obj, "level", ">", tank["min_level"]+ interval*rand_coefs[2])
#         cond2 = controls.ValueCondition(source_obj, "level", "<", tank["min_level"]+ interval*rand_coefs[1])
#
#         target_obj = controls_dict[name][:_then_actions][1][:_target_obj]
#         target_obj_name = target_obj[:name]
#         act1 = controls.ControlAction(target_obj, "status", 0)
#         act2 = controls.ControlAction(target_obj, "status", 1)
#
#         control1 = controls.Control(cond1, act1)
#         control2 = controls.Control(cond2, act2)
#
#         wn[:add_control]("$name-1", control1)
#         wn[:add_control]("$name-2", control2


text:
## Linear algebra API for adjacency matrix
using LinearAlgebra

function adjacency_matrix(adj::AbstractMatrix, T::DataType=eltype(adj))
    m, n = size(adj)
    (m == n) || throw(DimensionMismatch("adjacency matrix is not a square matrix: ($m, $n)"))
    T.(adj)
end

"""
    degrees(g[, T; dir=:out])

Degree of each vertex. Return a vector which contains the degree of each vertex in graph `g`.

# Arguments
- `g`: should be a adjacency matrix, `SimpleGraph`, `SimpleDiGraph` (from LightGraphs) or `SimpleWeightedGraph`, `SimpleWeightedDiGraph` (from SimpleWeightedGraphs).
- `T`: result element type of degree vector; default is the element type of `g` (optional).
- `dir`: direction of degree; should be `:in`, `:out`, or `:both` (optional).

# Examples
```jldoctest
julia> using GeometricFlux

julia> m = [0 1 1; 1 0 0; 1 0 0];

julia> GeometricFlux.degrees(m)
3-element Array{Int64,1}:
 2
 1
 1

```
"""
function degrees(adj::AbstractMatrix, T::DataType=eltype(adj); dir::Symbol=:out)
    if issymmetric(adj)
        d = vec(sum(adj, dims=1))
    else
        if dir == :out
            d = vec(sum(adj, dims=1))
        elseif dir == :in
            d = vec(sum(adj, dims=2))
        elseif dir == :both
            d = vec(sum(adj, dims=1)) + vec(sum(adj, dims=2))
        else
            throw(DomainError(dir, "invalid argument, only accept :in, :out and :both"))
        end
    end
    d
end

"""
    degree_matrix(g[, T; dir=:out])

Degree matrix of graph `g`. Return a matrix which contains degrees of each vertex in its diagonal.
The values other than diagonal are zeros.

# Arguments
- `g`: should be a adjacency matrix, `FeaturedGraph`, `SimpleGraph`, `SimpleDiGraph` (from LightGraphs) or `SimpleWeightedGraph`, `SimpleWeightedDiGraph` (from SimpleWeightedGraphs).
- `T`: result element type of degree vector; default is the element type of `g` (optional).
- `dir`: direction of degree; should be `:in`, `:out`, or `:both` (optional).

# Examples
```jldoctest
julia> using GeometricFlux

julia> m = [0 1 1; 1 0 0; 1 0 0];

julia> GeometricFlux.degree_matrix(m)
3×3 SparseArrays.SparseMatrixCSC{Int64,Int64} with 3 stored entries:
  [1, 1]  =  2
  [2, 2]  =  1
  [3, 3]  =  1
```
"""
function degree_matrix(adj::AbstractMatrix, T::DataType=eltype(adj); dir::Symbol=:out)
    d = degrees(adj, T, dir=dir)
    return SparseMatrixCSC(T.(diagm(0=>d)))
end

"""
    inv_sqrt_degree_matrix(g[, T; dir=:out])

Inverse squared degree matrix of graph `g`. Return a matrix which contains inverse squared degrees of each vertex in its diagonal.
The values other than diagonal are zeros.

# Arguments
- `g`: should be a adjacency matrix, `FeaturedGraph`, `SimpleGraph`, `SimpleDiGraph` (from LightGraphs) or `SimpleWeightedGraph`, `SimpleWeightedDiGraph` (from SimpleWeightedGraphs).
- `T`: result element type of degree vector; default is the element type of `g` (optional).
- `dir`: direction of degree; should be `:in`, `:out`, or `:both` (optional).
"""
function inv_sqrt_degree_matrix(adj::AbstractMatrix, T::DataType=eltype(adj); dir::Symbol=:out)
    d  = inv.(sqrt.(degrees(adj, T, dir=dir)))
    return Diagonal(d)
end

"""
    laplacian_matrix(g[, T; dir=:out])

Laplacian matrix of graph `g`.

# Arguments
- `g`: should be a adjacency matrix, `FeaturedGraph`, `SimpleGraph`, `SimpleDiGraph` (from LightGraphs) or `SimpleWeightedGraph`, `SimpleWeightedDiGraph` (from SimpleWeightedGraphs).
- `T`: result element type of degree vector; default is the element type of `g` (optional).
- `dir`: direction of degree; should be `:in`, `:out`, or `:both` (optional).
"""
function laplacian_matrix(adj::AbstractMatrix, T::DataType=eltype(adj); dir::Symbol=:out)
    degree_matrix(adj, T, dir=dir) - SparseMatrixCSC(T.(adj))
end

"""
    normalized_laplacian(g[, T; selfloop=false])

Normalized Laplacian matrix of graph `g`.

# Arguments
- `g`: should be a adjacency matrix, `FeaturedGraph`, `SimpleGraph`, `SimpleDiGraph` (from LightGraphs) or `SimpleWeightedGraph`, `SimpleWeightedDiGraph` (from SimpleWeightedGraphs).
- `T`: result element type of degree vector; default is the element type of `g` (optional).
- `selfloop`: adding self loop while calculating the matrix (optional).
"""
function normalized_laplacian(adj::AbstractMatrix, T::DataType=eltype(adj); selfloop::Bool=false)
    selfloop && (adj += I)
    inv_sqrtD = inv_sqrt_degree_matrix(adj, T, dir=:both)
    T.(I - inv_sqrtD * adj * inv_sqrtD)
end

@doc raw"""
    scaled_laplacian(adj::AbstractMatrix[, T::DataType])

Scaled Laplacien matrix of graph `g`,
defined as ``\hat{L} = \frac{2}{\lambda_{max}} L - I`` where ``L`` is the normalized Laplacian matrix.

# Arguments
- `g`: should be a adjacency matrix, `FeaturedGraph`, `SimpleGraph`, `SimpleDiGraph` (from LightGraphs) or `SimpleWeightedGraph`, `SimpleWeightedDiGraph` (from SimpleWeightedGraphs).
- `T`: result element type of degree vector; default is the element type of `g` (optional).
"""
function scaled_laplacian(adj::AbstractMatrix, T::DataType=eltype(adj))
    T(2. / eigmax(adj)) * normalized_laplacian(adj, T) - I
end

function neighbors(adj::AbstractMatrix, T::DataType=eltype(adj))
    n = size(adj,1)
    @assert n == size(adj,2) "adjacency matrix is not a square matrix."
    A = (adj .!= zero(T))
    if !issymmetric(adj)
        A = A .| A'
    end
    indecies = collect(1:n)
    ne = Vector{Int}[indecies[view(A, :, i)] for i = 1:n]
    return ne
end


text:
function runSimulation(p::precipParam)
  doInitialPrinting(p)
  if p.runID == 1
    mainDriver(p::precipParam)
  elseif p.runID == 2
    pureDiffusionOneChemicalTestDriver(p)
  elseif p.runID == 3
    pureAdvectionOneChemicalTestDriver(p)
  elseif p.runID == 4
    AdvectionDiffusionOneChemicalTestDriver(p)
  elseif p.runID == 5
    MultiphaseBrinkmanTestDriver(p)
  elseif p.runID == 6
    DiffusionAndReactionTwoChemicalTestDriver(p)
  elseif p.runID == 7
    pureReactionTestDriver(p)
  else
    println("runID of ",p.runID," not recognized.")
    println("runID 1-7 required")
  end

  nothing
end

function mainDriver(p::precipParam)
  timer, mesh, soln, tsoln, prob = initializeData(p)

  # initial brinkman
  mpbparam = BrinkmanMPParam(1.0,mpbFriction(soln.θf.u),soln.θf.u)
  sol = solve(prob.fluid,mesh.fMesh,mpbparam)
  tsoln.fluid.u = sol.u
  tsoln.fluid.v = sol.v
  tsoln.fluid.p = sol.p

  intializeTimeSteppingTimer!(timer)
  doPrinting!(timer,p,mesh,tsoln)
  while timer.simulationTime < timer.finalTime
    #Ninterp = 80
    #interpVals = zeros(Ninterp,6)
    #surf = mySurface([0.5,0.37],[0.5,-0.37])

    # reaction
    computeReaction!(p,tsoln)
    soln.θs.u = tsoln.θs.u
    soln.θf.u = tsoln.θf.u

    θfn = normalize_thetaf(tsoln.θf.u)
    peAnormalized = p.PeA./(θfn .+ 0.0 .*(1 .- θfn))        # 0.0 -> 0.5
    peBnormalized = p.PeB./(θfn .+ 0.0 .*(1 .- θfn))        # 0.0 -> 0.25
    ψAParam = AdvDiffParam(tsoln.fluid.u,tsoln.fluid.v,peAnormalized)
    ψBParam = AdvDiffParam(tsoln.fluid.u,tsoln.fluid.v,peBnormalized)
    ψCParam = AdvDiffParam(tsoln.fluid.u,tsoln.fluid.v,p.PeC./θfn)

    # advection-diffusion
    M = generateMassMatrixWScalar(mesh.sMesh,θfn)

    StiffA = GenerateSystem(mesh.sMesh,prob.ψA,ψAParam)
    ApplyBC!(StiffA,mesh.sMesh,prob.ψA,ψAParam,prob.ψA.OperatorType)

    StiffB = GenerateSystem(mesh.sMesh,prob.ψB,ψBParam)
    ApplyBC!(StiffB,mesh.sMesh,prob.ψB,ψBParam,prob.ψB.OperatorType)

    StiffC = GenerateSystem(mesh.sMesh,prob.ψC,ψCParam)
    ApplyBC!(StiffC,mesh.sMesh,prob.ψC,ψCParam,prob.ψC.OperatorType)

    if timer.stepIndex < 5
      soln.ψA.u = updateChemical(M,StiffA,tsoln.ψA.u,p;solver=:backwardeuler)
      soln.ψB.u = updateChemical(M,StiffB,tsoln.ψB.u,p;solver=:backwardeuler)
      soln.ψC.u = updateChemical(M,StiffC,tsoln.ψC.u,p;solver=:backwardeuler)
    else
      soln.ψA.u = updateChemical(M,StiffA,tsoln.ψA.u,p;solver=p.integrator)
      soln.ψB.u = updateChemical(M,StiffB,tsoln.ψB.u,p;solver=p.integrator)
      soln.ψC.u = updateChemical(M,StiffC,tsoln.ψC.u,p;solver=p.integrator)
    end

    # scrub scalars so they don't go below zero.
    for i=1:length(soln.ψA.u)
      if soln.ψA.u[i] < 0
        soln.ψA.u[i] = 0.0
      end
      if soln.ψB.u[i] < 0
        soln.ψB.u[i] = 0.0
      end
      if soln.ψC.u[i] < 0
        soln.ψC.u[i] = 0.0
      end
    end

    # brinkman
    # normalize thetaf
    mpbparam = BrinkmanMPParam(1.0,mpbFriction(θfn),θfn)
    sol = solve(prob.fluid,mesh.fMesh,mpbparam)
    soln.fluid.u = sol.u
    soln.fluid.v = sol.v
    soln.fluid.p = sol.p

    # inteprolate values
    #xvals = zeros(length(mesh.sMesh.xy))
    #yvals = zeros(length(mesh.sMesh.xy))
    #for i=1:length(mesh.sMesh.xy)
    #  xvals[i] = mesh.sMesh.xy[i].x
    #  yvals[i] = mesh.sMesh.xy[i].y
    #end
    #vmag = sqrt.(soln.fluid.u.^2 + soln.fluid.v.^2)
    #interpVals[:,1] = SurfaceInterp(mesh.sMesh,xvals,surf,Ninterp-1)
    #interpVals[:,2] = SurfaceInterp(mesh.sMesh,yvals,surf,Ninterp-1)
    #interpVals[:,3] = SurfaceInterp(mesh.sMesh,soln.ψA.u,surf,Ninterp-1)
    #interpVals[:,4] = SurfaceInterp(mesh.sMesh,soln.ψB.u,surf,Ninterp-1)
    #interpVals[:,5] = SurfaceInterp(mesh.sMesh,soln.θs.u,surf,Ninterp-1)
    #interpVals[:,6] = SurfaceInterp(mesh.fMesh,vmag,surf,Ninterp-1)

    tsoln = deepcopy(soln)
    if isPrintIndex(timer)
      # save interpolated values
      #savefile = string(p.SaveFolder,"/interpVals_",timer.printIndex,".jld2")
      #save(savefile, "time", timer.simulationTime, "interpVals", interpVals)

      # ordinary printing
      doPrinting!(timer,p,mesh,soln)
    end
    updateSimulationTime!(timer)
    updateStepIndex!(timer)
  end

  printSimulationFinishedMessage(timer)

  nothing
end

function normalize_thetaf(θf::Vector{T}) where T<:Real
  N=length(θf)
  θf_normalized = zeros(T,N)
  TOL = 0.1

  for i=1:N
    if θf[i] < TOL
      θf_normalized[i] = TOL
    else
      θf_normalized[i] = θf[i]
    end
  end

  return θf_normalized
end

function pureDiffusionOneChemicalTestDriver(p::precipParam)
    timer, mesh, soln, tsoln, prob = initializeData(p)

    ψAParam = PoissonParam(p.κA)

    for i=1:length(soln.θf.u)
      soln.θf.u[i] = 0.1
    end

    M = generateMassMatrixWScalar(mesh.sMesh,soln.θf.u)
    StiffA = GenerateSystem(mesh.sMesh,prob.ψA,ψAParam)
    ApplyBC!(StiffA,mesh.sMesh,prob.ψA,ψAParam,prob.ψA.OperatorType)

    intializeTimeSteppingTimer!(timer)
    doPrinting!(timer,p,mesh,soln)
    while timer.simulationTime < timer.finalTime
      if timer.stepIndex < 5
        soln.ψA.u = updateChemical(M,StiffA,tsoln.ψA.u,p;solver=:backwardeuler)
      else
        soln.ψA.u = updateChemical(M,StiffA,tsoln.ψA.u,p;solver=p.integrator)
      end

      tsoln = deepcopy(soln)
      if isPrintIndex(timer)
        doPrinting!(timer,p,mesh,soln)
      end
      updateSimulationTime!(timer)
      updateStepIndex!(timer)
    end

    printSimulationFinishedMessage(timer)

    nothing
end

function pureAdvectionOneChemicalTestDriver(p::precipParam)
  println("pure advection test driver isn't programmed yet")
  nothing
end

function AdvectionDiffusionOneChemicalTestDriver(p::precipParam)
  println("advection-diffusion test driver isn't programmed yet")
  nothing
end

function MultiphaseBrinkmanTestDriver(p::precipParam)
  function computeNorms(N)
    start = time()
    # load mesh
    mesh = squareMeshFluid([-2,2,-1,1],N)

    α1(x,y) = 1.0
    α2(x,y) = 0.5+0.4*cos(3*x*y)
    α3(x,y) = 0.5+0.4*sin(7*x*y)

    xm = [i.x for i in mesh.xy]
    ym = [i.y for i in mesh.xy]
    α1arr = [α1(xm[i],ym[i]) for i=1:length(mesh.xy)]
    α2arr = [α2(xm[i],ym[i]) for i=1:length(mesh.xy)]
    α3arr = [α3(xm[i],ym[i]) for i=1:length(mesh.xy)]

    param  = BrinkmanMPParam(α1arr,α2arr,α3arr)

    OperatorType = :BrinkmanMP2D

    dNodes = Dirichlet(:left,:top,:bottom,:right)
    Nodes = [dNodes]

    u(x,y) = -x^4*y^2
    v(x,y) = 4x^3*y^3/3

    dUBCf = Dirichlet( (x,y) -> u(x,y) )
    dVBCf = Dirichlet( (x,y) -> v(x,y) )

    d2dudx(x,y) = -12*x^2*y^2; d2dudy(x,y) = -2*x^4
    d2dvdx(x,y) = 8*x*y^3;     d2dvdy(x,y) = 8*x^3*y
    lapu(x,y) = d2dudx(x,y) + d2dudy(x,y)
    lapv(x,y) = d2dvdx(x,y) + d2dvdy(x,y)

    dpdx(x,y) = 3*x^2*y^3; dpdy(x,y) = 3*x^3*y^2

    Fx(x,y) = -α1(x,y)*lapu(x,y) + α2(x,y)*u(x,y) + α3(x,y)*dpdx(x,y)
    Fy(x,y) = -α1(x,y)*lapv(x,y) + α2(x,y)*v(x,y) + α3(x,y)*dpdy(x,y)

    ffx   = Forcing(Fx)
    ffy   = Forcing(Fy)
    bcfun = [dUBCf,dVBCf,ffx,ffy]

    prob = Problem(mesh,Nodes,bcfun,OperatorType)
    sol = solve(prob,mesh,param)

    elapsed = time()-start

    # compute condition number of operator matrix
    #LinOp = GenerateSystem(mesh,prob,param)
    #ApplyBC!(LinOp,mesh,prob,param,OperatorType)
    #κ = cond(Array(LinOp.Op),2)
    κ = 0.5

    # compute mesh h
    h = hCalc(mesh)

    # generate exact solution
    Uexact(x,y) = dUBCf.f(x,y)
    Vexact(x,y) = dVBCf.f(x,y)
    Pexact(x,y) = x^3*y^3
    UexactArr = [Uexact.(mesh.xy[i].x,mesh.xy[i].y) for i=1:length(mesh.xy)]
    VexactArr = [Vexact.(mesh.xy[i].x,mesh.xy[i].y) for i=1:length(mesh.xy)]
    PexactArr = [Pexact.(mesh.xyp[i].x,mesh.xyp[i].y) for i=1:length(mesh.xyp)]

    # compute difference
    velErr = sqrt.((sol.u - UexactArr).^2 + (sol.v -VexactArr).^2)
    presErr = sol.p - PexactArr

    # compute Li norm
    L1v   = DomainNorm(mesh.xy,mesh.cm,velErr;normID="1")
    L2v   = DomainNorm(mesh.xy,mesh.cm,velErr;normID="2")
    Linfv = DomainNorm(mesh.xy,mesh.cm,velErr;normID="Inf")
    L1p   = DomainNorm(mesh.xyp,mesh.cmp,presErr;normID="1")
    L2p   = DomainNorm(mesh.xyp,mesh.cmp,presErr;normID="2")
    Linfp = DomainNorm(mesh.xyp,mesh.cmp,presErr;normID="Inf")

    # plot solution
    sD = ScalarData(sol.p,α1arr,α2arr,α3arr)
    sN = ScalarNames("pressure","alpha_1","alpha_2","alpha_3")
    vD = VectorData([sol.u,sol.v])
    vN = VectorNames("velocity")
    fn = Path("solution_var")
    vtksave(mesh,sD,sN,vD,vN,fn)

    # plot solution
    sD = ScalarData(PexactArr)
    sN = ScalarNames("pressure")
    vD = VectorData([UexactArr,VexactArr])
    vN = VectorNames("velocity")
    fn = Path("exact_var")
    vtksave(mesh,sD,sN,vD,vN,fn)

    return κ,h,L1v,L2v,Linfv,L1p,L2p,Linfp,elapsed
  end

  Narr     = [4,4,8,16,32,64]#,96]
  N        = length(Narr)
  harr     = zeros(N)
  L1arrV   = zeros(N)
  L2arrV   = zeros(N)
  LInfarrV = zeros(N)
  L1arrP   = zeros(N)
  L2arrP   = zeros(N)
  LInfarrP = zeros(N)
  timearr  = zeros(N)
  κarr     = zeros(N)

  for i=1:N
    n = Narr[i]
    κarr[i],harr[i],L1arrV[i],L2arrV[i],LInfarrV[i],
      L1arrP[i],L2arrP[i],LInfarrP[i],timearr[i] =
      computeNorms(n)
      println("completed N=$(n)")
  end

  κarr     = κarr[2:end]
  harr     = harr[2:end]
  L1arrV   = L1arrV[2:end]
  L2arrV   = L2arrV[2:end]
  LInfarrV = LInfarrV[2:end]
  L1arrP   = L1arrP[2:end]
  L2arrP   = L2arrP[2:end]
  LInfarrP = LInfarrP[2:end]
  timearr  = timearr[2:end]

  # export output to *.jld file
  run(`mkdir -p data`)

  save("data/TEMP_mpb2D_validation.jld", "harr",harr,
                                          "VL1arr", L1arrV,
                                          "VL2arr", L2arrV,
                                          "VLInfarr",LInfarrV,
                                          "PL1arr", L1arrP,
                                          "PL2arr", L2arrP,
                                          "PLInfarr",LInfarrP,
                                          "timearr",timearr,
                                          "κarr",κarr)

  nothing
end

function DiffusionAndReactionTwoChemicalTestDriver(p::precipParam)
  timer, mesh, soln, tsoln, prob = initializeData(p)

  ψAParam = PoissonParam(p.κA)
  ψBParam = PoissonParam(p.κB)
  ψCParam = PoissonParam(p.κC)

  M = generateMassMatrix(mesh.sMesh)

  StiffA = GenerateSystem(mesh.sMesh,prob.ψA,ψAParam)
  ApplyBC!(StiffA,mesh.sMesh,prob.ψA,ψAParam,prob.ψA.OperatorType)

  StiffB = GenerateSystem(mesh.sMesh,prob.ψB,ψBParam)
  ApplyBC!(StiffB,mesh.sMesh,prob.ψB,ψBParam,prob.ψB.OperatorType)

  StiffC = GenerateSystem(mesh.sMesh,prob.ψC,ψCParam)
  ApplyBC!(StiffC,mesh.sMesh,prob.ψC,ψCParam,prob.ψC.OperatorType)

  intializeTimeSteppingTimer!(timer)
  doPrinting!(timer,p,mesh,soln)
  while timer.simulationTime < timer.finalTime
    println("i1=",timer.stepIndex," ",computeTotalMass(p,mesh.sMesh,soln))
    if timer.stepIndex < 5
      soln.ψA.u = updateChemical(M,StiffA,tsoln.ψA.u,p;solver=:backwardeuler)
      soln.ψB.u = updateChemical(M,StiffB,tsoln.ψB.u,p;solver=:backwardeuler)
      soln.ψC.u = updateChemical(M,StiffC,tsoln.ψC.u,p;solver=:backwardeuler)
    else
      soln.ψA.u = updateChemical(M,StiffA,tsoln.ψA.u,p;solver=p.integrator)
      soln.ψB.u = updateChemical(M,StiffB,tsoln.ψB.u,p;solver=p.integrator)
      soln.ψC.u = updateChemical(M,StiffC,tsoln.ψC.u,p;solver=p.integrator)
    end
    println("i2=",timer.stepIndex," ",computeTotalMass(p,mesh.sMesh,soln))
    computeReaction!(p,soln)
    println("i3=",timer.stepIndex," ",computeTotalMass(p,mesh.sMesh,soln))
    println()

    tsoln = deepcopy(soln)
    if isPrintIndex(timer)
      doPrinting!(timer,p,mesh,soln)
    end
    updateSimulationTime!(timer)
    updateStepIndex!(timer)
  end

  printSimulationFinishedMessage(timer)

  nothing
end

function pureReactionTestDriver(p::precipParam)
  timer, mesh, soln, tsoln, prob = initializeData(p)

  randomizeScalarSolutions!(soln)

  intializeTimeSteppingTimer!(timer)
  doPrinting!(timer,p,mesh,soln)
  while timer.simulationTime < timer.finalTime
    computeReaction!(p,soln)

    if isPrintIndex(timer)
      doPrinting!(timer,p,mesh,soln)
    end

    tsoln = deepcopy(soln)

    # step time
    updateSimulationTime!(timer)
    updateStepIndex!(timer)
  end

  printSimulationFinishedMessage(timer)

  nothing
end

function initializeData(p::precipParam)
  timer = initializeTimer(p)
  mesh = initializeMeshes(p)
  soln = initializeSolution(p)
  tsoln = initializeSolution(p)
  prob = defineProblems(p,mesh)

  return timer, mesh, soln, tsoln, prob
end

function mpbFriction(θf::Vector{T}) where T<:Real
  N = length(θf)
  friction = zeros(N)

  K=0.5
  n=2
  ξstar = 30.0
  θstar = 0.5
  h = 1000.0#computeH(ξstar,θstar,K,n)

  for i=1:N
    friction[i] = h*frictionKernel(θf[i],K,n)
  end

  return friction
end

function frictionKernel(θf,K,n)
  θs = (1.0-θf)
  return θs^n/(K^n + θs^n)
end

function computeH(ξstar,θstar,K,n)
  return ξstar/frictionKernel(θstar,K,n)
end


text:
module BatteriesRechargePlatform

using ..Ahorn, Maple

@mapdef Entity "batteries/recharge_platform" BatteryRechargePlatform(x::Integer, y::Integer)

const placements = Ahorn.PlacementDict(
    "Battery Recharge Platform" => Ahorn.EntityPlacement(
        BatteryRechargePlatform
    )
)

sprite = "batteries/recharge_platform/base0"

function Ahorn.selection(entity::BatteryRechargePlatform)
    x, y = Ahorn.position(entity)

    return Ahorn.getSpriteRectangle(sprite, x, y - 4)
end

Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::BatteryRechargePlatform, room::Maple.Room) = Ahorn.drawSprite(ctx, sprite, 0, -4)

end

max_stars_repo_name: PenguinOwl/Batteries

max_issues_repo_name: PenguinOwl/Batteries

max_forks_repo_name: PenguinOwl/Batteries