text:
#!/usr/bin/env julia
include("seir_var.jl")
push!(PGFPlotsX.CUSTOM_PREAMBLE, raw"\usepackage{fontspec}")
push!(PGFPlotsX.CUSTOM_PREAMBLE, raw"\setmainfont{Asana-Math}")
pgfsave("../../images/seir_var_susp.svg",
    seir_var_susp(2.5, .2, .1, .25, 1e-5, 365, false) |> seirplot)
pgfsave("../../images/seir_var_k_r730.svg",
    seirset(2.5, .2, .1, -5:7, 1e-5, 730) |> seirset_rendsplot)

max_stars_repo_name: klpn/static-dust

max_issues_repo_name: klpn/static-dust

max_forks_repo_name: klpn/static-dust


text:
# Gets the logs messages associated with a AWSBatch BatchJob as a single string
function log_messages(job::BatchJob; retries=5, wait_interval=7)
    i = 0
    events = log_events(job)
    # Retry if no logs are found
    while i < retries && (events === nothing || isempty(events))
        i += 1
        sleep(wait_interval)
        events = log_events(job)
    end

    events === nothing && return ""
    return join([string(event.timestamp, "  ", event.message) for event in events], '\n')
end

function report(io::IO, job::BatchJob)
    println(io, "Job ID: $(job.id)")
    print(io, "Status: $(status(job))")

    reason = status_reason(job)
    reason !== nothing && print(io, " ($reason)")
    println(io)

    log_str = log_messages(job)
    return !isempty(log_str) && println(io, '\n', log_str)
end

report(job::BatchJob) = sprint(report, job)

function job_duration(job::BatchJob)
    d = describe(job)
    if haskey(d, "createdAt") && haskey(d, "stoppedAt")
        Millisecond(d["stoppedAt"] - d["createdAt"])
    else
        nothing
    end
end

function job_runtime(job::BatchJob)
    d = describe(job)
    if haskey(d, "startedAt") && haskey(d, "stoppedAt")
        Millisecond(d["stoppedAt"] - d["startedAt"])
    else
        nothing
    end
end

function time_str(secs::Real)
    @sprintf("%02d:%02d:%02d", div(secs, 3600), rem(div(secs, 60), 60), rem(secs, 60))
end

time_str(seconds::Second) = time_str(Dates.value(seconds))
time_str(p::Period) = time_str(floor(p, Second))
time_str(::Nothing) = "N/A"  # Unable to determine duration

function register_job_definition(job_definition::AbstractDict)
    output = Batch.register_job_definition(
        job_definition["jobDefinitionName"], job_definition["type"], job_definition
    )

    return output["jobDefinitionArn"]
end

function submit_job(;
    job_name::AbstractString,
    job_definition::AbstractString,
    job_queue::AbstractString=STACK["WorkerJobQueueArn"],
    node_overrides::Dict=Dict(),
    retry_strategy::Dict=Dict(),
)
    options = Dict{String,Any}()

    if !isempty(node_overrides)
        options["nodeOverrides"] = node_overrides
    end

    if !isempty(retry_strategy)
        options["retryStrategy"] = retry_strategy
    end

    print(JSON.json(options, 4))

    output = Batch.submit_job(job_definition, job_name, job_queue, options)
    return BatchJob(output["jobId"])
end

function describe_compute_environment(compute_environment::AbstractString)
    # Equivalent to running the following on the AWS CLI
    # ```
    # aws batch describe-compute-environments
    #   --compute-environments $(STACK["ComputeEnvironmentArn"])
    #   --query computeEnvironments[0]
    # ```
    output = Batch.describe_compute_environments(
        Dict("computeEnvironments" => [compute_environment])
    )

    details = if !isempty(output["computeEnvironments"])
        output["computeEnvironments"][1]
    else
        nothing
    end

    return details
end

function wait_finish(job::BatchJob; timeout::Period=Minute(20))
    timeout_secs = Dates.value(Second(timeout))

    info(LOGGER, "Waiting for AWS Batch job to finish (~5 minutes)")
    # TODO: Disable logging from wait? Or at least emit timestamps
    wait(job, [AWSBatch.FAILED, AWSBatch.SUCCEEDED]; timeout=timeout_secs)  # TODO: Support timeout as Period
    duration = job_duration(job)
    runtime = job_runtime(job)
    info(LOGGER, "Job duration: $(time_str(duration)), Job runtime: $(time_str(runtime))")

    return nothing
end

ext: jl


text:
#InspectDR: Event handlers for user input of Gtk layer
#-------------------------------------------------------------------------------

#==Types
===============================================================================#
const KeyMap = Dict{Int, Function}

mutable struct KeyBindings
	nomod::KeyMap
	shiftmod::KeyMap
	ctrlmod::KeyMap
	altmod::KeyMap
end
KeyBindings() = KeyBindings(KeyMap(), KeyMap(), KeyMap(), KeyMap())

#==Module-level definitions
===============================================================================#
const keybindings = KeyBindings()

#==Helper functions
===============================================================================#
function setstate_normal(pwidget::PlotWidget)
	gdk_window_set_cursor(pwidget.canvas, CURSOR_DEFAULT)
	pwidget.state = ISNormal()
end
ignoreuserinput(pwidget::PlotWidget) = nothing

#Handle mousepress if a CtrlElement was clicked:
function handleevent_mousepress(pwidget::PlotWidget, ::Type{CtrlElement}, x::Float64, y::Float64)
	istrip = hittest(pwidget, x, y)

	if istrip > 0
		if handleevent_mousepress(pwidget, pwidget.markers, istrip, x, y)
			return true
		end

		#TODO: test for other control points here....
	end

	return false
end

#Set focus on active strip under x,y coordinates:
function focus_strip(pwidget::PlotWidget, x::Float64, y::Float64)
	istrip = hittest(pwidget, x, y)
	if istrip > 0
		pwidget.curstrip = istrip
	end
	nothing
end

#==Default event handlers
===============================================================================#
handleevent_mousepress(::InputState, pwidget::PlotWidget, event::Gtk.GdkEventButton) =
	nothing
handleevent_mouserelease(::InputState, pwidget::PlotWidget, event::Gtk.GdkEventButton) =
	nothing
function handleevent_mousemove(::InputState, pwidget::PlotWidget, event::Gtk.GdkEventMotion)
	handleevent_plothover(pwidget, event)
	nothing
end

#Event handler for plot widget keypress
#-------------------------------------------------------------------------------
function handleevent_keypress(::InputState, pwidget::PlotWidget, event::Gtk.GdkEventKey)
#event.event_type: GDK_KEY_PRESS/GDK_KEY_RELEASE
#@show event.event_type, event.state, event.keyval, event.length, event.string, event.hardware_keycode, event.group, event.flags
	global keybindings

	#Find bound functions & execute:
	if modifiers_pressed(event.state, 0)
		fn = get(keybindings.nomod, event.keyval, ignoreuserinput)
		fn(pwidget)
	elseif modifiers_pressed(event.state, MODIFIER_SHIFT)
		fn = get(keybindings.shiftmod, event.keyval, ignoreuserinput)
		fn(pwidget)
	elseif modifiers_pressed(event.state, MODIFIER_CTRL)
		fn = get(keybindings.ctrlmod, event.keyval, ignoreuserinput)
		fn(pwidget)
	elseif modifiers_pressed(event.state, MODIFIER_ALT)
		fn = get(keybindings.altmod, event.keyval, ignoreuserinput)
		fn(pwidget)
	end
#	@show event.keyval
end

function handleevent_mousescroll(::InputState, pwidget::PlotWidget, event::Gtk.GdkEventScroll)
	if modifiers_pressed(event.state, 0)
		if GdkScrollDirection.GDK_SCROLL_UP == event.direction
			pan_up(pwidget)
		elseif GdkScrollDirection.GDK_SCROLL_DOWN == event.direction
			pan_down(pwidget)
		end
	elseif modifiers_pressed(event.state, MODIFIER_SHIFT)
		if GdkScrollDirection.GDK_SCROLL_UP == event.direction
			pan_left(pwidget)
		elseif GdkScrollDirection.GDK_SCROLL_DOWN == event.direction
			pan_right(pwidget)
		end
	elseif modifiers_pressed(event.state, MODIFIER_CTRL)
		if GdkScrollDirection.GDK_SCROLL_UP == event.direction
			zoom_in(pwidget, event.x, event.y)
		elseif GdkScrollDirection.GDK_SCROLL_DOWN == event.direction
			zoom_out(pwidget, event.x, event.y)
		end
	end
end
function handleevent_scalechanged(::InputState, pwidget::PlotWidget)
	if !scalectrl_enabled(pwidget); return; end #Should not really happen
	scalectrl_apply(pwidget)
end

#==State: Normal input state
===============================================================================#
function handleevent_mousepress(::ISNormal, pwidget::PlotWidget, event::Gtk.GdkEventButton)
#	@show event.state, event.button, event.event_type
	focus_strip(pwidget, event.x, event.y)
	focus(pwidget.widget) #In case not in focus

	if 3==event.button
		boxzoom_setstart(pwidget, event.x, event.y) #Changes state
	elseif 1==event.button
		if modifiers_pressed(event.state, MODIFIER_SHIFT)
			mousepan_setstart(pwidget, event.x, event.y) #Changes state
		elseif !modifiers_pressed(event.state, MODIFIERS_SUPPORTED) #Un-modified
			handleevent_mousepress(pwidget, CtrlElement, event.x, event.y)
		end
	end
end

#==Key bindings
===============================================================================#
function _setdefaults(bnd::KeyBindings)
	bnd.nomod = KeyMap(
		GdkKeySyms.Escape => setstate_normal,
		GdkKeySyms.Up => pan_up,
		GdkKeySyms.Down => pan_down,
		GdkKeySyms.Left => pan_left,
		GdkKeySyms.Right => pan_right,
		GdkKeySyms_Plus => zoom_in,
		GdkKeySyms_Minus => zoom_out,

		Int('r') => addrefmarker,
		Int('d') => addΔmarker,
	)
	bnd.shiftmod = KeyMap(
		Int('D') => addΔmarkerref,
	)
	bnd.ctrlmod = KeyMap(
		Int('f') => zoom_full,
		Int('h') => zoom_hfull,
		Int('v') => zoom_vfull,
	)
	bnd.altmod = KeyMap(
	)
end

#Last line


text:
function always_assert(cond::Bool,
                       msg::String = "")
    if !cond
        throw(AlwaysAssertionError(msg))
    end
    return nothing
end


text:
using PtFEM

old = pwd()
ProjDir = joinpath(dirname(@__FILE__), "..", "examples", "4 Static Equilibrium")
#cd(ProjDir)

include(joinpath(ProjDir, "p4.1.1.jl"))

#=
data = Dict(
  # Rod(nels, np_types, nip, finite_element(nod, nodof))
  :struc_el => Rod(4, 1, 1, Line(2, 1)),
  :properties => [1.0e5;],
  :x_coords => range(0, stop=1, length=5),
  :support => [(1, [0])],
  :loaded_nodes => [(1,[-0.625]),(2,[-1.25]),(3,[-1.25]),(4,[-1.25]),(5,[-0.625])]
)
data |> display
println()

@time m = FE4_1(data)
println()

println("Displacements:")
m.displacements |> display
println()

println("Actions:")
m.actions |> display
println()
=#

@show m.loads
println()

m1 = fromSkyline(m.kv, m.kdiag)
@show m1
println()

@show m.kdiag
println()

@show m.kv
println()

ms = skyline2sparse(m.kv, m.kdiag)
@show ms
println()

loads = zeros(m.neq+1)
if :loaded_nodes in keys(data)
  for i in 1:size(data[:loaded_nodes], 1)
    loads[m.nf[:, data[:loaded_nodes][i][1]]+1] = data[:loaded_nodes][i][2]
  end
end

@show loads
println()

@show m1 \ loads[2:end]
println()

@show ms \ loads[2:end]
println()
cd(old)


text:
module LatticeSweeper

using LinearAlgebra: diagm, svd
using Printf: @printf
using Statistics: mean

using Arpack: eigs
using LinearMaps: LinearMap
using OffsetArrays: OffsetArray
using TensorOperations: @tensor

export
    S_vn,
    println_result,

    MPS,
    MPO,
    compress!,

    SweepSchedule,
    SweepOutput,
    SweepOutputFile,
    SweepOutputDynamic,
    dmrg!

include("util.jl")

include("mps.jl")
include("mpo.jl")
include("contraction.jl")

include("settings.jl")
include("state.jl")
include("history.jl")
include("output.jl")

include("dmrg.jl")

end


text:
function sortselected(a::AbstractVector{<:Real}, s::AbstractVector{<:Integer})
    sel = unique(sort(s))
    if sel[1] < 1 || length(a) < sel[end]
        throw(BoundsError())
    end
    b = collect(copy(a))
    b[sel] = sort(b[sel])
    return b
end

a = [7, 6, 5, 4, 3, 2, 1, 0]
sel = [7, 2, 8]
b = sortselected(a, sel)

println("Original: $a\n\tsorted on $sel\n -> sorted array: $b")

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using InfiniteLinearAlgebra, LinearAlgebra, BandedMatrices, ArrayLayouts, Test

@testset "infinite-cholesky" begin
    S = Symmetric(BandedMatrix(0 => 1:∞, 1=> Ones(∞)))
    b = [1; zeros(∞)]
    @test cholesky(S) \ b ≈ qr(S) \ b ≈ S \ b

    # go past adaptive
    b = [randn(10_000); zeros(∞)]
    @test cholesky(S) \ b ≈ qr(S) \ b ≈ S \ b

    @testset "singularly perturbed" begin
        ε = 0.0001
        S = Symmetric(BandedMatrix(0 => 2 .+ ε*(1:∞), 1=> -Ones(∞)))
        b = [1; zeros(∞)]
        @test cholesky(S) \ b ≈  S \ b
    end

    @testset "long bandwidths" begin
        S = Symmetric(BandedMatrix(0 => 1:∞, 50=> Ones(∞)))
        b = [1; zeros(∞)]
        @test cholesky(S) \ b ≈ qr(S) \ b ≈ S \ b
    end

    @testset "powers" begin
        b = [1; zeros(∞)]
        @test cholesky(S^2) \ b ≈ qr(S^2) \ b ≈ S^2 \ b

        @test cholesky(S^2).U[1:100,1:100] ≈ cholesky(Symmetric((S^2)[1:100,1:100])).U
    end
end


text:
function check_permutation(str1 :: ASCIIString, str2 :: ASCIIString)
  if length(str1) != length(str2) # permutations must have the same length
    return false
  end

  d :: Dict{Char, Int} = Dict{Char, Int}()

  for ch in str1      # count characters in the first string
    if haskey(d, ch)
      d[ch] += 1
    else
      d[ch] = 1
    end
  end

  for ch in str2      # discount characters in the second string
    if haskey(d, ch)
      if d[ch] == 0   # if there is negative count to appear
        return false  # return false
      end
      d[ch] -= 1
    else              # if there is new characters
      return false    # return false
    end
  end

  return true
end

# Time  : O(n)
# Space : O(n)

using Base: Test

@test check_permutation("abcd", "bcda") == true
@test check_permutation("AbCd", "bdAC") == true
@test check_permutation("abcd", "bcdda") == false
@test check_permutation("abcd", "beda") == false

max_stars_repo_name: Kiandr/crackingcodinginterview

max_issues_repo_name: Kiandr/crackingcodinginterview

max_forks_repo_name: Kiandr/crackingcodinginterview


text:
# include("TrainingPages.jl")

HSD = "HarmonieSystemDocumentation"

Configuration = [
    "Experiment" => "$HSD/ConfigureYourExperiment.md",
    "Platform" => "$HSD/PlatformConfiguration.md",
    "Namelist" => "$HSD/Namelists.md",
    "Update Namelist" => "$HSD/UpdateNamelists.md",
    "Model Domain" => "$HSD/ModelDomain.md",
    "Vertical Grid" => "$HSD/VerticalGrid.md",
    "Climate Simulation" => "$HSD/ClimateSimulation.md",
    "Hires topography" => "$HSD/How_to_use_hires_topography.md"
]

ECFLOW_TASKS = [
    "Climate" => "$HSD/ClimateGeneration.md",
    "Observations" => [
        "$HSD/UseofObservation.md",
        "$HSD/ObservationOperators.md",
        "$HSD/ObservationHowto/GNSS.md",
        "$HSD/ObservationHowto/Amv.md",
        "$HSD/RadarData.md",
        "Preprocessing" => [
         # "ObservationPreprocessing.md" ,
            "$HSD/ObservationPreprocessing/ObservationData.md",
            "Bator" => "$HSD/ObservationPreprocessing/Bator.md",
            "Oulan" => "$HSD/ObservationPreprocessing/Oulan.md",
            "Cope" => "$HSD/ObservationPreprocessing/Cope.md"
         ]
    ],
    "Bator" => "$HSD/ObservationPreprocessing/Bator.md",
    "Oulan" => "$HSD/ObservationPreprocessing/Oulan.md",
    "Boundaries" => "$HSD/BoundaryFilePreparation.md",
    "Screening" => "$HSD/Screening.md",
    "Surface Analysis" => "$HSD/Analysis/SurfaceAnalysis.md",
    "DFI" => "$HSD/DigitialFilterInitialization.md",
    "Forecast" => "$HSD/Forecast.md",
    "Post Processing" => [
        "Obsmon" => "$HSD/PostPP/Obsmon.md",
        "Diagnostics" => "$HSD/PostPP/Diagnostics.md",
        "xtool" => "$HSD/PostPP/xtool.md",
        "GL interpolation" => "$HSD/PostPP/gl/Interpolation.md",
        "Verification" => [
            "$HSD/PostPP/Verification.md",
            "$HSD/PostPP/Extract4verification.md"
        ],
        "Fullpos" => "$HSD/PostPP/Fullpos.md",
        "FileConversion" => "$HSD/PostPP/FileConversions.md",
        "gl_grib_api" => "$HSD/PostPP/gl_grib_api.md",
    ]
]

pages = [
    # "Overview" => ["index" => ["index.md"]],

        "Running at ECMWF and local" => [
            "Running at ECMWF" => "$HSD/Harmonie-mSMS.md",
            "Fast start on cca" => "$HSD/Fast_start_on_cca.md",
            "Quick start local" => "$HSD/QuickStartLocal.md",
        ],
        "Installation" => [
            "Installation" => "$HSD/Installation.md",
            "General Software Requirements" => "$HSD/General.md",
            "Makeup" => "$HSD/Build_with_makeup.md",
            "Stand Alone ODB" => "$HSD/StandaloneOdb.md",
            "Redhat7" => "$HSD/Redhat7Install.md",
            "Centos6" => "$HSD/Centos6Install.md"
        ],
        "The Harmonie Script" => "$HSD/TheHarmonieScript.md",
        "Profiling and Traceback" => "$HSD/DrHook.md",
        "File Formats" => "$HSD/FileFormats.md",
        "Job Scheduling" => [
            "ECFlow" => "$HSD/ECFLOW.md",
            "mXCdp" => "$HSD/scripts/mXCdp.md",
            "mSMS" => "$HSD/scripts/mSMS.md"
        ],
        "Testbed and Mitraillette" => [
            "Testbed" => "$HSD/Evaluation/HarmonieTestbed.md",
            "Mitraillette" => "$HSD/Evaluation/Mitraillette.md"
        ],
        "Directories  and Binaries" => [
            "Directories" => "$HSD/Content.md",
            "Binaries" => "$HSD/Binaries.md"
        ],
        "Configuration" => Configuration,
        "ECFlow tasks" => ECFLOW_TASKS,
        "Other" => [

            "$HSD/Scalability_and_Refactoring.md",
            "$HSD/HarmonieBenchMark.md",
            "$HSD/MFaccess.md",
            "$HSD/Phasing.md"

        ]

]


text:
using JLD2, Base.Test

function writeloop(f, sz)
    for i = 1:sz
        write(f, UInt8(6))
    end
end

# Force growing an MmapIO
name, io = mktemp()
f = JLD2.MmapIO(name, true, true, true)
close(io)
sz = JLD2.MMAP_GROW_SIZE+5
writeloop(f, sz)
truncate(f.f, position(f))
close(f)

@test filesize(name) == sz
f = open(name)
rd = read(f, UInt8, sz)
@test all(rd .== UInt8(6))
close(f)

f = JLD2.MmapIO(name, true, true, true)
write(f, rd)
truncate(f.f, position(f))
close(f)

@test filesize(name) == sz
f = open(name)
rd = read(f, UInt8, sz)
@test all(rd .== UInt8(6))
close(f)

@test JLD2.size_flag(1) === UInt8(0)
@test JLD2.size_flag(256) === UInt8(1)
@test JLD2.size_flag(65536) === UInt8(2)
@test JLD2.size_flag(UInt64(4294967296)) === UInt8(3)

buf = IOBuffer()
JLD2.write_size(buf, 175)
seek(buf, 0)
@test read(buf, UInt8) == 175
seek(buf, 0)
@test JLD2.read_size(buf, UInt8(0)) == 175
seek(buf, 0)

JLD2.write_size(buf, 42580)
seek(buf, 0)
@test read(buf, UInt16) == 42580
seek(buf, 0)
@test JLD2.read_size(buf, UInt8(1)) == 42580
seek(buf, 0)

JLD2.write_size(buf, 3804306107)
seek(buf, 0)
@test read(buf, UInt32) == 3804306107
seek(buf, 0)
@test JLD2.read_size(buf, UInt8(2)) == 3804306107
seek(buf, 0)

JLD2.write_size(buf, 7832227080891617460)
seek(buf, 0)
@test read(buf, UInt64) == 7832227080891617460
seek(buf, 0)
@test JLD2.read_size(buf, UInt8(3)) == 7832227080891617460
seek(buf, 0)

@test JLD2.size_size(1) === 1
@test JLD2.size_size(256) === 2
@test JLD2.size_size(65536) === 4
@test JLD2.size_size(UInt64(4294967296)) === 8


text:
function memset!(_ptr::Ptr, len, _v)
    ptr = Ptr{UInt8}(_ptr)
    v = UInt8(_v)
    for i in Base.OneTo(len)
        unsafe_store!(ptr, v, i)
    end
end

function memcmp(_ptrA::Ptr, _ptrB::Ptr, len)
    ptrA = Ptr{UInt8}(_ptrA)
    ptrB = Ptr{UInt8}(_ptrB)
    for i in Base.OneTo(len)
        (unsafe_load(ptrA, i) != unsafe_load(ptrB, i)) && return false
    end
    return true
end

@testset "Testing Operation of Cache Manager" begin
    Local = CachedArrays.Local
    Remote = CachedArrays.Remote
    PoolType = CachedArrays.PoolType
    getid = CachedArrays.getid

    GC.gc(true)
    minallocation = 12
    manager = CachedArrays.CacheManager(
        CachedArrays.AlignedAllocator(),
        CachedArrays.AlignedAllocator();
        localsize = 2^20,
        remotesize = 2^22,
        minallocation = minallocation,
    )

    # Notes:
    # 2^12 == 4096
    # 2^20 == 1_048_576
    # 2^22 == 4_194_304

    @test length(manager.map) == 0
    @test !haskey(manager.map, 1)
    A = CachedArrays.direct_alloc(PoolType{Local}(), manager, 500_000, getid(manager))
    @test isa(A, CachedArrays.Object)
    @test CachedArrays.pool(A) == Local
    @test getid(CachedArrays.metadata(A)) == 1
    @test haskey(manager.map, 1)

    # Allocation of B - set the data region of B
    B = CachedArrays.direct_alloc(PoolType{Local}(), manager, 500_000, getid(manager))
    @test isa(B, CachedArrays.Object)
    @test CachedArrays.pool(B) == Local
    @test getid(CachedArrays.metadata(B)) == 2
    @test haskey(manager.map, 2)

    # Length of the underlying block should be larger than or equal the requested size.
    len = length(CachedArrays.metadata(B))
    @test len >= 500_000
    @test iszero(mod(sizeof(CachedArrays.metadata(B)), 2^minallocation))
    memset!(pointer(B), length(CachedArrays.metadata(B)), 0xff)
    all_0xff = fill(0xff, len)
    @test memcmp(pointer(B), pointer(all_0xff), len)

    # Local pool should be full, further allocations should fail.
    id = getid(manager)
    @test_throws CachedArrays.AllocationException CachedArrays.direct_alloc(
        PoolType{Local}(),
        manager,
        500_000,
        id,
    )
    @test CachedArrays.check(manager)

    lock(CachedArrays.alloc_lock(manager)) do
        ptr = CachedArrays.unsafe_alloc_direct(PoolType{Local}(), manager, 500_000, id)
        @test ptr === nothing
    end
    @test CachedArrays.check(manager)

    lock(CachedArrays.alloc_lock(manager)) do
        # Okay, now we allocate a remote block for B, link as siblings, and set the
        # primary to the remote.
        block_local = CachedArrays.metadata(B)
        ptr_remote = CachedArrays.unsafe_alloc_direct(
            PoolType{Remote}(),
            manager,
            length(block_local),
            getid(block_local),
        )
        block_remote = CachedArrays.unsafe_block(ptr_remote)
        @test getid(block_remote) == getid(block_local)

        # Test that `copyto!` works correctly.
        @test false == memcmp(
            ptr_remote,
            CachedArrays.datapointer(block_local),
            length(block_local),
        )
        CachedArrays.copyto!(block_remote, block_local, manager)
        @test true == memcmp(
            ptr_remote,
            CachedArrays.datapointer(block_local),
            length(block_local),
        )

        CachedArrays.link!(block_remote, block_local)
        @test CachedArrays.getsibling(block_remote) === block_local
        @test CachedArrays.getsibling(block_local) === block_remote
        old_ptr = CachedArrays.unsafe_setprimary!(manager, block_local, block_remote)
        @test CachedArrays.unsafe_block(old_ptr) === block_local

        # Load the updated block from the backedge to ensure that we did in fact
        # update the backedge.
        map = CachedArrays.getmap(manager)
        new_block = CachedArrays.unsafe_block(unsafe_load(last(map[getid(block_remote)])))
        @test block_remote === new_block
    end
    @test CachedArrays.check(manager)

    # unsafe_free and marking.
    lock(CachedArrays.alloc_lock(manager)) do
        block = CachedArrays.metadata(A)
        CachedArrays.unsafe_free(A)
        @test CachedArrays.isnull(CachedArrays.unsafe_pointer(A))
        @test CachedArrays.candrain(manager)
        @test CachedArrays.isqueued(block)

        # See that `setprimary!` fails.
        ptr_remote = CachedArrays.unsafe_alloc_direct(
            PoolType{Remote}(),
            manager,
            length(block),
            getid(block),
        )
        block_remote = CachedArrays.unsafe_block(ptr_remote)
        CachedArrays.link!(block, block_remote)
        @test CachedArrays.unsafe_setprimary!(manager, block, block_remote) === nothing

        # Now, if we drain the manager, both the local and remote blocks should be
        # cleaned because we linked them.
        cleaned = CachedArrays.unsafe_cleanup!(manager, CachedArrays.getid(block))
        @test cleaned == true
        @test CachedArrays.isfree(block)
        @test CachedArrays.isfree(block_remote)
    end
    @test CachedArrays.check(manager)
end

@testset "Testing Manager Cleanup 1" begin
    GC.gc(true)
    @test length(CachedArrays.GlobalManagers) == 1
    manager = only(CachedArrays.GlobalManagers)
    @test CachedArrays.cangc(manager)
    @test CachedArrays.check(manager)
end

@testset "Testing Manager Cleanup 2" begin
    v = CachedArrays.gc_managers()
    @test v == 1
end


text:
using LinearAlgebra, DiffEqOperators, Random, Test
################################################################################
# Test 2d extension
################################################################################

#Create Array
n = 8
m = 15
A = rand(n,m)

#Create atomic BC
q1 = RobinBC((1.0, 2.0, 3.0), (0.0, -1.0, 2.0), 0.1, 4.0)
q2 = PeriodicBC(Float64)

BCx = vcat(fill(q1, div(m,2)), fill(q2, m-div(m,2)))  #The size of BCx has to be all size components *except* for x
BCy = vcat(fill(q1, div(n,2)), fill(q2, n-div(n,2)))

Qx = MultiDimBC{1}(BCx)
Qy = MultiDimBC{2}(BCy)

Ax = Qx*A
Ay = Qy*A

@test size(Ax)[1] == size(A)[1]+2
@test size(Ay)[2] == size(A)[2]+2

for j in 1:m
    @test Ax[:, j]  == Array(BCx[j]*A[:, j])
end
for i in 1:n
    @test Ay[i,:] == Array(BCy[i]*A[i,:])
end

################################################################################
# Test 3d extension
################################################################################

#Create Array
n = 8
m = 11
o = 12
A = rand(n,m, o)

#Create atomic BC
q1 = RobinBC((1.0, 2.0, 3.0), (0.0, -1.0, 2.0), 0.1, 4.0)
q2 = PeriodicBC(Float64)

BCx = vcat(fill(q1, (div(m,2), o)), fill(q2, (m-div(m,2), o)))  #The size of BCx has to be all size components *except* for x
BCy = vcat(fill(q1, (div(n,2), o)), fill(q2, (n-div(n,2), o)))
BCz = fill(Dirichlet0BC(Float64), (n,m))

Qx = MultiDimBC{1}(BCx)
Qy = MultiDimBC{2}(BCy)
Qz = Dirichlet0BC{3}(Float64, size(A)) #Test the other constructor

Q1 = (Qx+Qy+Qz)
Q2 = compose(Qx,Qy,Qz) #test addition combinations
@test_broken Q1 == Q2 #This fails
@test all(Q1.BCs .== Q2.BCs) # but this passes so it does actually work
@test_broken Qtmp = Qx + Qz
@test_skip Qz+Qx+Qy == Qy+Qtmp
Ax = Qx*A
Ay = Qy*A
Az = Qz*A

@test size(Ax)[1] == size(A)[1]+2
@test size(Ay)[2] == size(A)[2]+2
@test size(Az)[3] == size(A)[3]+2
for j in 1:m, k in 1:o
    @test Ax[:, j, k] == Array(BCx[j, k]*A[:, j, k])
end
for i in 1:n, k in 1:o
    @test Ay[i, :, k] == Array(BCy[i, k]*A[i, :, k])
end
for i in 1:n, j in 1:m
    @test Az[i, j, :] == Array(BCz[i, j]*A[i, j, :])
end

#test compositions to higher dimension
for N in 2:7
    sizes = rand(4:7, N)
    A = rand(sizes...)

    Q1_N = Neumann0BC(Float64, Tuple(ones(N)), 3.0, size(A))

    Q = compose(Q1_N...)

    A1_N = Q1_N.*fill(A, N)

    A_extended = Q*A
    @test Array(A_extended) == Array(compose(A1_N...))
end


text:
@testset "Construction" begin

end

@testset "Conversion" begin

end

@testset "Evaluation" begin

end

@testset "Arithmetic" begin

end

@testset "Elementwise operations" begin

end

@testset "Divrem" begin

end

@testset "Comparisons" begin

end

@testset "Derivatives and ntegrals" begin

end

@testset "roots" begin

end

@testset "fitting" begin

end


text:
include("../roblib.jl")

function mark_ladder(r)
    num_steps = through_rectangles_into_angle(r, (Down, Left))

    cells_in_raw = mark_and_enumerate_if_possible!(r, Right)
    cells_to_mark = cells_in_raw - 1

    while (isborder(r, Up) == false) || (isborder(r, Left) == false)
        movements_if_possible!(r, Left)
        if (isborder(r, Up) == false)
            move!(r, Up)
        end
        cells_to_not_mark = cells_in_raw - get_num_movements_if_possible!(r, Right)
        movements_if_possible!(r, Left, cells_in_raw)
        cells_remained = move_and_check(r, Right, cells_to_mark)
        if ((cells_remained - cells_to_not_mark) > 0)
            putmarkers_if_possible!(r, Right, cells_remained - cells_to_not_mark)
        end
        cells_to_mark -= 1
    end

    through_rectangles_into_angle(r, (Down, Left))
    movements_if_possible!(r, (Up, Right), reverse(num_steps))
end

function move_and_check(r, side, count)
    putmarker!(r)
    while (count!=0 && isborder(r, side) == false)
        move!(r,side)
        putmarker!(r)
        count -= 1
    end

    return count
end


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

mutable struct RunCommandParameterDefinition <: SwaggerModel
    name::Any # spec type: Union{ Nothing, String } # spec name: name
    type::Any # spec type: Union{ Nothing, String } # spec name: type
    defaultValue::Any # spec type: Union{ Nothing, String } # spec name: defaultValue
    required::Any # spec type: Union{ Nothing, Bool } # spec name: required

    function RunCommandParameterDefinition(;name=nothing, type=nothing, defaultValue=nothing, required=false)
        o = new()
        validate_property(RunCommandParameterDefinition, Symbol("name"), name)
        setfield!(o, Symbol("name"), name)
        validate_property(RunCommandParameterDefinition, Symbol("type"), type)
        setfield!(o, Symbol("type"), type)
        validate_property(RunCommandParameterDefinition, Symbol("defaultValue"), defaultValue)
        setfield!(o, Symbol("defaultValue"), defaultValue)
        validate_property(RunCommandParameterDefinition, Symbol("required"), required)
        setfield!(o, Symbol("required"), required)
        o
    end
end # type RunCommandParameterDefinition

const _property_map_RunCommandParameterDefinition = Dict{Symbol,Symbol}(Symbol("name")=>Symbol("name"), Symbol("type")=>Symbol("type"), Symbol("defaultValue")=>Symbol("defaultValue"), Symbol("required")=>Symbol("required"))
const _property_types_RunCommandParameterDefinition = Dict{Symbol,String}(Symbol("name")=>"String", Symbol("type")=>"String", Symbol("defaultValue")=>"String", Symbol("required")=>"Bool")
Base.propertynames(::Type{ RunCommandParameterDefinition }) = collect(keys(_property_map_RunCommandParameterDefinition))
Swagger.property_type(::Type{ RunCommandParameterDefinition }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_RunCommandParameterDefinition[name]))}
Swagger.field_name(::Type{ RunCommandParameterDefinition }, property_name::Symbol) =  _property_map_RunCommandParameterDefinition[property_name]

function check_required(o::RunCommandParameterDefinition)
    (getproperty(o, Symbol("name")) === nothing) && (return false)
    (getproperty(o, Symbol("type")) === nothing) && (return false)
    true
end

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


text:
function mdl_Zload!(
    voltage_r::T,
    voltage_i::T,
    current_r::AbstractArray{T},
    current_i::AbstractArray{T},
    static_device::PSY.PowerLoad,
) where {T <: Real}
    #Load squared voltage magnitude at steady state
    Vmag_sq = PSY.get_magnitude(PSY.get_bus(static_device))^2
    #Vmag_sq = 1.0

    #Load device parameters
    P = PSY.get_active_power(static_device)
    Q = PSY.get_reactive_power(static_device)

    #Compute impedance and load: Z = |V|^2/conj(S)
    #Z_load =  Vmag_sq/(P-Q*1im) #in pu
    #I = -(voltage_r[1] + voltage_i[1]*1im)/Z_load #in pu flowing out
    #current_r[1] += real(I)
    #current_i[1] += imag(I)

    #Update current
    #This model creates an equivalent RL/RC circuit based on steady state voltage
    current_r[1] += -(1.0 / Vmag_sq) * (voltage_r * P + voltage_i * Q) #in system pu flowing out
    current_i[1] += -(1.0 / Vmag_sq) * (voltage_i * P - voltage_r * Q) #in system pu flowing out

    return
end


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

using Markdown
using InteractiveUtils

# ╔═╡ 3f352894-81a6-11eb-3fd7-0333460c112f
using DifferentialEquations, StaticArrays

# ╔═╡ 9838c218-83c0-11eb-31fa-67cd1c50ff16
using ForwardDiff: derivative

# ╔═╡ e88be666-81a7-11eb-12ed-67300d719779
using Plots, LaTeXStrings

# ╔═╡ 36337070-83c4-11eb-11ce-9ff923e2e1b7
using HTTP, CSV

# ╔═╡ e5301ace-83c7-11eb-00cc-5fcafe9dde00
using Interpolations

# ╔═╡ 708f50ea-81a6-11eb-1cd8-477238b61afe
md"# Solve TOV equation for neutron star
Md Arif Shaikh, Postdoc, ICTS-TIFR, March 10 2021
"

# ╔═╡ 929b5bca-81a6-11eb-3793-dffff4210e82
md"## Newtonian Case:
First we solve for Newtonian stellar structure equations"

# ╔═╡ a9558cf2-81a6-11eb-0afc-9b1133f43885
function StellarStrunctureNewt(u, params, r)
	m, ρ = u
	G, K, ρ0 = params
	pofρ(ρ) = K * abs(ρ)^(5/3) * (1 - (ρ0/abs(ρ))^(1/3))
	dpdρ = derivative(pofρ, ρ)
	# dpdρ = K * abs(ρ)^(2/3)*(5 - 4 * (ρ0/abs(ρ))^(1/3))/3
	dpdr = - G * m * ρ / r^2
	dρdr = dpdr/dpdρ
	dmdr = 4 * pi * r^2 * ρ
	@SVector [dmdr, dρdr]
end

# ╔═╡ b3e63b26-81a6-11eb-3536-cfc5c3cd70e2
md"Define the constants"

# ╔═╡ ba55efcc-81a6-11eb-3d16-43035b98aa9f
begin
	G = 6.6743 * 1e-8
	μe = 2
	Z = 6
	ρ0 = 0.4 * Z^2 * μe
	K = 1e13 * μe^(-5/3)
end

# ╔═╡ cf0ca384-81a6-11eb-3be4-dbd44e9b276c
md"Define the callback condition"

# ╔═╡ dbd5633a-81a6-11eb-2c33-e1555b27a642
condition(u, r, integrator) = u[2] - ρ0

# ╔═╡ f2042eca-81a6-11eb-12c3-abe61c337809
md" Stop the integration when we hit the surface of the star, i.e  when $\rho$ becomes equal to $\rho_0$"

# ╔═╡ e6737ce6-81a6-11eb-0e25-d3304d934ea7
affect!(integrator) = terminate!(integrator)

# ╔═╡ 06468d24-81a7-11eb-150f-1b0190fb30cc
cb = ContinuousCallback(condition, affect!)

# ╔═╡ 0f690b34-81a7-11eb-2a6a-b5e92506ba43
md"Now we define Newtonian ODE problem."

# ╔═╡ 73a5c95c-81a7-11eb-3011-2f51c93c10b2
begin
	u0 = @SVector [0, 1e8]  # initial values
	rspan = [1e-10, 1e10]  # range of integration
	params = [G, K, ρ0]  # parameters in the ODE
	probNewt = ODEProblem(StellarStrunctureNewt, u0, rspan, params)
end

# ╔═╡ b13fef66-81a7-11eb-38ff-491f793cde05
md"Solution"

# ╔═╡ b7f2a512-81a7-11eb-1e24-ddd9b8bd0f90
solNewt = solve(probNewt, callback=cb)

# ╔═╡ d50944bc-81a7-11eb-2f41-c1de1d753953
md"Plot the mass vs radius"

# ╔═╡ fced62ba-81a7-11eb-3e7f-4b95cf280ef5
begin
	Msun = 1.989 * 1e33
	# the solution of DifferentialEquations solver returns an interpolation of the solution. So we can choose approriate time axis and plot a smooth curve.
	trange = LinRange(solNewt.t[1], solNewt.t[end], 1000)
	plot(trange ./ 1e5, solNewt(trange)[1, :] ./Msun, xlabel=L"\textrm{radius [km]}", 			ylabel=L"$\textrm{Mass }[M_\odot]$", label=L"$\rho_c = 10^{8}$", 					legend=:bottomright)
end

# ╔═╡ 3001faa6-81a8-11eb-2c91-43731d360e7c
md"Now we solve for Mass of the star vs radius of the star for different values of the central density."

# ╔═╡ 534fd89a-81a8-11eb-156a-e3cf27733d96
begin
	ρcs = 10 .^ LinRange(6, 9, 25)
	MassesNewt = zeros(size(ρcs))
	radiiNewt = zeros(size(ρcs))
end

# ╔═╡ 64eea234-81a8-11eb-032a-15ffc27dcac3
md"We now solve the Newtonian structure equation for different central densities.
We want to save only the events and therefore we modify the callback function slightly."

# ╔═╡ 7932ec20-81a8-11eb-2a03-c70394f794d5
cbMassRadiiNewt = ContinuousCallback(condition, affect!, save_positions=(false,true))

# ╔═╡ 8e840ac6-81a8-11eb-32aa-2b7dfb553f1d
for (idx, ρc) in enumerate(ρcs)
	u0 = @SVector [0.0, ρc]
	probMassRadiiNewt = ODEProblem(StellarStrunctureNewt, u0, rspan, params)
	# since we want to save only the events we add few more otions in our solver
	solMassRadiiNewt = solve(probMassRadiiNewt, callback=cbMassRadiiNewt, 	save_everystep=false, save_start=false, save_end=false)
	# Save the events
	MassesNewt[idx] = solMassRadiiNewt.u[end][1]
	radiiNewt[idx] = solMassRadiiNewt.t[end]
end

# ╔═╡ f0c763b6-81a8-11eb-32ba-0304d36c7fe6
scatter(radiiNewt ./ 1e5, MassesNewt ./ Msun, zcolor=log10.(ρcs),
	xlabel=L"$\textrm{radius [km]}$",
	ylabel=L"$\textrm{Mass }[M_\odot]$",
	label=L"$\textrm{central density}$")

# ╔═╡ 36e57fc4-81a9-11eb-21fe-654f42968a86
md"## Solve TOV

We solve solve the TOV equation in this section using the same equation of state as the earlier section."

# ╔═╡ fd2bffbe-81a9-11eb-294b-3f09a5cb31aa
md"Write down the TOV equations with $G$ and $c$ restored.

$$m\to Gm/c^2$$

$$\rho\to G\rho/c^2$$

$$p\to Gp/c^4$$

So the TOV equations become

$$\frac{dp}{dr} = -\frac{Gm\rho}{r^2}\left[1 + \frac{1}{c^2}\frac{p}{\rho}\right]\left[1 + \frac{1}{c^2}\frac{4\pi r^3 p}{m}\right]\left[1 - \frac{2Gm}{c^2r}\right]^{-1}$$"

# ╔═╡ 1b29c084-81aa-11eb-1788-4946c721fc33
function tov(u, params, r)
    m, ρ = u
	G, C, K, ρ0 = params
    pofρ(ρ) = K * abs(ρ)^(5/3) * (1 - (ρ0/abs(ρ))^(1/3))
    dmdr = 4 * pi * r^2 * ρ
    dpdr = - (G * m * ρ / r^2) * (1 + (pofρ(ρ)/(C^2 * ρ))) * (1 + (4* pi * r^3 * pofρ(ρ)/(C^2 * m))) * (1 - (2 * G * m/(C^2 * r)))^(-1)
    # dpdρ = K * abs(ρ)^(2/3) * (5 - 4 * (ρ0/abs(ρ))^(1/3)) / 3
	dpdρ = derivative(pofρ, ρ)
    dρdr = dpdr/dpdρ

	@SVector [dmdr, dρdr]
end

# ╔═╡ eccb8058-81aa-11eb-3bd5-ffb497b470f4
md"First solve for a single central density and compare with the Newtonian solution"

# ╔═╡ 3728e3e8-81ab-11eb-0ffb-f1a9360d6ee9
begin
	c = 1e10
	u0TOV = @SVector [1e-10, u0[2]]
	paramsTOV = [G, c, K, ρ0]
	ProbTOV = ODEProblem(tov, u0TOV, rspan, paramsTOV)
end

# ╔═╡ aed24844-81ab-11eb-0a02-4f3c1af140ff
solveTOV = solve(ProbTOV, callback=cb)

# ╔═╡ 712b33ba-81ac-11eb-2b07-eb59d01b946d
md"Plot both Newtonian and TOV solution"

# ╔═╡ 80b123b4-81ac-11eb-2c41-a3e4826ef900
begin
	p = plot(xlabel=L"$r \textrm{ [km]}$", ylabel=L"$m [M_\odot]$")
	plot!(p, trange ./ 1e5, solNewt(trange)[1, :] ./Msun, label=L"$\textrm{Newtonian}$")
	trangeTOV = LinRange(solveTOV.t[1], solveTOV.t[end], 1000)
	plot!(p, trangeTOV ./ 1e5, solveTOV(trangeTOV)[1, :] ./Msun, label=L"$\textrm{TOV}$", legend=:bottomright)
end

# ╔═╡ c0a059ba-81ad-11eb-0d10-c510ec43175e
md"As before we can vary the central density and plot mass vs radius"

# ╔═╡ d04ac76a-81ad-11eb-36b2-df50ae9fd5e1
begin
	MassesTOV = zeros(size(ρcs))
	radiiTOV = zeros(size(ρcs))
end

# ╔═╡ 049d4e8e-81ae-11eb-3563-e986cf2198cd
for (idx, ρc) in enumerate(ρcs)
	u0 = @SVector [1e-10, ρc]
	probMassRadiiTOV = ODEProblem(tov, u0, rspan, paramsTOV)
	# since we want to save only the events we add few more otions in our solver
	solMassRadiiTOV = solve(probMassRadiiTOV, callback=cbMassRadiiNewt, 	save_everystep=false, save_start=false, save_end=false)
	# Save the events
	MassesTOV[idx] = solMassRadiiTOV.u[end][1]
	radiiTOV[idx] = solMassRadiiTOV.t[end]
end

# ╔═╡ 4cfe94a8-81ae-11eb-2099-39022c741d8e
md"Plot mass vs radii for both Newtonian and TOV"

# ╔═╡ 5dce1b8c-81ae-11eb-17fc-dbd92026d51b
begin
	pTOVvsNewt = plot(xlabel=L"$r \textrm{ [km]}$", ylabel=L"$m [M_\odot]$")
	scatter!(pTOVvsNewt, radiiNewt ./ 1e5, MassesNewt ./ Msun, zcolor=log10.(ρcs), label=L"$\textrm{Newtonian}$")
	scatter!(pTOVvsNewt, radiiTOV ./ 1e5, MassesTOV ./ Msun, zcolor=log10.(ρcs), markershape =:square, label=L"$\textrm{TOV}$")
end

# ╔═╡ Cell order:
# ╟─708f50ea-81a6-11eb-1cd8-477238b61afe
# ╠═3f352894-81a6-11eb-3fd7-0333460c112f
# ╠═9838c218-83c0-11eb-31fa-67cd1c50ff16
# ╟─929b5bca-81a6-11eb-3793-dffff4210e82
# ╠═a9558cf2-81a6-11eb-0afc-9b1133f43885
# ╟─b3e63b26-81a6-11eb-3536-cfc5c3cd70e2
# ╠═ba55efcc-81a6-11eb-3d16-43035b98aa9f
# ╟─cf0ca384-81a6-11eb-3be4-dbd44e9b276c
# ╠═dbd5633a-81a6-11eb-2c33-e1555b27a642
# ╟─f2042eca-81a6-11eb-12c3-abe61c337809
# ╠═e6737ce6-81a6-11eb-0e25-d3304d934ea7
# ╠═06468d24-81a7-11eb-150f-1b0190fb30cc
# ╟─0f690b34-81a7-11eb-2a6a-b5e92506ba43
# ╠═73a5c95c-81a7-11eb-3011-2f51c93c10b2
# ╟─b13fef66-81a7-11eb-38ff-491f793cde05
# ╠═b7f2a512-81a7-11eb-1e24-ddd9b8bd0f90
# ╟─d50944bc-81a7-11eb-2f41-c1de1d753953
# ╠═e88be666-81a7-11eb-12ed-67300d719779
# ╠═fced62ba-81a7-11eb-3e7f-4b95cf280ef5
# ╟─3001faa6-81a8-11eb-2c91-43731d360e7c
# ╠═534fd89a-81a8-11eb-156a-e3cf27733d96
# ╟─64eea234-81a8-11eb-032a-15ffc27dcac3
# ╠═7932ec20-81a8-11eb-2a03-c70394f794d5
# ╠═8e840ac6-81a8-11eb-32aa-2b7dfb553f1d
# ╠═f0c763b6-81a8-11eb-32ba-0304d36c7fe6
# ╟─36e57fc4-81a9-11eb-21fe-654f42968a86
# ╟─fd2bffbe-81a9-11eb-294b-3f09a5cb31aa
# ╠═1b29c084-81aa-11eb-1788-4946c721fc33
# ╟─eccb8058-81aa-11eb-3bd5-ffb497b470f4
# ╠═3728e3e8-81ab-11eb-0ffb-f1a9360d6ee9
# ╠═aed24844-81ab-11eb-0a02-4f3c1af140ff
# ╟─712b33ba-81ac-11eb-2b07-eb59d01b946d
# ╠═80b123b4-81ac-11eb-2c41-a3e4826ef900
# ╟─c0a059ba-81ad-11eb-0d10-c510ec43175e
# ╟─d04ac76a-81ad-11eb-36b2-df50ae9fd5e1
# ╠═049d4e8e-81ae-11eb-3563-e986cf2198cd
# ╟─4cfe94a8-81ae-11eb-2099-39022c741d8e
# ╠═5dce1b8c-81ae-11eb-17fc-dbd92026d51b
# ╟─18e64b4c-81af-11eb-0757-47bd91603d06
# ╟─1a190bd0-83c4-11eb-278f-81628d9f818f
# ╠═36337070-83c4-11eb-11ce-9ff923e2e1b7
# ╠═d5358d76-83c6-11eb-1cbb-5b8da156e2ce
# ╟─0989088c-83c7-11eb-1eff-b532874bb091
# ╠═e5301ace-83c7-11eb-00cc-5fcafe9dde00
# ╠═c7147610-83c8-11eb-0ac6-730563f1a04a
# ╠═7e88240e-83c9-11eb-2248-15d4da1ac62d
# ╠═6b709b54-83c8-11eb-13c1-0da2153e0aac
# ╠═8a3d09b8-83c9-11eb-2110-83bd65528e6a
# ╠═b4666814-83cc-11eb-27c1-2714f6d02638
# ╟─461500c2-83c9-11eb-3bce-9f3af623761e
# ╠═b1faf5b2-83c8-11eb-19d2-bf79e95dbe55
# ╠═ac04e2e4-83c9-11eb-0848-a5509d39b134
# ╠═c3757474-83c9-11eb-0084-e75e47fca316
# ╟─f4f39308-83c9-11eb-2cfd-71335de54d22
# ╠═1a787080-83ca-11eb-36aa-67de913dcb52
# ╠═dae3f700-83c9-11eb-2d90-454aa644f10e
# ╠═9bdc6b36-83ca-11eb-11fc-3336722450b5
# ╠═ed06ccc8-83ce-11eb-0d48-1b0321a68567
# ╠═f5dd09d2-83ce-11eb-0e08-cda8fb30830a
# ╠═b7cd8fa2-8729-11eb-1c6a-c5742b33b554
# ╠═b8a145b6-83ca-11eb-0418-1d91b983a3e5


text:
### Data structures ###
abstract type OptData <: InputData end

"""
     OptModelCEP
-model::JuMP.Model
-info::Array{String}
-set::Dict{String,Array}
"""
struct OptModelCEP
  model::JuMP.Model
  info::Array{String}
  set::Dict{String,Array}
end

"""
     OptVariable
-`data::Array` - includes the optimization variable output in  form of an array
-`axes_names::Array{String,1}`` - includes the names of the different axes and is equivalent to the sets in the optimization formulation
-`axes::Tuple` - includes the values of the different axes of the optimization variables
-`type::String` - defines the type of the variable being cv - cost variable - dv -design variable - ov - operating variable - sv - slack variable
"""
struct OptVariable{T,N,Ax,L<:NTuple{N,Dict}} <: AbstractArray{T,N}
    data::Array{T,N}
    axes::Ax
    lookup::L
    axes_names::Array{String,1}
    type::String
end

"""
      OptResult{status::Symbol,objective::Float64,variables::Dict{String,Any},sets::Dict{String,Array},opt_config::Dict{String,Any},opt_info::Dict{String,Any}}
- `status`: Symbol about the solution status of the model in normal cases `:OPTIMAL`
- `objective`: Value of the objective function
- `variables`: Dictionary with each OptVariable as an entry
- `sets`: Dictionary with each set as an entry
- `opt_config`: The configuration of the model setup - for more detail see tye `run_opt` documentation that sets the `opt_config` up
- `opt_info`: Holds information about the model. E.g. `opt_info["model"]` contains the exact equations used in the model.
"""
struct OptResult
 status::Symbol
 objective::Float64
 variables::Dict{String,Any}
 sets::Dict{String,Array}
 opt_config::Dict{String,Any}
 opt_info::Dict{String,Any}
end

"""
     OptDataCEP{region::String, costs::OptVariable, techs::OptVariable, nodes::OptVariable, lines::OptVariabl} <: OptData
- `region::String`:          name of state or region data belongs to
- `costs::OptVariable`:    costs[tech,node,year,account,impact] - Number
- `techs::OptVariable`:    techs[tech] - OptDataCEPTech
- `nodes::OptVariable`:    nodes[tech, node] - OptDataCEPNode
- `lines::OptVarible`:     lines[tech, line] - OptDataCEPLine
instead of USD you can also use your favorite currency like EUR
"""
struct OptDataCEP <: OptData
   region::String
   costs::OptVariable
   techs::OptVariable
   nodes::OptVariable
   lines::OptVariable
end

#struct LatLon() adapted from Package Geodesy.jl: Copyright (c) 2014-2016: Ted Steiner, Sean Garborg, Yeesian Ng, Andy Ferris, Andrew Smith

#Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:

#The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.
"""
    LatLon(lat, lon)
    LatLon(lat = ϕ, lon = Θ)
Latitude and longitude co-ordinates. *Note:* assumes degrees not radians
"""
struct LatLon{T <: Number}
    lat::T
    lon::T
end
LatLon(lat::Number, lon::Number) = LatLon(promote(lat, lon)...)
LatLon(;lat=NaN,lon=NaN) = LatLon(lat,lon) # Constructor that is independent of storage order
Base.show(io::IO, ll::LatLon) = print(io, "LatLon(lat=$(ll.lat)°, lon=$(ll.lon)°)")
Base.isapprox(ll1::LatLon, ll2::LatLon; atol = 1e-6, kwargs...) = isapprox(ll1.lat, ll2.lat; atol = 180*atol/6.371e6, kwargs...) & isapprox(ll1.lon, ll2.lon; atol = 180*atol/6.371e6, kwargs...) # atol in metres (1μm)

"""
     OptDataCEPNode{name::String,value::Number,lat::Number,lon::Number} <: OptData
- `name`
- `power_ex` existing capacity [MW or MWh (tech_e)]
- `power_lim` capacity limit [MW or MWh (tech_e)]
- `region`
- `latlon` hold geolocation information [°,°]
"""
struct OptDataCEPNode <: OptData
  name::String
  power_ex::Number
  power_lim::Number
  region::String
  latlon::LatLon
end

"""
     OptDataCEPLine{name::String,node_start::String,node_end::String,reactance::Number,resistance::Number,power::Number,circuits::Int,voltage::Number,length::Number} <: OptData
- `name`
- `node_start` Node where line starts
- `node_end` Node where line ends
- `reactance`
- `resistance` [Ω]
- `power_ex`: existing power limit [MW]
- `power_lim`: limit power limit [MW]
- `circuits` [-]
- `voltage` [V]
- `length` [km]
- `eff` [-]
"""
struct OptDataCEPLine <: OptData
  name::String
  node_start::String
  node_end::String
  reactance::Number
  resistance::Number
  power_ex::Number
  power_lim::Number
  circuits::Int
  voltage::Number
  length::Number
  eff::Number
end

"""
     OptDataCEPTech{name::String,categ::String,sector::String,eff::Number,time_series::String,lifetime::Number,financial_lifetime::Number,discount_rate::Number, annuityfactor::Number} <: OptData
- `name`
- `categ`: the category of this technology (is it storage, transmission or generation)
- `sector`: sector of the technology (electricity or heat)
- `eff`: efficiency of this technologies conversion [-]
- `time_series`: time_series name for availability
- `lifetime`: product lifetime [a]
- `financial_lifetime`: financial time to break even [a]
- `discount_rate`: discount rate for technology [a]
- `annuityfactor`: annuity factor, important for cap-costs [-]
"""
struct OptDataCEPTech <: OptData
  name::String
  categ::String
  sector::String
  eff::Number
  time_series::String
  lifetime::Number
  financial_lifetime::Number
  discount_rate::Number
  annuityfactor::Number
end

"""
  is_in(k::Symbol,table::DataFrame,alt_value::Any)
is Symbol `k` in `table`? Lookup value if true, return `alt_value` if false
"""
function is_in(k::Symbol,table::DataFrame,alt_value::Any)
  if k in names(table)
    return table[k][1]
  else
    @warn "$k not provided in $(repr(table))"
    return alt_value
  end
end

"""
     Scenario{descriptor::String,clust_res::ClustResult,opt_res::OptResult}
-`descriptor::String`
-`clust_res::ClustResult`
-`opt_res::OptResult`
"""
struct Scenario
 descriptor::String
 clust_res::ClustResult
 opt_res::OptResult
end


text:
module FerriteMeshParser
using Ferrite
# Convenience when debugging
const DEBUG_PARSE = false

# Mesh formats
struct AutomaticMeshFormat end
struct AbaqusMeshFormat end

# Exceptions
struct UndetectableMeshFormatError <: Exception
    filename::String
end
Base.showerror(io::IO, e::UndetectableMeshFormatError) = println(io, "Couldn't automatically detect mesh format in $(e.filename)")

struct UnsupportedElementType <: Exception
    elementtype::String
end
Base.showerror(io::IO, e::UnsupportedElementType) = println(io, "The element type \"$(e.elementtype)\" is not supported or given in user_elements")

struct InvalidFileContent <: Exception
    msg::String
end
Base.showerror(io::IO, e::InvalidFileContent) = println(io, e.msg)

include("rawmesh.jl")
include("elements.jl")
include("reading_utils.jl")
include("abaqusreader.jl")
include("gridcreator.jl")

"""
    function get_ferrite_grid(filename; meshformat=AutomaticMeshFormat(), user_elements=Dict{String,DataType}(), generate_facesets=true)

    Greate a `Ferrite.Grid` by reading in the file specified by `filename`.

    ## Optional arguments (default)
    * `meshformat` (`FerriteMeshParser.AutomaticMeshFormat()`): Which format the mesh
      is given in, normally automatically detected by the file extension
    * `user_elements` (`Dict{String,DataType}()`): Used to add extra elements not supported,
      might require separate cell constructor.
    * `generate_facesets` (`true`): Should facesets be detected automatically from all nodesets?

"""
function get_ferrite_grid(filename; meshformat=AutomaticMeshFormat(), user_elements::Dict{String,DataType}=Dict{String,DataType}(), generate_facesets::Bool=true)
    detected_format = detect_mesh_format(filename, meshformat)
    mesh = read_mesh(filename, detected_format)
    checkmesh(mesh)
    grid = create_grid(mesh, detected_format, user_elements)
    generate_facesets && generate_facesets!(grid)
    return grid
end

"""
    create_faceset(grid::Ferrite.AbstractGrid, nodeset::Set{Int}, cellset::Union{Nothing,Set{Int}}=nothing)

Find the faces in the grid for which all nodes are in `nodeset`. Return them as a `Set{FaceIndex}`.
A `cellset` can be given to only look only for faces amongst those cells to speed up the computation.
Otherwise the search is over all cells.

This function is normally only required if calling `get_ferrite_grid` with `generate_facesets=false`.
The created `faceset` can be added to the grid as `merge!(getfacesets(grid), Dict("facesetkey" => faceset`
"""
function create_faceset(grid::Ferrite.AbstractGrid, nodeset::Set{Int}, ::Nothing=nothing)
    faceset = Set{FaceIndex}()
    for (cellid, cell) in enumerate(getcells(grid))
        if any(map(n-> n ∈ nodeset, cell.nodes))
            # check actual faces
            for (faceid, face) in enumerate(Ferrite.faces(cell))
                if all(map(n -> n ∈ nodeset, face))
                    push!(faceset, FaceIndex(cellid, faceid))
                end
            end
        end
    end
    return faceset
end

function create_faceset(grid::Ferrite.AbstractGrid, nodeset::Set{Int}, cellset::Set{Int})
    faceset = Set{FaceIndex}()
    cellrange = collect(cellset)
    for (i, cell) in enumerate(getcells(grid)[cellrange])
        cellid = cellrange[i]
        for (faceid, face) in enumerate(Ferrite.faces(cell))
            if all(map(n -> n ∈ nodeset, face))
                push!(faceset, FaceIndex(cellid, faceid))
            end
        end
    end
    return faceset
end

detect_mesh_format(_, meshformat) = meshformat
function detect_mesh_format(filename, ::AutomaticMeshFormat)
    if endswith(filename, ".inp")
        return AbaqusMeshFormat()
    else
        throw(UndetectableMeshFormatError(filename))
    end
end

export get_ferrite_grid, create_faceset

end


text:
module InvestorsExchange

export
    # download functions
    download_file, query_hist_filenames,
    # PCAP-parsing functions
    assemble_trade_report_frame, read_trade_report_messages

include("download.jl")
include("read_pcap.jl")

end


text:
#] add OpenStreetMapX DelimitedFiles CSV LightGraphs DataFrames Statistics Distributions HypothesisTests SparseArrays Serialization Parameters LinearAlgebra StatsBase PyCall Conda
#] precompile

#println("Using $(Threads.nthreads()) threads!")
#using Revise

using Pkg
Pkg.activate(".")
using SignalBroadcastingSim
using Random
using OpenStreetMapX
using DelimitedFiles, CSV
using LightGraphs
using DataFrames
using Statistics
using Distributions, HypothesisTests

using Distributed
addprocs(Sys.CPU_THREADS)
println("Number of workers $(nworkers())")

@everywhere using Pkg
@everywhere Pkg.activate(".")
@everywhere using SignalBroadcastingSim
@everywhere using Random
@everywhere using OpenStreetMapX
@everywhere using DelimitedFiles, CSV
@everywhere using LightGraphs
@everywhere using DataFrames
@everywhere using Statistics
@everywhere using Distributions, HypothesisTests

m = OpenStreetMapX.get_map_data("maps", "torontoF.osm")
#discretize_m = 25
#p = ModelParams(discretize_m=discretize_m)
#s = Simulation(p,m)

dfstepstats = DataFrame()
dfendstats = DataFrame()
for jump_infect in [false, true]
    for n_agents in vcat(10:30:100, 150:50:200, 250:150:1000, 2000:1000:10000, 12000:2000:20000)
        for discretize_m in [50]
            p = ModelParams(n_agents=n_agents, discretize_m=discretize_m, jump_infect=jump_infect)
            s = Simulation(p,m; store_map = false)
            dump(p)
            simres = simulate_all!(s; reps=40*96, max_idle=50000)

			CSV.write("zombiecar2_sweep_res5_step_$(n_agents)_$(jump_infect)_$(discretize_m).csv", simres.stepstats)
			CSV.write("zombiecar2_sweep_res5_end_$(n_agents)_$(jump_infect)_$(discretize_m).csv", simres.endstats)
			append!(dfstepstats, simres.stepstats)
            append!(dfendstats, simres.endstats)
        end
    end
    CSV.write("zombiecar2_sweep_res5_step$(jump_infect).csv", dfstepstats)
    CSV.write("zombiecar2_sweep_res5_end$(jump_infect).csv", dfendstats)
end

CSV.write("zombiecar2_sweep_res4_step.csv", dfstepstats)
CSV.write("zombiecar2_sweep_res4_end.csv", dfendstats)

function doSim(m;jump_infect=true)
    Random.seed!(0);
    p = ModelParams(n_agents=500, discretize_m=75.0, jump_infect=jump_infect)

    s = Simulation(p,m)

    plot_sim(s, "SimStart$jump_infect.html")

    for i in 1:10
        step!(s)
    end
    plot_sim(s, "SimStep$(jump_infect)10.html")
    for i in 11:20
        step!(s)
    end
    plot_sim(s, "SimStep$(jump_infect)20.html")

    for i in 21:100
        step!(s)
    end
    plot_sim(s, "SimStep$(jump_infect)100.html")
    for i in 101:200
        step!(s)
    end
    plot_sim(s, "SimStep$(jump_infect)200.html")
    for i in 201:500
        step!(s)
    end
    plot_sim(s, "SimStep$(jump_infect)500.html")
end

#doSim(m, jump_infect=true)
#doSim(m, jump_infect=false)


text:
# Apply fitted trees or forests

function apply_tree(leaf::Leaf, x::DataArray)
    if leaf.majority == nothing
        error("leaf.majority == nothing for $x")
    end
    return leaf.majority
end

function apply_tree(tree::Node, x::DataArray)
    if tree.split_value == nothing                  # true at leaf node, maybe??
        error("WOW!!! Actually made it where tree.split_value == nothing!!!")
        return apply_tree(tree.left, x)
    elseif isna(x[tree.col_idx])
        n_surr = count_surrogates(tree)

        if n_surr > 0               # most often, we do have surrogates

            for i = 1:n_surr
                idx, val = tree.surrogates[i]
                if idx ≠ 0 && !isna(x[idx])

                    if x[idx] < val
                        return apply_tree(tree.left, x)
                    else
                        return apply_tree(tree.right, x)
                    end
                # If there are no surrogates for this predictor we flip a
                # coin to decide left/right.
                elseif idx == 0 || (isna(x[idx]) && i == n_surr)
                    go_left = bitrand(1)[1]
                    if go_left
                        return apply_tree(tree.left, x)
                    else
                        return apply_tree(tree.right, x)
                    end
                end
            end
        # There are corner cases where no initialized surrogates exist
        elseif n_surr == 0
            warn("No surrogates, so we flip a coin for observation: $x")
            go_left = bitrand(1)[1]
            if go_left
                return apply_tree(tree.left, x)
            else
                return apply_tree(tree.right, x)
            end
        end
    # in the simplest case, we have observed value
    elseif x[tree.col_idx] < tree.split_value
        return apply_tree(tree.left, x)
    else
        return apply_tree(tree.right, x)
    end
end

function apply_tree(tree::LeafOrNode, X::DataFrame)
    n = size(X, 1)
    predictions = Array(Any, n)             # this should not be type Any
    for i = 1:n
        # Careful here: conversion to DataArray does not drop dimension,
        # so X is a 1-row matrix; but this this works for now.
        predictions[i] = apply_tree(tree, DataArray(X[i, :])')
        if predictions[i] == nothing
            error("Prediction for row $i == nothing")
        end
    end
    return predictions
end

# This version is dispatched when OOB importance is being calculated
function apply_tree(tree::LeafOrNode, X::DataFrame, y::Vector)
    n = size(X, 1)
    predictions = similar(y)
    for i = 1:n
        # Careful here: conversion to DataArray does not drop dimension,
        # so X is a 1-row matrix; but this this works for now.
        predictions[i] = apply_tree(tree, DataArray(X[i, :])')
        if predictions[i] == nothing
            error("Prediction for row $i == nothing")
        end
    end

    mse_inc_dict = Dict{Int, Float64}()
    preds_used = predictors_used(tree)
    mse_predictions = mean_squared_error(y, predictions)

    for j in preds_used
        X_perm = copy(X)                        # inefficient: find a better way!

        X_perm[:, j] = shuffle(X_perm[:, j])    # some evidence this could be causing randomly occuring seg faults
        permute_predictions = similar(y)

        for i = 1:n
            # Careful here: conversion to DataArray does not drop dimension,
            # so X is a 1-row matrix; but this this works for now.
            permute_predictions[i] = apply_tree(tree, DataArray(X_perm[i, :])')
            if permute_predictions[i] == nothing
                error("Prediction for row $i == nothing")
            end
        end
        mse_inc_dict[j] = mean_squared_error(y, permute_predictions) - mse_predictions
    end
    return (predictions, mse_inc_dict)
end

# Functions for applying the fitted random forest

function apply_forest(forest::Ensemble, X::Matrix)
    n = size(X, 1)
    y_hat = Array(Any, n)
    for i in 1:n
        if VERSION < v"0.5.0-dev"
            y_hat[i] = apply_forest(forest, squeeze(X[i, :], 1))
        else
            y_hat[i] = apply_forest(forest, X[i, :])
        end
    end
    if eltype(y_hat) <: Float64
        return float(y_hat)
    else
        return y_hat
    end
end

function apply_forest(forest::Ensemble, x::Vector)
    ntrees = length(forest)
    votes = Array(Any, ntrees)

    for i in 1:ntrees
        votes[i] = apply_tree(forest.trees[i], x)
    end

    if eltype(votes) <: AbstractFloat
        return mean(votes)
    else
        return majority_vote(votes)
    end
end


text:
include("./src/ShellModel.jl")
using .ShellModel
#include("./src/lanczos_methods.jl")
#include("./src/transit.jl")
#include("./src/input_int_snt.jl")
#include("./src/eigenvector_continuation.jl")

"""
Arguments for main: digonalize the model-space Hamiltonian
sntf:       input interaction file (.snt)
target_nuc: target nucleus
num_ev:     number of eigenstates to be evaluated
target_J    J*2 for the target states 0=>0, 1=>2, 3/2=>3, ...
(Optional)
save_wav=false  # to specify the option to save wavefunctions
is_show = true  # to show elapsed time & allocations
lm = 100        # number of Lanczos vectors
ls = 15         # number of vectors to be used for Thick-Restart
tol= 1.e-6      # tolerance for convergence check in the Lanczos methods
"""
function run(num_ev=10)
    Hs = ["./snts/usdb.snt"]
    intf = "./snts/random_input/usdb.int"

    target_nuc = "Si28"
    tJNs = [ [0,2],[2,2],[4,1],[6,4],[8,3]]

    Erefs = [ [-135.70000,-130.72100], [-127.37200,-126.20400],[-133.92100],
              [-129.42400,-127.90100,-127.11200,-125.49100],
              [-131.08200,-128.81200,-126.53600] ]

    Eexact = [ [-135.861,-131.024], [ -127.669,-126.278],[-133.929],
               [-129.531,-127.852,-126.826,-126.443],
               [-131.254,-128.856,-126.471] ]

    errors = [[0.06764981972995088,0.11807824551232216],[0.14625584538718783,0.28744120626794256],
              [0.05034078995433333],
              [0.10487511604722499,0.23692392369497384, 0.21450202957453257, 0.31933638080241167],
              [0.06909094120663895,0.14593063441913046, 0.26155402077326073]]
    #bnum = 10000 ; itnum = 100000 + bnum
    bnum = 2000; itnum = 10000 + bnum
    #bnum = 10; itnum = 100 + bnum
    #bnum = 50000; itnum = 500000
    solveEC_UQ(Hs,target_nuc,tJNs,Erefs,errors;
               intf=intf,itnum_MCMC=itnum,burnin=bnum,var_M=6.e-2,
               num_replica=1)
    #intf=intf,itnum_MCMC=10000,burnin=1000,var_M=3.e-2)
    #plot_from_history(tJNs,Erefs)
end
run()


text:
"""
    filmpressure!(pressure, height, γ, θ, n, m, hmin, hcrit)

Calculation of the capillary pressure which is given by `` p = - γ∇²h+ Π(h) ``.

# Arguments

- `pressure :: Array{Number,2}`: Array that store the result of the compuation
- `height :: Array{Number,2}`: Height field ``h(\\mathbf{x},t)``
- `γ <: Number`: Forcing strenght due to surface tension
- `θ <: Number`: Equilibrium contact angle
- `n :: Int`: Larger power law exponent for `` Π(h) ``
- `m :: Int`: Smaller power law exponent for `` Π(h) ``
- `hmin <: Number`: Parameter of `` Π(h) ``, in fact `` Π(hmin) = 0 ``
- `hcrit <: Number`: Numerical stabilizer for case `` h(\\mathbf{x},t) \\ll hmin ``

# Mathematics

The capillary pressure ``p_{\\text{cap}}`` is the centeral angle to match our model with the thin film equation.
It consists of two parts, first being the laplace pressure `` \\nabla^2 h `` and second being the derivative of the disjoining pontential `` \\Pi(h) ``

`` p_{\\text{cap}} = -\\gamma \\nabla^2 h + \\Pi(h). ``

For the laplacian term we use the same nine point discretization as in `Swlabe.∇²f!`.
`` \\Pi(h) `` on the other hand is given by

`` \\Pi(h) = \\kappa(\\theta)f(h), ``

where `` \\kappa(\\theta) `` is simply a measure for the **Hamaker constant** and given as

`` \\kappa(\\theta) = \\gamma(1- \\cos(\\theta))\\frac{(n-1)(m-1)}{(n-m)h_{\\text{min}}}.``

For `` f(h) `` one can use various forms, a very common however is the power law given by

`` f(h) = \\bigg[\\bigg(\\frac{h_{\\text{min}}}{h}\\bigg)^n - \\bigg(\\frac{h_{\\text{min}}}{h}\\bigg)^m\\bigg]. ``

# Examples

```jldoctest
julia> using Swalbe, Test

julia> h = reshape(collect(1.0:25.0),5,5) # A dummy height field
5×5 Matrix{Float64}:
 1.0   6.0  11.0  16.0  21.0
 2.0   7.0  12.0  17.0  22.0
 3.0   8.0  13.0  18.0  23.0
 4.0   9.0  14.0  19.0  24.0
 5.0  10.0  15.0  20.0  25.0

julia> pressure = zeros(5,5); θ = 0.0; # Fully wetting substrate

julia> Swalbe.filmpressure!(pressure, h, θ) # default γ = 0.01

julia> result = [30.0 5.0 5.0 5.0 -20;
                 25.0 0.0 0.0 0.0 -25.0;
                 25.0 0.0 0.0 0.0 -25.0;
                 25.0 0.0 0.0 0.0 -25.0;
                 20.0 -5.0 -5.0 -5.0 -30.0];

julia> for i in eachindex(result)
           @test result[i] .≈ -100 .* pressure[i] atol=1e-12
       end
```

# References

- [Peschka et al.](https://www.pnas.org/content/116/19/9275)
- [Craster and Matar](https://journals.aps.org/rmp/abstract/10.1103/RevModPhys.81.1131)
- [Derjaguin and Churaev](https://www.sciencedirect.com/science/article/abs/pii/0021979778900565)

"""
function filmpressure!(output, f, γ, θ, n, m, hmin, hcrit)
    # Straight elements j+1, i+1, i-1, j-1
    hip = circshift(f, (1,0))
    hjp = circshift(f, (0,1))
    him = circshift(f, (-1,0))
    hjm = circshift(f, (0,-1))
    # Diagonal elements
    hipjp = circshift(f, (1,1))
    himjp = circshift(f, (-1,1))
    himjm = circshift(f, (-1,-1))
    hipjm = circshift(f, (1,-1))
    # Disjoining pressure part
    κ = (1 .- cospi.(θ)) .* (n - 1) .* (m - 1) ./ ((n - m) * hmin)

    output .= -γ .* ((2/3 .* (hjp .+ hip .+ him .+ hjm)
                   .+ 1/6 .* (hipjp .+ himjp .+ himjm .+ hipjm)
                   .- 10/3 .* f) .+ κ .* (power_broad.(hmin./(f .+ hcrit), n)
                                       .- power_broad.(hmin./(f .+ hcrit), m)))
    return nothing
end

function filmpressure!(output, f, dgrad, γ, θ, n, m, hmin, hcrit)
    hip, hjp, him, hjm, hipjp, himjp, himjm, hipjm = viewneighbors(dgrad)
    # Straight elements j+1, i+1, i-1, j-1
    circshift!(hip, f, (1,0))
    circshift!(hjp, f, (0,1))
    circshift!(him, f, (-1,0))
    circshift!(hjm, f, (0,-1))
    # Diagonal elements
    circshift!(hipjp, f, (1,1))
    circshift!(himjp, f, (-1,1))
    circshift!(himjm, f, (-1,-1))
    circshift!(hipjm, f, (1,-1))
    #=
    Disjoining pressure part:
    1. Constant part due to angle, n, m, hmin
    2. Part due to the powerlaw
    =#
    output .= -γ .* ((1 .- cospi.(θ)) .* (n - 1) .* (m - 1) ./ ((n - m) * hmin)
                 .* (power_broad.(hmin./(f .+ hcrit), n)
                  .- power_broad.(hmin./(f .+ hcrit), m)) )

    output .-= γ .* (2/3 .* (hjp .+ hip .+ him .+ hjm)
                   .+ 1/6 .* (hipjp .+ himjp .+ himjm .+ hipjm)
                   .- 10/3 .* f)
    return nothing
end
# Standard usage parameters
function filmpressure!(output, f, θ)
    # Straight elements j+1, i+1, i-1, j-1
    hip = circshift(f, (1,0))
    hjp = circshift(f, (0,1))
    him = circshift(f, (-1,0))
    hjm = circshift(f, (0,-1))
    # Diagonal elements
    hipjp = circshift(f, (1,1))
    himjp = circshift(f, (-1,1))
    himjm = circshift(f, (-1,-1))
    hipjm = circshift(f, (1,-1))
    # Disjoining pressure part
    κ = (1 .- cospi.(θ)) .* 16/0.6

    output .= -0.01 .* ((2/3 .* (hjp .+ hip .+ him .+ hjm)
                   .+ 1/6 .* (hipjp .+ himjp .+ himjm .+ hipjm)
                   .- 10/3 .* f) .+ κ .* (power_broad.(0.1./(f .+ 0.05), 9)
                                       .- power_broad.(0.1./(f .+ 0.05), 3)))
    return nothing
end
# One dim implementation
function filmpressure!(output::Vector, f, dgrad, γ, θ, n, m, hmin, hcrit)
    hip, him = viewneighbors_1D(dgrad)
    # Straight elements j+1, i+1, i-1, j-1
    circshift!(hip, f, 1)
    circshift!(him, f, -1)
    #=
    Disjoining pressure part:
    1. Constant part due to angle, n, m, hmin
    2. Part due to the powerlaw
    =#
    output .= -γ .* ((1 .- cospi.(θ)) .* (n - 1) .* (m - 1) ./ ((n - m) * hmin)
                 .* (power_broad.(hmin./(f .+ hcrit), n)
                  .- power_broad.(hmin./(f .+ hcrit), m)) )

    output .-= γ .* (hip .- 2 .* f .+ him)
    return nothing
end

function filmpressure!(output::Vector, f, dgrad, rho, γ, θ, n, m, hmin, hcrit; Gamma=0.0)
    hip, him = viewneighbors_1D(dgrad)
    # Straight elements j+1, i+1, i-1, j-1
    circshift!(hip, f, 1)
    circshift!(him, f, -1)
    #=
    Disjoining pressure part:
    1. Constant part due to angle, n, m, hmin
    2. Part due to the powerlaw
    =#
    output .= - (γ .+ Gamma .* rho) .* ((1 .- cospi.(θ)) .* (n - 1) .* (m - 1) ./ ((n - m) * hmin)
                 .* (power_broad.(hmin./(f .+ hcrit), n)
                  .- power_broad.(hmin./(f .+ hcrit), m)) )

    output .-= (γ .+ Gamma .* rho) .* (hip .- 2 .* f .+ him)
    return nothing
end

"""
    power_broad(arg, n)

Computes `arg` to the power `n`.

Actually this is useful because the `^` operator is much slower.
Same thing I learned about the `pow` function in **C**, * yes it does what you want, but it is slow as fuck *.

# Examples
```jldoctest
julia> using Swalbe, Test

julia> Swalbe.power_broad(3, 3)
27

julia> Swalbe.power_broad.([2.0 5.0 6.0], 2) # Use the broadcasting operator `.`
1×3 Matrix{Float64}:
 4.0  25.0  36.0

```

See also: [`filmpressure`](@ref)
"""
function power_broad(arg::Float64, n::Int)
    temp = 1.0
    for i = 1:n
        temp *= arg
    end
    return temp
end

function power_broad(arg::Float32, n::Int)
    temp = 1.0f0
    for i = 1:n
        temp *= arg
    end
    return temp
end

function power_broad(arg::Int, n::Int)
    temp = 1
    for i = 1:n
        temp *= arg
    end
    return temp
end


text:
module eFEMpart

using Reexport

@reexport using eFEM
@reexport using StokesParticles

export update_particles_noFEMcl!,
       update_particles!,
       generate_femap,
       interpFGT!

export scratch_data

include("struct_defs.jl")
include("barycentric_interpolation.jl")
include("seepage_force.jl")
include("update_particles.jl")
include("fem_sp_coupling.jl")

end # module

max_stars_repo_name: pseastham/eFEMpart

max_issues_repo_name: pseastham/eFEMpart


text:
using DelimitedFiles

# slist = [
    #     "BFFFBBFRRR",
    #     "FFFBBBFRRR",
    #     "BBFFBBFRLL"
    # ]

function to_rowcol(pass, nrow, ncol)
    r1 = 0
    r2 = nrow - 1
    c1 = 0
    c2 = ncol - 1
    # println(seat)
    for d in pass[1:c2]
        dr = r2 - r1

        if d == 'F'
            r2 = floor(r2 - dr/2)
        elseif d == 'B'
            r1 = ceil(r1 + dr/2)
        end
    end

    for d in pass[c2+1:10]
        dc = c2 - c1
        if d == 'L'
            c2 = floor(c2 - dc/2)
        elseif d == 'R'
            c1 = ceil(c1 + dc/2)
        end
    end

    return((r1, c1, (r1*8 + c1)))
end

function process_seat(chart, r, c)
    r = Int(r) + 1
    c = Int(c) + 1

    chart[r,c] = false
    return(chart)
end

function find_seat(chart, nrows)
    for r in 1:nrows
        # println(chart[i,:])
        row = chart[r,:]
        if count(row) == 1
            c = findfirst(x->x==true, row)
            return(r-1 , c-1)
        end
    end
end

nrow = 128
ncol = 8

global chart = trues((nrow, ncol))

slist = readdlm("./seats.dat")

global max_id = 0
for pass in slist
    x,y,id = to_rowcol(pass, nrow, ncol)
    global chart = process_seat(chart, x,y)
    global max_id = max(id, max_id)
end

println("Maximum ", max_id)
myr, myc = find_seat(chart, nrow)
println(myr*8 + myc)


text:
Base.parent(x::ComponentArray) = getfield(x, :data)

Base.size(x::ComponentArray) = size(getdata(x))
ArrayInterface.size(A::ComponentArray) = ArrayInterface.size(parent(A))

Base.elsize(x::Type{<:ComponentArray{T,N,A,Axes}}) where {T,N,A,Axes} = Base.elsize(A)

Base.axes(x::ComponentArray) = axes(getdata(x))

Base.reinterpret(::Type{T}, x::ComponentArray, args...) where T = ComponentArray(reinterpret(T, getdata(x), args...), getaxes(x))

# Cats
# TODO: Make this a little less copy-pastey
function Base.hcat(x::AbstractComponentVecOrMat, y::AbstractComponentVecOrMat)
    ax_x, ax_y = second_axis.((x,y))
    if reduce((accum, key) -> accum || (key in keys(ax_x)), keys(ax_y); init=false) || getaxes(x)[1] != getaxes(y)[1]
        return hcat(getdata(x), getdata(y))
    else
        data_x, data_y = getdata.((x, y))
        ax_y = reindex(ax_y, size(x,2))
        idxmap_x, idxmap_y = indexmap.((ax_x, ax_y))
        axs = getaxes(x)
        return ComponentArray(hcat(data_x, data_y), axs[1], Axis((;idxmap_x..., idxmap_y...)), axs[3:end]...)
    end
end

second_axis(ca::AbstractComponentVecOrMat) = getaxes(ca)[2]
second_axis(::ComponentVector) = FlatAxis()

# Are all these methods necessary?
# TODO: See what we can reduce down to without getting ambiguity errors
Base.vcat(x::ComponentVector, y::AbstractVector) = vcat(getdata(x), y)
Base.vcat(x::AbstractVector, y::ComponentVector) = vcat(x, getdata(y))
function Base.vcat(x::ComponentVector, y::ComponentVector)
    if reduce((accum, key) -> accum || (key in keys(x)), keys(y); init=false)
        return vcat(getdata(x), getdata(y))
    else
        data_x, data_y = getdata.((x, y))
        ax_x, ax_y = getindex.(getaxes.((x, y)), 1)
        ax_y = reindex(ax_y, length(x))
        idxmap_x, idxmap_y = indexmap.((ax_x, ax_y))
        return ComponentArray(vcat(data_x, data_y), Axis((;idxmap_x..., idxmap_y...)))
    end
end
function Base.vcat(x::AbstractComponentVecOrMat, y::AbstractComponentVecOrMat)
    ax_x, ax_y = getindex.(getaxes.((x, y)), 1)
    if reduce((accum, key) -> accum || (key in keys(ax_x)), keys(ax_y); init=false) || getaxes(x)[2:end] != getaxes(y)[2:end]
        return vcat(getdata(x), getdata(y))
    else
        data_x, data_y = getdata.((x, y))
        ax_y = reindex(ax_y, size(x,1))
        idxmap_x, idxmap_y = indexmap.((ax_x, ax_y))
        return ComponentArray(vcat(data_x, data_y), Axis((;idxmap_x..., idxmap_y...)), getaxes(x)[2:end]...)
    end
end
Base.vcat(x::CV...) where {CV<:AdjOrTransComponentArray} = ComponentArray(reduce(vcat, map(y->getdata(y.parent)', x)), getaxes(x[1]))
Base.vcat(x::ComponentVector, args...) = vcat(getdata(x), getdata.(args)...)
Base.vcat(x::ComponentVector, args::Vararg{AbstractVector{T}, N}) where {T,N} = vcat(getdata(x), getdata.(args)...)

function Base.hvcat(row_lengths::Tuple{Vararg{Int}}, xs::AbstractComponentVecOrMat...)
    i = 1
    idxs = UnitRange{Int}[]
    for row_length in row_lengths
        i_last = i + row_length - 1
        push!(idxs, i:i_last)
        i = i_last + 1
    end
    rows = [reduce(hcat, xs[idx]) for idx in idxs]
    return vcat(rows...)
end

function Base.permutedims(x::ComponentArray, dims)
    axs = getaxes(x)
    return ComponentArray(permutedims(getdata(x), dims), map(i->axs[i], dims)...)
end

## Indexing
Base.IndexStyle(::Type{<:ComponentArray{T,N,<:A,<:Axes}}) where {T,N,A,Axes} = IndexStyle(A)

# Since we aren't really using the standard approach to indexing, this will forward things to
# the correct methods
Base.to_indices(x::ComponentArray, i::Tuple) = i
Base.to_indices(x::ComponentArray, i::Tuple{Vararg{Union{Integer, CartesianIndex}, N}}) where N = i
Base.to_indices(x::ComponentArray, i::Tuple{Vararg{Int64}}) where N = i
Base.to_index(x::ComponentArray, i) = i

# Get ComponentArray index
Base.@propagate_inbounds Base.getindex(x::ComponentArray, idx::CartesianIndex) = getdata(x)[idx]
Base.@propagate_inbounds Base.getindex(x::ComponentArray, idx::FlatIdx...) = getdata(x)[idx...]
Base.@propagate_inbounds function Base.getindex(x::ComponentArray, idx::FlatOrColonIdx...)
    axs = map((ax, i) -> getindex(ax, i).ax, getaxes(x), idx)
    axs = remove_nulls(axs...)
    return ComponentArray(getdata(x)[idx...], axs...)
end
Base.@propagate_inbounds Base.getindex(x::ComponentArray, ::Colon) = getdata(x)[:]
Base.@propagate_inbounds Base.getindex(x::ComponentArray, ::Colon...) = x
@inline Base.getindex(x::ComponentArray, idx...) = getindex(x, toval.(idx)...)
@inline Base.getindex(x::ComponentArray, idx::Val...) = _getindex(getindex, x, idx...)

# Set ComponentArray index
Base.@propagate_inbounds Base.setindex!(x::ComponentArray, v, idx::FlatOrColonIdx...) = setindex!(getdata(x), v, idx...)
Base.@propagate_inbounds Base.setindex!(x::ComponentArray, v, ::Colon) = setindex!(getdata(x), v, :)
@inline Base.setindex!(x::ComponentArray, v, idx...) = setindex!(x, v, toval.(idx)...)
@inline Base.setindex!(x::ComponentArray, v, idx::Val...) = _setindex!(x, v, idx...)

# Explicitly view
Base.@propagate_inbounds Base.view(x::ComponentArray, idx::ComponentArrays.FlatIdx...) = view(getdata(x), idx...)
Base.@propagate_inbounds Base.view(x::ComponentArray, idx...) = _getindex(view, x, toval.(idx)...)

Base.@propagate_inbounds Base.maybeview(x::ComponentArray, idx::ComponentArrays.FlatIdx...) = Base.maybeview(getdata(x), idx...)
Base.@propagate_inbounds Base.maybeview(x::ComponentArray, idx...) = _getindex(Base.maybeview, x, toval.(idx)...)

# Generated get and set index methods to do all of the heavy lifting in the type domain
@generated function _getindex(index_fun, x::ComponentArray, idx...)
    ci = getindex.(getaxes(x), getval.(idx))
    inds = map(i -> i.idx, ci)
    axs = map(i -> i.ax, ci)
    axs = remove_nulls(axs...)
    return :(Base.@_inline_meta; ComponentArray(index_fun(getdata(x), $inds...), $axs...))
end

@generated function _setindex!(x::ComponentArray, v, idx...)
    ci = getindex.(getaxes(x), getval.(idx))
    inds = map(i -> i.idx, ci)
    return :(Base.@_inline_meta; setindex!(getdata(x), v, $inds...))
end

## Linear Algebra
Base.pointer(x::ComponentArray{T,N,A,Axes}) where {T,N,A<:DenseArray,Axes} = pointer(getdata(x))

Base.unsafe_convert(::Type{Ptr{T}}, x::ComponentArray{T,N,A,Axes}) where {T,N,A,Axes} = Base.unsafe_convert(Ptr{T}, getdata(x))

Base.strides(x::ComponentArray) = strides(getdata(x))
ArrayInterface.strides(A::ComponentArray) = ArrayInterface.strides(parent(A))
for f in [:device, :stride_rank, :contiguous_axis, :contiguous_batch_size, :dense_dims]
    @eval ArrayInterface.$f(::Type{ComponentArray{T,N,A,Axes}}) where {T,N,A,Axes} = ArrayInterface.$f(A)
end

Base.stride(x::ComponentArray, k) = stride(getdata(x), k)
Base.stride(x::ComponentArray, k::Int64) = stride(getdata(x), k)

ArrayInterface.lu_instance(jac_prototype::ComponentArray) = ArrayInterface.lu_instance(getdata(jac_prototype))

ArrayInterface.parent_type(::Type{ComponentArray{T,N,A,Axes}}) where {T,N,A,Axes} = A

for f in [*, \, /]
    op = nameof(f)
    @eval begin
        function Base.$op(A::ComponentVecOrMat, B::ComponentMatrix)
            C = $op(getdata(A), getdata(B))
            ax1 = getaxes(A)[1]
            ax2 = getaxes(B)[2]
            return ComponentArray(C, (ax1, ax2))
        end
        function Base.$op(A::ComponentVecOrMat, b::ComponentVector)
            c = $op(getdata(A), getdata(b))
            ax1 = getaxes(A)[1]
            return ComponentArray(c, ax1)
        end
        function Base.$op(Aᵀ::ComponentMatrix{T}, B::AdjOrTransComponentVecOrMat{T}) where {T}
            Cᵀ = $op(getdata(Aᵀ), getdata(B))
            ax1 = getaxes(Aᵀ)[1]
            ax2 = getaxes(B)[2]
            return ComponentArray(Cᵀ, ax1, ax2)
        end
    end
    for (adjfun, AdjType) in zip([adjoint, transpose], [Adjoint, Transpose])
        adj = nameof(adjfun)
        Adj = nameof(AdjType)
        @eval begin
            function Base.$op(aᵀ::$Adj{T,<:ComponentVector}, B::ComponentMatrix) where {T}
                cᵀ = parent($op(getdata(aᵀ), getdata(B)))
                ax2 = getaxes(B)[2]
                return $adj(ComponentArray(cᵀ, ax2))
            end
            function Base.$op(Aᵀ::$Adj{T,<:ComponentMatrix}, B::AdjOrTransComponentVecOrMat) where {T}
                Cᵀ = $op(getdata(Aᵀ), getdata(B))
                ax1 = getaxes(Aᵀ)[1]
                ax2 = getaxes(B)[2]
                return ComponentArray(Cᵀ, ax1, ax2)
            end
            function Base.$op(Aᵀ::$Adj{T,<:ComponentMatrix}, B::ComponentMatrix{T}) where {T}
                Cᵀ = $op(getdata(Aᵀ), getdata(B))
                ax1 = getaxes(Aᵀ)[1]
                ax2 = getaxes(B)[2]
                return ComponentArray(Cᵀ, ax1, ax2)
            end
            function Base.$op(Aᵀ::$Adj{T,<:ComponentMatrix}, b::ComponentVector) where {T}
                cᵀ = $op(getdata(Aᵀ), getdata(b))
                ax1 = getaxes(Aᵀ)[1]
                return ComponentArray(cᵀ, ax1)
            end
        end
    end
end

# While there are some cases where these were faster, it is going to be almost impossible to
# to keep up with method ambiguity errors due to other array types overloading *, /, and \.
# Leaving these here and commented out for now, but will delete them later.

# # Avoid slower fallback
# for f in [:(*), :(/), :(\)]
#     @eval begin
#         # The normal stuff
#         Base.$f(x::ComponentArray, y::AbstractArray) = $f(getdata(x), y)
#         Base.$f(x::AbstractArray, y::ComponentArray) = $f(x, getdata(y))
#         Base.$f(x::ComponentArray, y::ComponentArray) = $f(getdata(x), getdata(y))

#         # A bunch of special cases to avoid ambiguous method errors
#         Base.$f(x::ComponentArray, y::AbstractMatrix) = $f(getdata(x), y)
#         Base.$f(x::AbstractMatrix, y::ComponentArray) = $f(x, getdata(y))

#         Base.$f(x::ComponentArray, y::AbstractVector) = $f(getdata(x), y)
#         Base.$f(x::AbstractVector, y::ComponentArray) = $f(x, getdata(y))
#     end
# end

# # Adjoint/transpose special cases
# for f in [:(*), :(/)]
#     @eval begin
#         Base.$f(x::Adjoint, y::ComponentArray) = $f(getdata(x), getdata(y))
#         Base.$f(x::Transpose, y::ComponentArray) = $f(getdata(x), getdata(y))

#         Base.$f(x::Adjoint{T,<:AbstractVector{T}}, y::ComponentVector) where T = $f(x, getdata(y))
#         Base.$f(x::Transpose{T,<:AbstractVector{T}}, y::ComponentVector) where T = $f(x, getdata(y))

#         Base.$f(x::Adjoint{T,<:AbstractVector{T}}, y::ComponentMatrix{T,A,Axes}) where {T,A,Axes} = $f(x, getdata(y))
#         Base.$f(x::Transpose{T,<:AbstractVector{T}}, y::ComponentMatrix{T,A,Axes}) where {T,A,Axes} = $f(x, getdata(y))

#         Base.$f(x::Adjoint{T,<:AbstractMatrix{T}}, y::ComponentVector) where {T} = $f(x, getdata(y))
#         Base.$f(x::Transpose{T,<:AbstractMatrix{T}}, y::ComponentVector) where {T} = $f(x, getdata(y))

#         Base.$f(x::ComponentArray, y::Adjoint{T,<:AbstractVector{T}}) where T = $f(getdata(x), y)
#         Base.$f(x::ComponentArray, y::Transpose{T,<:AbstractVector{T}}) where T = $f(getdata(x), y)

#         Base.$f(x::ComponentArray, y::Adjoint{T,<:ComponentVector}) where T = $f(getdata(x), getdata(y))
#         Base.$f(x::ComponentArray, y::Transpose{T,<:ComponentVector}) where T = $f(getdata(x), getdata(y))

#         # There seems to be a new method in Julia > v.1.4 that specializes on this
#         Base.$f(x::Adjoint{T,<:AbstractVector{T}}, y::ComponentVector{T,A,Axes}) where {T<:Number,A,Axes} = $f(x, getdata(y))
#         Base.$f(x::Transpose{T,<:AbstractVector{T}}, y::ComponentVector{T,A,Axes}) where {T<:Number,A,Axes} = $f(x, getdata(y))

#         Base.$f(x::Adjoint{T,<:AbstractVector{T}}, y::ComponentVector{T,A,Axes}) where {T<:Real,A,Axes} = $f(getdata(x), getdata(y))
#         Base.$f(x::Transpose{T,<:AbstractVector{T}}, y::ComponentVector{T,A,Axes}) where {T<:Real,A,Axes} = $f(getdata(x), getdata(y))
#     end
# end


text:
#
# This file is a part of MolecularGraph.jl
# Licensed under the MIT License http://opensource.org/licenses/MIT
#

export
    Canvas, Color,
    singlebond!,
    wedged!,
    dashedwedged!,
    wavesingle!,
    doublebond!,
    clockwisedouble!,
    counterdouble!,
    crossdouble!,
    triplebond!,
    atomsymbolright!,
    atomsymbolcenter!,
    atomsymbolleft!,
    atom_annotation!

struct Color
    r::Int
    g::Int
    b::Int
end

abstract type Canvas end

"""
    setbond!(
        canvas::Canvas, bondorder::Int, bondnotation::Int, coords::Segment2D,
        ucolor::Color, vcolor::Color, u_visible::Bool, v_visible::Bool
    )

Interface for bond drawing.
"""
function setbond! end

"""
    setatomright!(
        canvas::Canvas, coords::Array{Float64,2}, atomsymbol::Symbol,
        color::Color, implicithconnected::Int, charge::Int
    )

Interface for atom drawing.
"""
function setatomright! end
function setatomcenter! end
function setatomleft! end

"""
    setatomnote!(
        canvas::Canvas, coords::Array{Float64,2}, text::String, color::Color,
        bgcolor::Color
    )

Interface for atom note drawing.
"""
function setatomnote! end


text:
export eigbounds, eigval, eigenvector

# determine lower and upper bounds for eigenvalues number k1:k2.
# Bisectional method for function bif, where bif(x) is the number of eigenvalues
# less than or equal to x.
# scale is an (under-) estimation for the maximal abolute value of eigenvectors. It is used
# to control interval size growth for smallest and largest ev.
function eigbounds(k1::Int, k2::Int, ws, scale::T; rtol=T(Inf), atol=T(Inf), rtolg=eps(T), atolg=eps(scale)) where T<:AbstractFloat
    n = k2 - k1 + 1
    lb = Vector{T}(undef, n)
    ub = Vector{T}(undef, n)
    eigbounds!(lb, ub, k1, k2, ws, scale, rtol, atol, rtolg, atolg)
end
function eigbounds!(lb::V, ub::V, k1::Int, k2::Int, ws, scale::T, rtol, atol, rtolg, atolg) where {T<:AbstractFloat,V<:AbstractVector{T}}
    n = k2 - k1 + 1
    length(lb) == n == length(ub) || throw(ArgumentError("length mismatch"))
    fill!(lb, -T(Inf))
    fill!(ub, T(Inf))
    scalep = abs(scale)
    scalen = -abs(scale)
    if isinf(atol) != isinf(rtol)
        atol = ifelse(isinf(atol), zero(atol), atol)
        rtol = ifelse(isinf(rtol), zero(rtol), rtol)
    end
    # find (first) index j with ub[j] - lb[j+1] > tol
    function findgap(rtol, atol)
        for j = 1:n-1
            dx = ub[j] - lb[j+1]
            if dx >= max(max(abs(ub[j]), abs(lb[j]), ub[j+1] - lb[j]) * rtol, atol)
                return j
            end
        end
        return 0
    end

    # find (first) index j with ub[j] - lb[j+1] > tol
    function findbad(rtol, atol)
        for j = 1:n
            dx = ub[j] - lb[j]
            if abs(dx) > max(max(abs(ub[j]), abs(lb[j])) * rtol, atol)
                return j
            end
        end
        return 0
    end

    # midpoint for finite a, b, otherwise extend towards infinity
    function midpoint(a, b)
        a + (b - a) / 2
    end

    setx!(lb, ub, 0.0, k1, ws)
    while isinf(lb[1])
        x = ub[1] + scalen
        scalen *= 2
        setx!(lb, ub, x, k1, ws)
    end
    while isinf(ub[n])
        x = lb[n] + scalep
        scalep *= 2
        setx!(lb, ub, x, k1, ws)
    end
    while ( j = findgap(rtolg, atolg) ) != 0
        x = midpoint(lb[j+1], ub[j])
        setx!(lb, ub, x, k1, ws)
    end
    while ( j = findbad(rtol, atol) ) != 0
        x = midpoint(lb[j], ub[j])
        setx!(lb, ub, x, k1, ws)
    end
    lb, ub
end

# evaluate bif(x) and improve bounds according to result.
function setx!(lb, ub, x, k1, ws)
    n = length(lb)
    k, = detderivates!(ws, x)
    k = k - k1 + 1
    for j =  max(k-k1+1, 0):n-1
        if x > lb[j+1]
            lb[j+1] = x
        else
            break
        end
    end
    for j = min(n, k):-1:1
        if x < ub[j]
            ub[j] = x
        else
            break
        end
    end
end

function laguerre1(η::T, ζ::T, n::Int, r) where T<:AbstractFloat
    disc = max(ζ * n - η^2, 0)
    sq = sqrt(disc * (n - r) / r)
    n / ( η + copysign(sq, η) )
end

function laguerre2(η::T, ζ::T, n::Int, r) where T<:AbstractFloat
    η / ζ
end

function eigval(A, B, k::Int, a::T, b::T, r=1) where T<:AbstractFloat
    ws = make_workspace(A, B)
    eigval!(ws, k, a, b, r)
end

function eigval!(ws::Workspace{T}, k::Int, a::S, b::S, r, check=true) where {S,T}
    n = size(ws.A, 1)
    a, b = promote(a, b, zero(real(T)))
    κa = k - 1; κb = k
    if check
        κa, = detderivates!(ws, a)
        κb, = detderivates!(ws, b)
        κa < k <= κb || throw(ArgumentError("eigenvalue($k) not in [$a,$b]"))
        r = max(κb - κa, 1)
        println("start: κa=$κa κb=$κb r=$r")
    end
    println("start: κa=$κa κb=$κb r=$r")
    x = a + (b - a) / 2
    x1 = real(T)(NaN)
    step = 0
    η0 = NaN
    while step < 100
        step += 1
        κ, η, ζ0 = detderivates!(ws, x)
        if κ < k
            a = x
            κa = κ
            r = max(κb - κa, 1)
        else
            b = x
            κb = κ
            r = max(κb - κa, 1)
        end
        println("$step: κ=$κ η=$η ζ=$ζ0 r=$r")
        x0 = x1
        x1 = x
        if κ < k && η <= 0 || κ >= k && η >= 0 # Newton points to ev[k]
            if isnan(η0)
                ζ = zero(η)
            else
                Δ = (inv(η0) - inv(η)) / (x0 - x)
                ζ = (1 - Δ) * η^2
            end
            if k - 1 <= κ <= k #
                η0 = η
                dx = laguerre1(η, ζ0, n, r)
                if dx * η >= 0.5
                    x = min(b, max(a, x1 - dx))
                    op = "Laguerre r=$r"
                elseif isnan(dx) && (dx = inv(η)) |> isfinite && abs(dx) <  (b - a) / 2
                    x = min(b, max(a, x1 - dx))
                    op = "Newton"
                else
                    x = a + (b - a) / 2
                    op = "bisect1"
                end
            else
                x = κ > k ? a + (b - a) / (κ - k + 1) : b - (b - a) / (k - κ)
                op = "division $κ / $k"
            end
        else
            x = a + (b - a) / 2
            op = "bisect2"
        end
        dx = x - x1
        println("step $step: $x $dx $op")
        convergence2(x, x1, x0) && break
    end
    x
end

"""
    convergence(x2, x1, x0)

Determines convergence of a series, given 3 succcessive elements.

returns `true`, if `x1 == x2`
returns `true`, if not monotonous and `|x2 - x1| >= |x1 - x0|`
returns `true`, if `|x2 - x1|^2 <= (|x1 - x0| - |x2 - x1|) * ϵ
returns `false` in all other cases
"""
function convergence2(x2, x1, x0)
    ϵ = eps(x2)
    dx2 = abs(x2 - x1)
    x1 == x2 && return true
    dx1 = abs(x1 - x0)
    dx2 >= dx1 && (x2 > x1) != (x1 > x0) && return true
    dx2^2 <= (dx1 - dx2) * ϵ && return true
    return false
end

function eigenvector(A, B, k, a, b, r)
    ws = make_workspace(A, B)
    x = eigval!(ws, k, a, b, r)
    L = zeros(size(ws.Q, 1), size(A, 1), 2)
    detderivates!(ws, x, L)
    _eigenvector(x, L)
end

function _eigenvector(x, L::Array{T}) where T
    q, n = size(L)
    v = Vector{T}(undef, n)
    w = Vector{T}(undef, n)
    v[n] = one(T)
    w[n] = zero(T)
    for i = n-1:-1:1
        s = zero(T)
        t = zero(T)
        for j = 2:min(q, n - i + 1)
            s -= v[i+j-1] * L[j,i,1]
            t -= v[i+j-1] * L[j,i,2] + w[i+j-1] * L[j,i,1]
        end
        v[i] = s
        w[i] = t
    end
    u = inv(sqrt(norm(v)))
    v .*= u
    w .*= u
    x, v, w
end


text:
# All functions related to reproduction.

"""
    reproduce!(patch, settings)

Reproduction of individuals in a patch (default function).
"""
function reproduce!(patch::Patch, settings::Dict{String, Any}) #TODO: refactor!
    #TODO This is one of the most compute-intensive functions in the model - optimise?
    for ind in patch.community
        ind.marked && continue # individual might not have established yet
        ind.size < ind.traits["repsize"] && continue
        metaboffs = settings["fertility"] * ind.size^(-1/4) * exp(-act/(boltz*patch.temp))
        noffs = rand(Poisson(metaboffs))
        noffs < 1 && continue
        partners = findmate(patch.community, ind, settings)
        if length(partners) < 1 && rand() < ind.traits["selfing"]
            partners = [ind]
        elseif length(partners) < 1
            continue
        end
        numpartners = Integer(round(ind.traits["numpollen"]))
        if numpartners == 0
            #FIXME happens regularly!
            # `numpollen` not handled specifically in `createtraits()`?
            simlog("Individual cannot reproduce, `numpollen` too low.", 'd')
            continue
        end
        for ptn in 1:numpartners
            partner = rand(partners)
            # subtract offspring mass from parent #XXX is that sensible, does it make a difference?
            parentmass = ind.size - noffs * ind.traits["seedsize"]
            if parentmass <= 0
                break
            else
                ind.size = parentmass
            end
            append!(patch.seedbank, createoffspring(noffs, ind, partner, settings["traitnames"]))
        end
    end
    simlog("Patch $(patch.id): $(length(patch.seedbank)) offspring", 'd')
end

"""
    greproduce!(patch, settings)

Reproduction of individuals in a patch with global mating.
"""
function greproduce!(patch::Patch, settings::Dict{String, Any})
    for ind in patch.community
        ind.marked && continue # individual might not have established yet
        ind.size < ind.traits["repsize"] && continue
        metaboffs = settings["fertility"] * ind.size^(-1/4) * exp(-act/(boltz*patch.temp))
        noffs = rand(Poisson(metaboffs))
        noffs < 1 && continue
        partners = findmate([(map(x -> x.community, world)...)...], ind, settings)
        if length(partners) < 1 && rand() < ind.traits["selfing"]
            partners = [ind]
        elseif length(partners) < 1
            continue
        end
        numpartners = Integer(round(ind.traits["numpollen"]))
        for ptn in 1:numpartners
            partner = rand(partners, 1)[1]
            parentmass = ind.size - noffs * ind.traits["seedsize"] # subtract offspring mass from parent
            if parentmass <= 0
                continue
            else
                ind.size = parentmass
            end
            append!(patch.seedbank, createoffspring(noffs, ind, partner, settings["traitnames"]))
        end
    end
    simlog("Patch $(patch.id): $(length(patch.seedbank)) offspring", 'd')
end

"""
    reproduce!(world, settings)

Carry out reproduction on all patches.
"""
function reproduce!(world::Array{Patch,1}, settings::Dict{String, Any})
    for patch in world
        if settings["globalmating"]
            greproduce!(patch, settings)
        elseif settings["mode"] == "zosterops"
            zreproduce!(patch, settings)
        else
            (patch.isisland || !settings["static"]) && reproduce!(patch, settings) # pmap(!,patch) ???
        end
    end
end

"""
    findmate(population, individual, settings)

Find a reproduction partner for the given individual in the given population.
"""
function findmate(population::AbstractArray{Individual, 1}, ind::Individual, settings::Dict{String, Any})
    #XXX This function is pretty expensive
    indstate = ind.marked
    ind.marked = true
    mates = Individual[]
    startidx = rand(eachindex(population))
    mateidx = startidx
    while true
        mate = population[mateidx]
        if !mate.marked && iscompatible(mate, ind, settings)
            push!(mates, mate)
            break
        end
        mateidx += 1
        mateidx > length(eachindex(population)) && (mateidx = 1)
        mateidx == startidx && break
    end
    ind.marked = indstate
    mates
end

"""
    createoffspring(noffs, individual, partner, traitnames, dimorphism)

The main reproduction function. Take two organisms and create the given number
of offspring individuals. Returns an array of individuals.
"""
function createoffspring(noffs::Integer, ind::Individual, partner::Individual, traitnames::Array{String, 1}, dimorphism::Bool=false)
    offspring = Individual[]
    for i in 1:noffs # pmap? this loop could be factorized!
        partnergenome = meiosis(partner.genome, false) # offspring have different genome!
        mothergenome = meiosis(ind.genome, true)
        (isempty(partnergenome) || isempty(mothergenome)) && continue
        genome = vcat(partnergenome,mothergenome)
        traits = gettraitdict(genome, traitnames)
        marked = true
        fitness = 0.0
        newpartner = 0
        newsize = ind.traits["seedsize"]
        newid = rand(UInt32) #XXX Is this collision-safe?
        sex = hermaphrodite
        if dimorphism
            rand(Bool) ? sex = male : sex = female
        end
        push!(offspring, Individual(ind.lineage, genome, traits, marked, fitness,
                                    fitness, newsize, sex, newpartner, newid))
    end
    offspring
end

max_stars_repo_name: lleiding/gemm

max_issues_repo_name: lleiding/gemm

max_forks_repo_name: lleiding/gemm