text:
module BioXP

include("structs.jl")
include("randomize.jl")
include("parser.jl")
include("format.jl")
include("dg.jl")

import JSON
using SparseArrays
# using .Structs, .Randomize, .Parser

export
# Functions
readexpansion,
readmaster,
readcompounds,
readids,
readkeyedids,
randomizecompounds,
expand,
find_minimal_seed_set,
simple_write_out,
formatbioxpoutput,
formatexpansion,
filter_reactions_by_dg,

# Structs
System,
Reaction,
Compound,
Expansion,

# Consts
Reactions,
Compounds,
IDs

## reaction_structs only give detailed info about reactions.
##   they are never supposed to be used as a comprehensive set of reactions
##   for the expansion process

function remove_rids_not_in_rstructs(
    rstructs::Reactions,
    rids::IDs)

    collect(String,intersect(rids,keys(rstructs)))
end

function matrixify_compounds(
    rstructs::Reactions,
    cids::IDs,
    rids::IDs)

    ## Reactons in columns, Compounds in rows. Rows, columns.
    R = spzeros(Int, length(cids),length(rids)) ## Reactants
    P = spzeros(Int, length(cids),length(rids)) ## Products

    cdict = Dict([(k,v) for (v,k) in enumerate(cids)])

    for (i,r) in enumerate(rids)
        for cpd in rstructs[r].left
            R[cdict[cpd],i] = 1
        end
        for cpd in rstructs[r].right
            P[cdict[cpd],i] = 1
        end
    end

    R, P
end

function matrixify_targets(
    cids::IDs,
    tids::IDs)

    biosystem_tids = intersect(tids,cids) ## ensures impossible to reach tids don't count

    sparsevec([Int(i in biosystem_tids) for i in cids])

end
"""
Return vector of compounds with 1 in position where seeds are present and 0 where seeds are absenst.
"""
function matrixify_seeds(
    sids::IDs,
    cids::IDs)

    sparsevec([Int(i in sids) for i in cids])
end

function identify_biosystem_compounds(
    rstructs::Reactions,
    rids::IDs)

    ## only look at rids which are in my rstructs

    IDs(unique(Iterators.flatten([vcat(rstructs[r].left,rstructs[r].right) for r in rids])))
end

"""
If we want to take into account dg values, dispatch on this function
"""
function expandmatrices(
    R::SparseMatrixCSC{Int64,Int64},
    P::SparseMatrixCSC{Int64,Int64},
    x::SparseVector{Int64,Int64};
    allowed_forward::Union{SparseVector{Bool,Int64},Nothing}=nothing,
    allowed_backward::Union{SparseVector{Bool,Int64},Nothing}=nothing)

    # initialize variables:
    RT = sparse(transpose(R))
    PT = sparse(transpose(P))

    br = sparsevec(sum(RT, dims=2)) # sum the rows of RT. Do I need the vec call here? yes, turns it into 1 element array
    bp = sparsevec(sum(PT, dims=2))

    # force unallowed reactions to be out of reach
    if !isnothing(allowed_forward)
        br = br + .~allowed_forward
    end
    if !isnothing(allowed_backward)
        bp = bp + .~allowed_backward
    end

    # find the total number of metabolites in the seed set
    k = sum(x);

    # initialize previous iteration count of metabolites in the network
    k0 = 0;

    # iteration 1 consistes of the seed set
    X = SparseVector{Int64,Int64}[x]

    # initialize reaction accumulation matrix
    Y = SparseVector{Int64,Int64}[]

    # while the metabolite set has not converged
    while k > k0
        # update previous number of metabolites
        k0 = k;

        # RT*x ==> represnts the number of present metabolites in the
        # network within each reaction; if this isequal to the total
        # number of metabolites in each reaction, then the reaction
        # is added to the network

        ## Forward reactions
        y = RT * x .== br
        ## Backward reactions
        yp = PT * x .== bp

        # P*y > 0 ==> represents the vector of reactions that produce
        # metabolite i. (i in 1:m).  If this is >0,
        # then that metabolite is producable

        ## Forward reactions
        xnew = P * y .> 0

        ## Backward reactions
        xnewp = R * yp .> 0

        # add to previous set of metabolites (only needed to retain seed set)

        ## Add forward and backward reactions
        x = (x .| xnew .| xnewp)
        y = SparseMatrixCSC{Int64,Int64}(y .| yp)

        # find new total number of metabolites in network
        k = sum(x);

        # append accumulation matricies
        push!(X, x)
        push!(Y, y)
    end
    X, Y
end

function expand(
    system::System;
    write_path::Union{String,Nothing}=nothing)

    expand(system.rstructs,
        system.rids,
        system.sids,
        system.tids,
        write_path = write_path,
        allowed_forward = system.allowed_forward,
        allowed_backward = system.allowed_backward)
end

function expand(
    rstructs::Reactions,
    rids::IDs,
    sids::IDs;
    tids::IDs=IDs(),
    write_path::Union{String,Nothing}=nothing,
    allowed_forward::Union{SparseVector{Bool,Int64},Nothing}=nothing,
    allowed_backward::Union{SparseVector{Bool,Int64},Nothing}=nothing)

    rids = remove_rids_not_in_rstructs(rstructs,rids)
    cids = identify_biosystem_compounds(rstructs,rids)

    (R, P) = matrixify_compounds(rstructs,cids,rids)
    x = matrixify_seeds(sids, cids)
    t = matrixify_targets(cids,tids)

    # X, Y = Vector{Int}[], Vector{Int}[] ## same as Vector{Vector{Int}}(),Vector{Vector{Int}}()
    X, Y = expandmatrices(R, P, x, allowed_forward = allowed_forward, allowed_backward = allowed_backward)

    if write_path !== nothing
        simple_write_out(write_path,x,t,cids,rids,X,Y)
    end

    x, t, cids, X, Y
end

"""
Return indices of seeds within the compounds vector. TEST
"""
function seed_indicies(sids::IDs, cids::IDs)
    # This is a generator, not an array. You can iterate over this thing exactly once
    # because it only stores the current state and what it needs to find the next state.
    [findfirst(isequal(c), cids) for c in sids]
end

function find_minimal_seed_set(
    system::System;
    write_path::Union{String,Nothing}=nothing,
    allowed_forward::Union{SparseVector{Bool,Int64},Nothing}=nothing,
    allowed_backward::Union{SparseVector{Bool,Int64},Nothing}=nothing)

    find_minimal_seed_set(system.rstructs,
        system.rids,
        system.sids,
        tids = system.tids,
        write_path = write_path,
        allowed_forward = allowed_forward,
        allowed_backward = allowed_backward)
end

"""
    find_minimal_seed_set(rstructs,rids,sids,tids,write_path)

Return:
- seeds (binary version) (these must be provided because their ordering matters for the algorithm)
- targets (binary version)
- compound ids
- reactions ids
- compounds (binary version) present at each timestep
- reactions (binary version) availabel at each timestep
"""
function find_minimal_seed_set(
    rstructs::Reactions,
    rids::IDs,
    sids::IDs;
    tids::IDs=IDs(),
    write_path::Union{String,Nothing}=nothing,
    allowed_forward::Union{SparseVector{Bool,Int64},Nothing}=nothing,
    allowed_backward::Union{SparseVector{Bool,Int64},Nothing}=nothing)

    rids = remove_rids_not_in_rstructs(rstructs,rids)
    cids = identify_biosystem_compounds(rstructs,rids)
    # println(cids)
    # println(sids)

    (R, P) = matrixify_compounds(rstructs,cids,rids)
    t = matrixify_targets(cids,tids)
    # X, Y = Vector{Int}[], Vector{Int}[]
    x = matrixify_seeds(sids, cids) ## This should be a vector of all 1s of length(cids)
    # println(x)
    # println(typeof(x))
    # println(ones(Int,length(cids)))
    # println(typeof(ones(Int,length(cids))))
    Vector(x) != ones(Int,length(cids)) && throw(DomainError("This should be a vector of all 1s of length(cids"))
    X, Y, x = loop_and_remove_seeds(sids,cids,x,t,R,P,allowed_forward = allowed_forward, allowed_backward = allowed_backward)

    if write_path !== nothing
        simple_write_out(write_path,x,t,cids,rids,X,Y)
    end

    x,t,cids,rids,X,Y

end

"""
    find_minimal_seed_set()

Return system variables after finding many minimal seed sets.
"""
function find_minimal_seed_set(
    rstructs::Reactions,
    rids::IDs,
    sid_sets::Vector{IDs};
    tids::IDs=IDs(),
    write_path::Union{String,Nothing}=nothing,
    allowed_forward::Union{SparseVector{Bool,Int64},Nothing}=nothing,
    allowed_backward::Union{SparseVector{Bool,Int64},Nothing}=nothing)

    rids = remove_rids_not_in_rstructs(rstructs,rids)
    cids = identify_biosystem_compounds(rstructs,rids)
    (R, P) = matrixify_compounds(rstructs,cids,rids)
    t = matrixify_targets(cids,tids)
    # X, Y = Vector{Int}[], Vector{Int}[]
    all_seed_results = Vector{}
    if isnothing(write_path) ## no parallel processing because of the push

        for (i,sids) in enumerate(sid_sets)
            x = matrixify_seeds(sids, cids) ## This should be a vector of all 1s of length(cids)
            Vector(x) !== ones(Int,length(cids)) && throw(DomainError("This should be a vector of all 1s of length(cids"))
            X, Y, x = loop_and_remove_seeds(sids,cids,x,t,R,P,allowed_forward = allowed_forward, allowed_backward = allowed_backward)
            push!(all_seed_results,(x,t,cids,rids,X,Y))
        end

    else ## Use parallel processing

        Threads.@threads for (i,sids) in collect(enumerate(sid_sets))
            x = matrixify_seeds(sids, cids) ## This should be a vector of all 1s of length(cids)
            # println(length(x))
            # println(length(cids))
            # println(x)
            ## I'm commenting this out for now because it's throwing even though it shouldn't be based on the above prints
            # x !== ones(Int,length(cids)) && throw(DomainError("This should be a vector of all 1s of length(cids"))
            X, Y, x = loop_and_remove_seeds(sids,cids,x,t,R,P,allowed_forward = allowed_forward, allowed_backward = allowed_backward)
            simple_write_out(joinpath(write_path,"$i.json"),x,t,cids,rids,X,Y)
        end

    end

    all_seed_results # This will be empty if write_path!=nothing

end

function find_minimal_seed_set(
    rstructs::Reactions,
    rids::IDs,
    sid_sets::Vector{Tuple{IDs,Int64}};
    tids::IDs=IDs(),
    write_path::Union{String,Nothing}=nothing,
    allowed_forward::Union{SparseVector{Bool,Int64},Nothing}=nothing,
    allowed_backward::Union{SparseVector{Bool,Int64},Nothing}=nothing)

    rids = remove_rids_not_in_rstructs(rstructs,rids)
    cids = identify_biosystem_compounds(rstructs,rids)
    (R, P) = matrixify_compounds(rstructs,cids,rids)
    t = matrixify_targets(cids,tids)
    # X, Y = Vector{Int}[], Vector{Int}[]
    all_seed_results = Vector{}
    if isnothing(write_path) ## no parallel processing because of the push

        for (i,sid_tup) in enumerate(sid_sets)
            x = matrixify_seeds(sid_tup[1], cids) ## This should be a vector of all 1s of length(cids)
            Vector(x) !== ones(Int,length(cids)) && throw(DomainError("This should be a vector of all 1s of length(cids"))
            X, Y, x = loop_and_remove_seeds(sid_tup[1],cids,x,t,R,P,allowed_forward = allowed_forward, allowed_backward = allowed_backward, skip_seed_indices = sid_tup[2])
            push!(all_seed_results,(x,t,cids,rids,X,Y))
        end

    else ## Use parallel processing

        Threads.@threads for (i,sid_tup) in collect(enumerate(sid_sets))
            x = matrixify_seeds(sid_tup[1], cids) ## This should be a vector of all 1s of length(cids)
            # println(length(x))
            # println(length(cids))
            # println(x)
            ## I'm commenting this out for now because it's throwing even though it shouldn't be based on the above prints
            # x !== ones(Int,length(cids)) && throw(DomainError("This should be a vector of all 1s of length(cids"))
            X, Y, x = loop_and_remove_seeds(sid_tup[1],cids,x,t,R,P,allowed_forward = allowed_forward, allowed_backward = allowed_backward, skip_seed_indices = sid_tup[2])
            simple_write_out(joinpath(write_path,"$i.json"),x,t,cids,rids,X,Y)
        end

    end

    all_seed_results # This will be empty if write_path!=nothing

end

function loop_and_remove_seeds(
    sids::IDs,
    cids::IDs,
    x::SparseVector{Int64,Int64},
    t::SparseVector{Int64,Int64},
    R::SparseMatrixCSC{Int64,Int64},
    P::SparseMatrixCSC{Int64,Int64};
    allowed_forward::Union{SparseVector{Bool,Int64},Nothing}=nothing,
    allowed_backward::Union{SparseVector{Bool,Int64},Nothing}=nothing,
    skip_seed_indices::Int=0)

    tT = sparse(transpose(t))
    sum_t = sum(t)
    ## Run 1 network expansion per seed variation

    X, Y = SparseVector{Int64,Int64}[], SparseVector{Int64,Int64}[]

    for i in seed_indicies(sids, cids)[(skip_seed_indices+1):end]
        x[i] = 0

        X, Y = expandmatrices(R, P, x, allowed_forward = allowed_forward, allowed_backward = allowed_backward) ## global needed to access the variables defined in-loop

        (tT * X[end]) != sum_t && (x[i] = 1) ## This is a short-circuit if statement
        # if (tT * X[end]) != sum_t
        #     x[i] = 1
        # end
    end

    X, Y, x

end

function simple_write_out(
    path::String,
    x::Union{Vector{Int},SparseVector{Int64,Int64}},
    t::Union{Vector{Int},SparseVector{Int64,Int64}},
    cids::IDs,
    rids::IDs,
    X::Union{Vector{Vector{Int}},Vector{SparseVector{Int64,Int64}}},
    Y::Union{Vector{Vector{Int}},Vector{SparseVector{Int64,Int64}}})

    # println("Writing out single network expansion...")
    # if !ispath(path)
    #     mkpath(path)
    # end

    data = Dict()
    data["x"] = x
    data["t"] = t
    data["compounds"] = cids
    data["reactions"] = rids
    data["X"] = X
    data["Y"] = Y

    open(path,"w") do f
        JSON.print(f, data, 2) #indent=2
    end
end

## Change the below to IOs?
# rstructs = readmaster("path/to/master.json")
# rids = readids("path/to/rids.json")
# sids = readids("path/to/sids.json")
# tids = readids("path/to/tids.json")
# path = "path/to/output/netexp/results.json"

# (x, t, cids, X, Y) = expand(rstructs,rids,sids,tids)

end

max_stars_repo_name: hbsmith/BioXP

max_issues_repo_name: hbsmith/BioXP

max_forks_repo_name: hbsmith/BioXP


text:
@testset "code snippets from docs" begin
    @testset "DATA.md - The Network Data Dictionary" begin
        network_data = PowerModels.parse_file("../test/data/case14.m")

        @test length(network_data["bus"]) == 14
        @test length(network_data["branch"]) == 20
    end

    @testset "README.md - Modifying Network Data" begin
        network_data = PowerModels.parse_file("../test/data/case3.m")

        result = run_opf(network_data, ACPPowerModel, IpoptSolver(print_level=0))

        @test result["status"] == :LocalOptimal
        @test isapprox(result["objective"], 5906.88; atol = 1e0)

        network_data["bus"]["3"]["pd"] = 0.0
        network_data["bus"]["3"]["qd"] = 0.0

        result = run_opf(network_data, ACPPowerModel, IpoptSolver(print_level=0))

        @test result["status"] == :LocalOptimal
        @test isapprox(result["objective"], 2937.16; atol = 1e0)
    end

    @testset "README.md - JuMP Model Inspection" begin
        pm = build_generic_model("../test/data/case3.m", ACPPowerModel, PowerModels.post_opf)

        #pretty print the model to the terminal
        #print(pm.model)

        @test MathProgBase.numlinconstr(pm.model) == 8
        @test MathProgBase.numquadconstr(pm.model) == 12
        @test MathProgBase.numconstr(pm.model) - MathProgBase.numlinconstr(pm.model) - MathProgBase.numquadconstr(pm.model) == 12
        @test MathProgBase.numvar(pm.model) == 28

        result = solve_generic_model(pm, IpoptSolver(print_level=0))

        @test result["status"] == :LocalOptimal
        @test isapprox(result["objective"], 5906.88; atol = 1e0)
    end
end

ext: jl


text:
# Code for Video 1.5:

function sum(x, y, z)
  return x + y + z
end

show(sum(2, 3, 4)) # 9
show(sum(2.1, 3.2, 4.3)) # 9.600000000000001
show(sum(true, false, true)) # 2

typeof(sum) # Function

function sum_prod(x, y)
    x + y, x * y
end

# sum("Iris", 7, true)
# ERROR: MethodError: `+` has no method matching +(::ASCIIString, ::Int64)

function sum(x:: Int64, y::Int64, z::Int64)
    return x + y + z
end

function sump(x:: Int64, y::Int64, z::Int64)
    p::Int64 = x * y * z
    return x + y + z + p
end

 sum(x, y, z) = x + y + z

function sumvar(args...)
    sum = 0
    for n in args
        sum += n
    end
    return sum
end

sumvar(2,3,4,100,-9) # 100
sumvar(2,3,4) # 9
sumvar() # 0

function sum(n1, n2, n3=5; n4=-9)
    return n1 + n2 + n3 + n4
end

sum(2, 3) # 1 (= 2 + 3 + 5 -9)
sum(2, 3, 4) # 0  (= 2 + 3 + 4 -9)
sum(2, 3, 4, n4=100) # 109 (= 2 + 3 + 4 + 100)

(x, y, z) -> x + y + z
# (anonymous function)

arr = rand(100)
# 100-element Array{Float64,1}:
#  0.0552819
#  0.319113
#  0.488649
#  ...

arr2 = map(x -> x^2, arr)
# 100-element Array{Float64,1}:
#  0.00305609
#  0.101833
#  0.238778

arr2 = [x^2 for x in arr]
# 100-element Array{Any,1}:
#  0.00305609
#  0.101833
#  0.238778

isodd(n) = (n % 2 != 0)

filter(isodd, round(Int64, rand(100) * 100))
# 50-element Array{Int64,1}:
#  57
#  99
#  23

filter(n -> n %2 != 0, round(Int64, rand(100) * 100))

map(collect(1:100)) do x
        7x^2 - 3.14x + sin(x) + 4
end
# 100-element Array{Float64,1}:
#      8.70147
#     26.6293
#     57.7211
#    102.683
# ...

max_stars_repo_path: Section 1/1_5.jl

max_stars_repo_name: lytemar/Julia-for-Data-Science-Video

max_issues_repo_path: Section 1/1_5.jl

max_issues_repo_name: lytemar/Julia-for-Data-Science-Video

max_forks_repo_path: Section 1/1_5.jl

max_forks_repo_name: lytemar/Julia-for-Data-Science-Video


text:
# Automatically generated using Clang.jl

const STBI_VERSION = 1
const stbi_uc = Cuchar
const stbi_us = UInt16

struct stbi_io_callbacks
    read::Ptr{Cvoid}
    skip::Ptr{Cvoid}
    eof::Ptr{Cvoid}
end

const stbi_write_func = Cvoid

const STBIR_ALPHA_CHANNEL_NONE = -1
const STBIR_FLAG_ALPHA_PREMULTIPLIED = 1 << 0
const STBIR_FLAG_ALPHA_USES_COLORSPACE = 1 << 1
const stbir_uint8 = UInt8
const stbir_uint16 = UInt16
const stbir_uint32 = UInt32

@cenum stbir_edge::UInt32 begin
    STBIR_EDGE_CLAMP = 1
    STBIR_EDGE_REFLECT = 2
    STBIR_EDGE_WRAP = 3
    STBIR_EDGE_ZERO = 4
end

@cenum stbir_filter::UInt32 begin
    STBIR_FILTER_DEFAULT = 0
    STBIR_FILTER_BOX = 1
    STBIR_FILTER_TRIANGLE = 2
    STBIR_FILTER_CUBICBSPLINE = 3
    STBIR_FILTER_CATMULLROM = 4
    STBIR_FILTER_MITCHELL = 5
end

@cenum stbir_colorspace::UInt32 begin
    STBIR_COLORSPACE_LINEAR = 0
    STBIR_COLORSPACE_SRGB = 1
    STBIR_MAX_COLORSPACES = 2
end

@cenum stbir_datatype::UInt32 begin
    STBIR_TYPE_UINT8 = 0
    STBIR_TYPE_UINT16 = 1
    STBIR_TYPE_UINT32 = 2
    STBIR_TYPE_FLOAT = 3
    STBIR_MAX_TYPES = 4
end


text:
module LocalStore
    include("store.jl")
end


text:
# Tests based on the Vector API tutorial found at http://gdal.org/ogr_apitut.html

GU.registerdrivers() do
    GU.read("data/point.geojson") do dataset
        @fact GU.nlayer(dataset) --> 1
        layer = GU.fetchlayer(dataset, 0)
        @fact GU.getname(layer) --> "OGRGeoJSON"
        layerbyname = GU.fetchlayer(dataset, "OGRGeoJSON")
        @fact layerbyname --> layer
        GU.resetreading(layer)

        featuredefn = GU.getlayerdefn(layer)
        @fact GU.nfield(featuredefn) --> 2
        fielddefn = GU.fetchfielddefn(featuredefn, 0)
        @fact GDAL.gettype(fielddefn) --> GDAL.OFTReal
        fielddefn = GU.fetchfielddefn(featuredefn, 1)
        @fact GDAL.gettype(fielddefn) --> GDAL.OFTString

        GU.fetchnextfeature(layer) do feature
            @fact GU.asdouble(feature, 0) --> roughly(2.0)
            @fact GU.asstring(feature, 1) --> "point-a"
        end
        GU.fetchnextfeature(layer) do feature # second feature
            @fact GU.asdouble(feature, 0) --> roughly(3.0)
            @fact GU.asstring(feature, 1) --> "point-b"

            geometry = GU.getgeom(feature)
            @fact GU.getgeomname(geometry) --> "POINT"
            @fact GU.getgeomtype(geometry) --> GDAL.wkbPoint
            @fact GU.nfield(featuredefn) --> 2
            @fact GU.getx(geometry, 0) --> roughly(100.2785)
            @fact GU.gety(geometry, 0) --> roughly(0.0893)
            @fact GU.getpoint(geometry, 0) --> (100.2785,0.0893,0.0)
        end
    end

    pointshapefile = "tmp/point_out"
    GU.create("$pointshapefile.shp", "ESRI Shapefile") do dataset
        layer = GU.createlayer(dataset, "point_out", GDAL.wkbPoint)
        fielddefn = GU.createfielddefn("Name", GDAL.OFTString)
        GU.setwidth(fielddefn, 32)
        GU.createfield(layer, fielddefn, true)
        GU.destroy(fielddefn)

        featuredefn = GU.getlayerdefn(layer)
        @fact GU.getname(featuredefn) --> "point_out"
        GU.createfeature(featuredefn) do feature
            GU.setfield(feature, GU.getfieldindex(feature, "Name"), "myname")
            GU.creategeom(GDAL.wkbPoint) do point
                GU.setpoint(point, 0, 100.123, 0.123)
                GU.setgeom(feature, point)
            end
        end
    end

    rm("$pointshapefile.dbf")
    rm("$pointshapefile.shp")
    rm("$pointshapefile.shx")
end


text:
"""
    VertexCrossover{S, PF, CF}
    VertexCrossover(crossover ;selection=FilterMutationAllowed(), pairgen=default_pairgen)
    VertexCrossover(crossover, deviation::Number; selection=FilterMutationAllowed())

Applies `crossover` to each pair of selected vertices from two `CompGraph`s.

Vertices to select from the first graph is determined by `selection` (default `FilterMutationAllowed`) while `pairgen`
(default [`default_pairgen`](@ref)) determines how to pair the selected vertices with vertices in the second graph.

The default pairing function will try to pair vertices which have similar relative topologial order within their graphs.
For instance, if the first graph has 5 vertices and the second has 10, it will pair vertex 2 from the first graph with
vertex 4 from the second (assuming they are of compatible type). The parameter `deviation` can be used to inject noise
in this process so that the pairing will randomly deviate where the magnitude of `deviation` sets how much and how often.

See also [`crossover`](@ref).
"""
struct VertexCrossover{S, PF, CF} <: AbstractCrossover{CompGraph}
    selection::S
    pairgen::PF
    crossover::CF
end
VertexCrossover(crossover ;selection=FilterMutationAllowed(), pairgen=default_pairgen) = VertexCrossover(selection, pairgen, crossover)
VertexCrossover(crossover, deviation::Number; selection=FilterMutationAllowed()) = VertexCrossover(selection, (vs1,vs2;ind1=1) -> default_pairgen(vs1, vs2, deviation;ind1=ind1), crossover)

(c::VertexCrossover)((g1,g2)::Tuple) = crossover(g1, g2; selection=SelectWithMutation(c.selection, c), pairgen=c.pairgen, crossoverfun=c.crossover)

"""
    CrossoverSwap{F1, F2, F3, S} <: AbstractCrossover{<:AbstractVertex}
    CrossoverSwap(pairgen, mergefun, strategy, selection)
    CrossoverSwap(;pairgen=default_inputs_pairgen, mergefun=default_mergefun, strategy=default_crossoverswap_strategy, selection=FilterMutationAllowed())
    CrossoverSwap(deviation::Number; mergefun=default_mergefun, strategy=default_crossoverswap_strategy, selection=FilterMutationAllowed())

Swap out a part of one graph with a part of another graph, making sure that the graphs do not become connected in the process.

More concretely, swaps a set of consecutive vertices `vs1` set of consecutive vertices `vs2` returning the swapped `v1` and
`v2` respectively if successful or `v1` and `v2` if not.

The last vertex in `vs1` is `v1` and the last vertex of `vs2` is `v2`. The other members of `vs1` and `vs2` are determined
by `pairgen` (default [`default_inputs_pairgen`](@ref)) and `selection` (default `FilterMutationAllowed`).

If a vertex `v` is not capable of having multiple inputs (determined by `singleinput(v) == true`), `vm = mergefun(vi)` where
`vi` is the input to `v` will be used instead of `v` and `v` will be added as the output of `vm` if necessary.

See also [`crossoverswap`](@ref).

Note: High likelyhood of large accuracy degradation after applying this mutation.
"""
struct CrossoverSwap{F1, F2, F3, S} <: AbstractCrossover{AbstractVertex}
    pairgen::F1
    mergefun::F2
    strategy::F3
    selection::S
end
CrossoverSwap(;pairgen=default_inputs_pairgen, mergefun=default_mergefun, strategy=default_crossoverswap_strategy, selection=FilterMutationAllowed()) = CrossoverSwap(pairgen, mergefun, strategy, selection)
CrossoverSwap(deviation::Number; mergefun=default_mergefun, strategy=default_crossoverswap_strategy, selection=FilterMutationAllowed()) = CrossoverSwap((vs1,vs2) -> default_inputs_pairgen(vs1, vs2, deviation), mergefun, strategy, selection)

(c::CrossoverSwap)((v1,v2)::Tuple) = crossoverswap(v1, v2; pairgen = c.pairgen, mergefun=c.mergefun, selection=SelectWithMutation(c.selection, c), strategy=c.strategy)

"""
    crossover(g1::CompGraph, g2::CompGraph; selection=FilterMutationAllowed(), pairgen=default_pairgen, crossoverfun=crossoverswap)

Perform crossover between `g1` and `g2` and return the two children `g1'` and `g2'`.

How the crossover is performed depends on `pairgen` and `crossoverfun` as well as `selection`.

`selection` is used to filter out at which vertices the crossoverpoints may be.

`pairgen` return indices of the potential crossover points to use out of the allowed crossover points. New crossover points will be drawn from pairgen until it returns `nothing`.

`crossoverfun` return the result of the crossover which may be same as inputs depending on implementation.
"""
function crossover(g1::CompGraph, g2::CompGraph; selection=FilterMutationAllowed(), pairgen=default_pairgen, crossoverfun=crossoverswap)

    sel(g) = select(selection, g)

    ninputs1 = length(g1.inputs)
    ninputs2 = length(g2.inputs)
    noutputs1 = length(g1.outputs)
    noutputs2 = length(g2.outputs)

    # pairgen api is very close to the iterator specification. Not sure if things would be easier if it was an iterator instead...
    vs1, vs2 = sel(g1), sel(g2)
    inds = pairgen(vs1, vs2, ind1 = 1)
    while !isnothing(inds)
        ind1, ind2 = inds

        v1, v2 = crossoverfun((vs1[ind1], vs2[ind2]))

        g1 = regraph(v1, ninputs1, noutputs1)
        g2 = regraph(v2, ninputs2, noutputs2)
        vs1, vs2 = sel(g1), sel(g2)

        # Graphs may be different now, so we need to reselect
        inds = pairgen(vs1, vs2, ind1 = ind1+1)
    end
    return g1, g2
end

regraph(v::AbstractVertex) = regraph(all_in_graph(v))
function regraph(vs)
    ins = filter(v -> isempty(inputs(v)), vs)
    outs = filter(v -> isempty(outputs(v)), vs)
    return CompGraph(ins, outs)
end
function regraph(x, ninputs, noutputs)
    g = regraph(x)
    @assert ninputs == length(g.inputs) "Incorrect number of inputs! Expected $ninputs but got $(length(g.inputs))"
    @assert noutputs == length(g.outputs) "Incorrect number of inputs! Expected $noutputs but got $(length(g.outputs))"
    return g
end

"""
    default_pairgen(vs1, vs2, deviation = 0.0; rng=rng_default, compatiblefun = sameactdims, ind1 = rand(rng, eachindex(vs1)))

Return integers `ind1` and `ind2` so that `vs1[ind1]` and `vs2[ind2]` are a suitable pair for crossover.

Input function `compatiblefun(v1,v2)` shall return `true` if `v1` and `v2` can be swapped and is used to determine the set of vertices to select from given.

From the set of compatible vertices in `vs2`, the one which has the smallest relative topologial distance from `vs1[ind2]` is selected. The parameter `devation` can be used to randomly deviate from this where larger magnitude means more deviation.
"""
function default_pairgen(vs1, vs2, deviation = 0.0; rng=rng_default, compatiblefun = sameactdims, ind1 = rand(rng, eachindex(vs1)))
    ind1 > length(vs1) && return nothing

    # Make sure the returned indices map to vertices which are compatible so that output from a convolutional layer does not suddenly end up being input to a dense layer.
    candidate_ind2s = filter(i2 -> compatiblefun(vs1[ind1], vs2[i2]), eachindex(vs2))

    isempty(candidate_ind2s) && return nothing

    order1 = relative_positions(vs1) .+ deviation .* randn(rng, length(vs1))
    order2 = relative_positions(vs2) .+ deviation .* randn(rng, length(vs2))

    ind2 = candidate_ind2s[argmin(abs.(order2[candidate_ind2s] .- order1[ind1]))]
    return ind1, ind2
end

"""
    default_inputs_pairgen(vs1, vs2, args...;kwargs...)

Same as [`default_pairgen`](@ref) except it also ensures that shape changes of feature maps are consistent between the pairs.

Feature map here refers to the shape of inputs to convolutional-type layers (Conv, Pooling) in dimensions other than the batch dimension or the channel dimension.

For example, for 2D convolutions, the arrays may have the shape WxHxCxB where C and B are the channel and batch dimensions respectively. The shape of the feature map in this case is then WxH.

More concretely, if `x` is the shape of a feature maps input to `vin1` and `f1(x)` is the shape of the feature maps output from `vout1` and `f2` describes the same relation between `vin2` and `vout2` then only vertices `vin2'` for which `f1(x) == f2(x) ∀ x` for a selected vertex `vin1` may be returned.

This function assumes that the last vertex in `vs1` and `vs2` are `vout1` and `vout2` respectively.

This prevents issues where graphs become inconsistent due to
1) Feature maps become zero sized
2) Feature maps of different sizes are merged (e.g. concatenated or element wise added)

Note that this is more strict than needed as a change in the feature maps size does not necessarily result in 1) or 2).
"""
function default_inputs_pairgen(vs1, vs2, args...;kwargs...)
    samedimsandshapes(vin1, vin2) = NaiveGAflux.sameactdims(vin1,vin2) && sameoutshape(vin1, last(vs1), vin2, last(vs2))
    return default_pairgen(vs1,vs2, args...;compatiblefun = samedimsandshapes, kwargs...)
end
# We also check all size terminating outputs (perhaps sufficient with one) since transparent layers (incorrectly in this case) just report their inputs dims, meaning that a flattening global pool will seem to have conv output dims
sameactdims(v1, v2) =  _sameactdims(v1, v2) && all(Iterators.product(findterminating(v1, outputs), findterminating(v2, outputs))) do (v1o, v2o)
    _sameactdims(v1o, v2o)
end
_sameactdims(v1, v2) = NaiveNASflux.actdim(v1) == NaiveNASflux.actdim(v2) && NaiveNASflux.actrank(v1) == NaiveNASflux.actrank(v2)

function sameoutshape(vin1, vout1, vin2, vout2)
    tr1 = squashshapes(shapetrace(vout1, vin1))
    tr2 = squashshapes(shapetrace(vout2, vin2))
    return isempty(Δshapediff(tr1,tr2))
end

relative_positions(arr) = collect(eachindex(arr) / length(arr))

"""
    crossoverswap((v1,v2)::Tuple; pairgen=default_inputs_pairgen, selection=FilterMutationAllowed(), kwargs...)

Perform [`crossoverswap!`](@ref) with `v1` and `v2` as output crossover points.

Inputs are selected from the feasible set (as determined by `separablefrom` and `selection`) through the supplied `pairgen` function.

Inputs `v1` and `v2` along with their entire graph are copied before operation is performed and originals are returned if operation is not successful.

Additional keyword arguments will be passed on to [`crossoverswap!`](@ref).
"""
crossoverswap((v1, v2)::Tuple; kwargs...) = crossoverswap(v1, v2; kwargs...)
function crossoverswap(v1, v2; pairgen=default_inputs_pairgen, selection=FilterMutationAllowed(), kwargs...)
    # This is highly annoying: crossoverswap! does multiple remove/create_edge! and is therefore very hard to revert should something go wrong with one of the steps.
    # To mitigate this a backup copy is used. It is however not easy to backup a single vertex as it is connected to all other vertices in the graph, meaning that the whole graph must be copied. Sigh...
    function copyvertex(v)
        g = regraph(v)
        vs = vertices(deepcopy(g))
        return vs[indexin([v], vertices(g))][], filter(vv -> isempty(inputs(vv)), vs)
    end
    v1c, ivs1 = copyvertex(v1)
    v2c, ivs2 = copyvertex(v2)

    vs1 = reverse!(select(selection, separablefrom(v1c, ivs1)))
    vs2 = reverse!(select(selection, separablefrom(v2c, ivs2)))

    # From the two sets of separable vertices, find two matching pairs to use as endpoints in input direction
    # Endpoint in output direction is already determined by sel1 and sel2 and is returned by separablefrom as the first element
    inds = pairgen(vs1, vs2)

    isnothing(inds) && return v1, v2

    ind1, ind2 = inds

    success1, success2 = crossoverswap!(vs1[ind1], vs1[end], vs2[ind2], vs2[end]; kwargs...)

    # Note: Success 1 and 2 are mapped to v2 and v1 respectively as success1 means we successfully inserted v1c into the graph which held v2c and vice versa
    v1ret = success2 ? v1c : v1
    v2ret = success1 ? v2c : v2

    # Note! Flip v1 and v2 for same reason as success1 and 2 are flipped in crossoverswap!
    return v2ret, v1ret
end

"""
    crossoverswap!(v1::AbstractVertex, v2::AbstractVertex) = crossoverswap!(v1,v1,v2,v2)
    crossoverswap!(vin1::AbstractVertex, vout1::AbstractVertex, vin2::AbstractVertex, vout2::AbstractVertex)

Swap vertices `vin1` to `vout1` with `vin2` and `vout2` so that `vin1` to `vin2` is placed in the same position of the graph as `vin2` to `vout2` and vice versa.

Vertices must come from different graphs.

This operation can fail, leaving one or both graphs in a corrupted state where evaluating them results in an error (typically a `DimensionMismatch` error).

Return a tuple `(success1, success2)` where `success1` is true if `vin2` and `vout2` was successfully swapped in to the graph which previously contained `vin1` and `vout1` and vice versa for `success2`.
"""
crossoverswap!(v1::AbstractVertex, v2::AbstractVertex; kwargs...) = crossoverswap!(v1,v1,v2,v2; kwargs...)

function crossoverswap!(vin1::AbstractVertex, vout1::AbstractVertex, vin2::AbstractVertex, vout2::AbstractVertex; mergefun = default_mergefun, strategy=default_crossoverswap_strategy)

    # If vin1 can only take a single input while vin2 has multiple inputs (or vice versa) we need to add a vertex before which merges the inputs
    vin1, vin2 = check_singleinput!(vin1, vin2, mergefun)

    # Beware: ix and ox are not the same type of thing!! Check the strip functions!
    i1, o1 = stripedges!(vin1, vout1)
    i2, o2 = stripedges!(vin2, vout2)

    # success1 mapped to vin2 and vout2 looks backwards, but remember that vin2 and vout2 are the new guys being inserted
    # in everything connected to i1 and o1
    # Returning success status instead of acting as a noop at failure is not very nice, but I could not come up with a way
    # which was 100% to revert a botched attempt and making a backup of vertices before doing any changes is not easy to deal
    # with for the receiver either
    success1 = addinedges!(vin2, i1, strategy)
    success1 &= success1 && addoutedges!(vout2, o1, strategy)

    success2 = addinedges!(vin1, i2, strategy)
    success2 &= success2 && addoutedges!(vout1, o2, strategy)

    return success1, success2
end

function check_singleinput!(v1, v2, mergefun)
    singleinputs = singleinput.([v1, v2])

    needmerge1 = singleinputs[1] && !all(singleinputs)
    needmerge2 = singleinputs[2] && !all(singleinputs)

    function addmerge!(v)
        # vs -> [v] means only add the new vertex between vi and v as vi could have other outputs
        insert!(inputs(v)[], vi -> mergefun(vi), vs -> [v])
        return inputs(v)[]
    end
    v1 = needmerge1 ? addmerge!(v1) : v1
    v2 = needmerge2 ? addmerge!(v2) : v2
    return v1, v2
end

function default_crossoverswap_strategy(utilityfun = NaiveNASlib.defaultutility)
    warnfailalign = FailAlignSizeWarn(msgfun = (vin,vout) -> "Failed to align sizes for vertices $(name(vin)) and $(name(vout)) for crossover. Attempt aborted!")
    alignstrat = TruncateInIndsToValid(WithUtilityFun(utilityfun, AlignNinToNout(;fallback=ΔSizeFailNoOp())))
    return PostAlign(alignstrat, fallback=warnfailalign)
end

"""
    CrossoverSizeDummy(startsize::Int)

Struct whose only purpose is to follow and apply size changes from NaiveNASlib during the crossover operation.

Acts a SizeInvariant operation except for when there are no inputs in which case the latest size seen is given.

This is a bit of a hack as when we do crossover we will end up with size transparent vertices with zero inputs.
NaiveNASlib thinks such things are invalid and does not attempt to handle them in any way.
It turns out that in the context of what we are doing in crossover, the only issue is that we end up with a vertex of
zero size and this means that NaiveNASlib creates zero variables to determine the new size which in turn means that any
size change is impossible. Just one single variable is enough though to create as many new neurons as needed, but we also
need to reflect any size changes made by NaiveNASlib to satisfy its (optional) size validation check.
"""
mutable struct CrossoverSizeDummy
    s::Int
end
NaiveNASlib.nin(s::CrossoverSizeDummy, t::NaiveNASlib.SizeInvariant, v::AbstractVertex) = isempty(inputs(v)) ? [s.s] : nin(identity, t, v)
NaiveNASlib.nout(s::CrossoverSizeDummy, t::NaiveNASlib.SizeInvariant, v::AbstractVertex) = isempty(outputs(v)) ? s.s : nout(identity, t, v)
NaiveNASlib.Δsize!(s::CrossoverSizeDummy, ins::AbstractVector, outs::AbstractVector) = s.s = length(outs)
NaiveNASlib.defaultutility(s::CrossoverSizeDummy) = 0

# Invariant vertex so that removal is trivial.
# Using size absorbing dummies leads to selection of input neurons which dont exist (e.g. take input nr 243 from a layer with 16 inputs).
# Using size stacking dummies leads to neurons being whiped out (e.g. when connecting a 16 neuron output to a 8 neuron input then all 8 are replaced with new neurons).
# We must ensure size is non-zero though, or else NaiveNASlib will not generate any variables for setting the size
dummyvertex(v) = invariantvertex(CrossoverSizeDummy(nout(v)), v; traitdecoration = named("$(name(v)).dummy"))

stripedges!(vin, vout) = stripinedges!(vin) ,stripoutedges!(vout)

function stripinedges!(v)
    # We would like to just remove all input edges using remove_edge!, but this messes up the order of inputs to the outputs of v if v is not the only input
    # Instead we add a  dummyvertex after each input to v and disconnect it (which preserves order as we know it has only one input)
    i = copy(inputs(v))

    for (ind, vi) in enumerate(i)

        # We do this instead of insert! as it is a bit too smart and recreates all edges if any vi == vj for i != j where vi and vj are input vertices to v and this throws addinedges! off.
        dummy = dummyvertex(vi)

        # Manually replace v with dummy as output for vi
        inputs(v)[ind] = dummy
        push!(outputs(dummy), v)
        deleteat!(outputs(vi), findall(vx -> vx == v, outputs(vi)))

        remove_edge!(vi, dummy, strategy = NoSizeChange())
    end
    return i
end

function addinedges!(v, ins, strat = default_crossoverswap_strategy)
    dummies = copy(inputs(v))
    outs = copy(dummies)

    #More outs than ins: We want to connect every input so lets just pad outs with v as this is what we want in the end
    while length(outs) < length(ins)
        push!(outs, v)
    end
    # More outs than ins: This means that v had more inputs in its old graph compared to the vertex it was swapped with.
    # No problem really as we only care about connecting the ins and extra dummies can be left hanging with no input before removal. However, this might leave v with larger nout than it really has if not removed before adding inputs.
    while length(ins) < length(dummies)
        # If we end up here outs and dummies are just copies of each other so we just pop the same element from both.
        pop!(outs)
        vrm = pop!(dummies)
        remove!(vrm, RemoveStrategy(ConnectNone(), NoSizeChange()))
    end
    create_edge_strat = (i == length(ins) ? strat() : NoSizeChange() for i in eachindex(ins))
    success = map((iv, ov, s) -> create_edge!(iv, ov; strategy = s), ins, outs, create_edge_strat) |> all
    remove_dummy_strat = (i == length(dummies) ? strat() : NoSizeChange() for i in eachindex(dummies))
    return success && map((dv,  s)-> remove!(dv, RemoveStrategy(s)), dummies, remove_dummy_strat) |> all
end

function stripoutedges!(v)
    # Similar story as stripinedges to avoid messing with the input order for outputs to v.
    insert!(v, dummyvertex, reverse)
    dummy = outputs(v)[]
    remove_edge!(v, dummy; strategy = NoSizeChange())
    return dummy
end

function addoutedges!(v, dummy, strat = default_crossoverswap_strategy)
    # Perhaps try to avoid large negative size changes here, e.g. by doing a Δnout/Δnin before creating the edge?
    # Should perhaps be baked into strat, but then one probably can't use the same strat for addinedges!
    success = create_edge!(v, dummy, strategy = strat())
    return success && remove!(dummy, RemoveStrategy(strat()))
end

"""
    separablefrom(v)

Return an array of vertices which may be separated from the graph.

More precisely, a connected component which is not connected to the graph can be created if
    1. All output edges from `v` are removed
    2. All input edges from a vertex `v'` in the returned array are removed

The disconncted component has `v'` as first vertex and `v` as last and contains all vertices in between them.

Note that output always contains `v`, i.e it is never empty.
"""
function separablefrom(v, forbidden = AbstractVertex[])
    # Rewrite in a guaranteed to be non-destrucive manner? LightGraphs?
    o = stripoutedges!(v)
    swappable = separablefrom(v, forbidden, AbstractVertex[])
    addoutedges!(v, o, NoSizeChange)
    return swappable
end

function separablefrom(v, forbidden, seen)
    push!(seen, v)
    ins = stripinedges!(v)
    seen_and_dummies = vcat(seen, inputs(v))
    ok = all(vv -> vv in seen_and_dummies && vv ∉ forbidden, all_in_graph(v))
    addinedges!(v, ins, NoSizeChange)
    swappable = mapreduce(vi -> separablefrom(vi, forbidden, seen), vcat, inputs(v), init=[])
    return ok ? vcat(v, swappable) : swappable
end

"""
    OptimizerCrossover{C} <: AbstractCrossover{FluxOptimizer}
    OptimizerCrossover()
    OptimizerCrossover(crossover)

Apply crossover between optimizers.

Type of crossover is determined by `crossover` (default `optimizerswap`) which when given a a tuple of two optimizers will return the result of the crossover operation as a tuple of optimizers.

Designed to be composable with most utility `AbstractMutation`s as well as with itself. For instance, the following seemingly odd construct will swap components of a `Flux.Optimiser` with a probability of `0.2` per component:

`OptimizerCrossover(MutationProbability(OptimizerCrossover(), 0.2))`

Compare with the following which either swaps all components or none:

`MutationProbability(OptimizerCrossover(), 0.2)`
"""
struct OptimizerCrossover{C} <: AbstractCrossover{FluxOptimizer}
    crossover::C
end
OptimizerCrossover() = OptimizerCrossover(optimizerswap)

"""
    LearningRateCrossover()

Return an `OptimizerCrossover` which will swap learning rates between optimizers but not change anything else.

Does not do anything if any of the optimizers don't have a learning rate (e.g. WeightDecay).
"""
LearningRateCrossover() = OptimizerCrossover(learningrateswap)

EitherIs{T} = Union{Tuple{T, Any}, Tuple{Any,T}}

(oc::OptimizerCrossover)(os) = oc.crossover(os)
(oc::OptimizerCrossover)(os::EitherIs{ShieldedOpt}) = os
function (oc::OptimizerCrossover)((o1,o2)::EitherIs{Flux.Optimiser})
    os1,o1re = optiter(o1)
    os2,o2re = optiter(o2)
    res = oc.crossover.(zip(os1,os2))
    os1n = (t[1] for t in res)
    os2n = (t[2] for t in res)
    return o1re(os1n..., os1[length(os2)+1:end]...), o2re(os2n..., os2[length(os1)+1:end]...)
end

optiter(o) = (o,), (os...) -> os[1]
optiter(o::Flux.Optimiser) = Tuple(o.os), (os...) -> Flux.Optimiser(os...)

optimizerswap((o1, o2)::Tuple) = o2,o1

learningrateswap((o1,o2)::Tuple) = (@set o1.eta = learningrate(o2)) , (@set o2.eta = learningrate(o1))
learningrateswap(os::EitherIs{ShieldedOpt}) = os
learningrateswap(os::EitherIs{WeightDecay}) = os


text:
# Script to compare the stacking reoutines in SAC and those
# in Beamforming.jl
#
# Requires Plots in addition to the packages in this project

import Pkg
Pkg.activate(joinpath(@__DIR__, "..", ".."))

using Revise
using DelimitedFiles: writedlm
import Beamforming, Seis, Plots, Seis.Plot
using Plots: plot

# Use UK network data that comes with SAC
t = Seis.sample_data(:array)
t.meta.file = joinpath.(@__DIR__, "stack." .* Seis.channel_code.(t))
Seis.write_sac.(t, t.meta.file)

sacdata = Dict{String,Seis.Trace{Float32,Vector{Float32},String}}()
seisdata = deepcopy(sacdata)

# Time limits of stack, relative to A marker
t1, t2 = -5, 10

# Stacking weights
weights = rand(length(t))
writedlm(joinpath(@__DIR__, "stack_weights.txt"), weights)

for stack_type in ("linear", "nthroot", "phaseweight"), n in 1:2:5
    stack_file = joinpath(@__DIR__, "$(stack_type)")
    if stack_type == "linear"
        n > 1 && continue
        nstring = ""
    else
        stack_file *= "-$n"
        nstring = "$n"
    end
    stack_file *= ".SHZ"

    open(pipeline(`sac`, stdin), "w") do f
        println(f, "sss")
        for (tt, w) in zip(t, weights)
            println(f, "as $(tt.meta.file) weight $w delay $(-tt.picks.A.time)")
        end
        println(f, """
            tw $t1 $t2
            sumstack type $stack_type $nstring
            ws $(stack_file)
            """)
    end

    sacdata[stack_file] = Seis.read_sac(stack_file)
    s = Beamforming.stack(t, (t1, t2), :A,
        method=Symbol(stack_type), n=n)
    s.meta.file = stack_file
    seisdata[stack_file] = s
    @show ≈(Seis.trace(Seis.normalise(s)),
        Seis.trace(Seis.normalise(sacdata[stack_file])), rtol=0.2)
end

sacdata = sort(collect(values(sacdata)), by=x->x.meta.file)
seisdata = sort(collect(values(seisdata)), by=x->x.meta.file)

p = plot(sacdata .|> Seis.normalise, label=last.(splitpath.(sacdata.meta.file)))
for i in 1:length(seisdata)
    s = seisdata[i] |> Seis.normalise
    plot!(p[i], Seis.times(s), Seis.trace(s), line=(:dash, :red))
end
p

max_stars_repo_name: anowacki/Beamforming

max_issues_repo_name: anowacki/Beamforming

max_forks_repo_name: anowacki/Beamforming


text:
function searchforexactheaderstring(header_string, file_strings)
    # find the next "loop_"
    ind = findfirst(xx->header_string==xx, file_strings)
    if typeof(ind) == Int # must not be empty.
        return ind
    end

    return 0
end

function searchforheaderstring(header_string::String, file_strings)
    # find the next "loop_"
    ind = findfirst(xx->occursin(header_string, xx), file_strings)
    if typeof(ind) == Int # must not be empty.
        return ind
    end

    return 0
end

function searchforheaderstring(header_strings::Vector{String}, file_strings)
    # find the next "loop_"
    ind = findfirst(xx->any( occursin(header_strings[i], xx) for i = 1:length(header_strings) ), file_strings)
    if typeof(ind) == Int # must not be empty.
        return ind
    end

    return 0
end

function findblock(file_strings,
                    header_phrase,
                    termination_phrase,
                    file_label::String)
    #
    ind = searchforheaderstring(header_phrase, file_strings)
    if !(ind > 0)
        msg_string = "Error processing $(header_phrase) in $(file_label)."
        println(msg_string)
        return Vector{String}(undef, 0)
    end

    text_buffer = file_strings[ind+1:end]

    # find block termination.
    stop_ind = searchforheaderstring(termination_phrase, text_buffer)
    if !(stop_ind > 0)
        msg_string = "Error processing $(termination_phrase) in $(file_label)."
        println(msg_string)
        return Vector{String}(undef, 0)
    end

    return text_buffer[1:stop_ind-1]
end
# config_path = "/home/roy/MEGAsync/inputs/NMR/configs/cs_config.txt"
# file_strings = readlines(config_path)
# header_phrase = "(-)-Norepinephrine"
# termination_phrase = "end"
# Q = findblock(file_strings, header_phrase, termination_phrase, config_path)


text:
include("datapipeline.jl")
include("densenet.jl")
include("utils.jl")
include("train.jl")
using CUDAnative

#--------------Hyperparameters-------------
model = get_densenet_model(169)
threshold = 0.5
lr = 0.0001
batch_size = 16
which_cat = "XR_HAND"
save_interval = 300
overwrite_save_interval = 3 # Overwrite the model after `x` saves. If `never` overwrite set it to `0`
model_save_path = "../saved_models/model_$(which_cat)" # Donot put the extension it is automatically managed
verbose = 1
epochs = 1
plot_save_path = "../saved_plots/plot_$(which_cat).png"
#------------------------------------------

data_dict = get_batched_images(which_cat, batch_size, path_t1 = "../MURA-v1.1/train_labeled_studies.csv", path_t2 = "../MURA-v1.1/train_image_paths.csv", path_v1 = "../MURA-v1.1/valid_labeled_studies.csv", path_v2 = "../MURA-v1.1/valid_image_paths.csv")

counts = get_count(data_dict)

weight_1 = Dict(i => counts[1][i]/(counts[1][i] + counts[2][i]) for i in ["train", "valid"])

weight_2 = Dict(i => counts[2][i]/(counts[1][i] + counts[2][i]) for i in ["train", "valid"])

loss(x_true, x_pred, cat) = mean(- weight_2[cat] * x_true .* CUDAnative.log.(x_pred) - weight_1[cat] * (1.0 - x_true) .* CUDAnative.log.(1.0 - x_pred))

accuracy(x_true, x_pred, threshold = threshold) = mean((x_pred .>= threshold) .== x_true)

preci(mat) = mat[2, 2] / (mat[2, 2] + mat[1, 2])

recall(mat) = mat[2, 2] / (mat[2, 2] + mat[2, 1])

function f1_score(mat)
  p = preci(mat)
  r = recall(mat)
  2 * p * r / (p + r)
end

cost_metric = Dict("train" => [], "valid" => [])
accuracy_metric = Dict("train" => [], "valid" => [])

opt = ADAM(params(model), lr)

for i in 1:epochs
  info("Starting Epoch $i")
  Flux.Optimise.@interrupts train_model()
  @save model_save_path*"_end_epoch.bson" model
  Flux.Optimise.@interrupts validate_model()
  info("Epoch $i Complete")
end

plot = plot_training(accuracy_metric, cost_metric)

savefig(plot, plot_save_path)


text:
@symbol_func function n_hat(cur_reactor::AbstractReactor)
  cur_nu_n = cur_reactor.nu_n

  cur_n = 0.5

  cur_n *= sqrt( cur_reactor.pi )

  cur_n *= gamma( cur_nu_n + 2.0 )

  cur_n /= gamma( cur_nu_n + 3/2 )

  cur_n *= cur_reactor.n_bar

  cur_n
end


text:
import LinearAlgebra.LAPACK: gebrd!

@testset "test svd $T by gesvd $interface" for T in [Float32, Float64, ComplexF32, ComplexF64], interface in ["CPU", "GPU"]

    # randomly generate a 2 by 2 matrix for testing
    size = 2
    matrixToTest = rand(T, size, size)

    # use the default Linear Algebra lib to calculate the right answer for testing
    right_answer = svd(matrixToTest).S
    S = right_answer

    # to test the GPU interface, one should convert the matrix data to cuda
    if interface == "GPU"
        matrixToTest = cu(matrixToTest)
    end

    # define the job U and job VT, which are required by the MAGMA lib
    jobu = 'A'
    jobvt = 'A'

    # initialize the MAGMA lib, serving as a necessary part before working
    magma_init()

    # call the basic (overloaded) wrapper gesvd! for gesvd subroutines
    result = magma_gesvd!(jobu,jobvt,matrixToTest)

    # in the result, the wanted answer lies in the second position
    s = result[2]

    # finalize the MAGMA lib, serving as a necessary part after working
    magma_finalize()

    # if S is approximately equal to s, we defined then it's alright
    @test S ≈ s

end

@testset "test svd $T by gesdd $interface" for T in [Float32, Float64, ComplexF32, ComplexF64], interface in ["CPU", "GPU"]

    # randomly generate a 2 by 2 matrix for testing
    size = 2
    matrixToTest = rand(T, size, size)

    # use the default Linear Algebra lib to calculate the right answer for testing
    right_answer = svd(matrixToTest).S
    S = right_answer

    # to test the GPU interface, one should convert the matrix data to cuda
    if interface == "GPU"
        matrixToTest = cu(matrixToTest)
    end

    # define the job MAGMA, which are required by the MAGMA lib
    job_magma = 'A'

    # initialize the MAGMA lib, serving as a necessary part before working
    magma_init()

    # call the basic (overloaded) wrapper gesdd! for gesvd subroutines
    result = magma_gesdd!(job_magma,matrixToTest)

    # in the result, the wanted answer lies in the second position
    s = result[2]

    # finalize the MAGMA lib, serving as a necessary part after working
    magma_finalize()

    # if S is approximately equal to s, we defined then it's alright
    @test S ≈ s

end

@testset "test svd $T by gebrd $interface" for T in [Float32, Float64, ComplexF32, ComplexF64], interface in ["CPU", "GPU"]

    # randomly generate a 2 by 2 matrix for testing
    size = 3
    matrixToTest = rand(T, size, size)
    matrixToTest_copy = copy(matrixToTest)

    # use the default Linear Algebra lib to calculate the right answer for testing
    right_answer = gebrd!(matrixToTest)

    # to test the GPU interface, one should convert the matrix data to cuda
    if interface == "GPU"
        matrixToTest_copy = cu(matrixToTest_copy)
    end

    # initialize the MAGMA lib, serving as a necessary part before working
    magma_init()

    # call the basic (overloaded) wrapper gesdd! for gesvd subroutines
    result = magma_gebrd!(matrixToTest_copy)

    # finalize the MAGMA lib, serving as a necessary part after working
    magma_finalize()

    # if S is approximately equal to s, we defined then it's alright
    for index in 1:length(result)
        result_piece = result[index]
        right_piece = right_answer[index]
        if index == 3
            right_piece = right_piece[1:length(right_piece) - 1]
        end
        @test result_piece ≈ right_piece
    end
end


text:
# ---
# title: 633. Sum of Square Numbers
# id: problem633
# author: Indigo
# date: 2021-06-27
# difficulty: Medium
# categories: Math
# link: <https://leetcode.com/problems/sum-of-square-numbers/description/>
# hidden: true
# ---
#
# Given a non-negative integer `c`, decide whether there're two integers `a` and
# `b` such that `a2 + b2 = c`.
#
#
#
# **Example 1:**
#
#
#
#     Input: c = 5
#     Output: true
#     Explanation: 1 * 1 + 2 * 2 = 5
#
#
# **Example 2:**
#
#
#
#     Input: c = 3
#     Output: false
#
#
# **Example 3:**
#
#
#
#     Input: c = 4
#     Output: true
#
#
# **Example 4:**
#
#
#
#     Input: c = 2
#     Output: true
#
#
# **Example 5:**
#
#
#
#     Input: c = 1
#     Output: true
#
#
#
#
# **Constraints:**
#
#   * `0 <= c <= 231 - 1`
#
#
## @lc code=start
using LeetCode

function judge_square_sum(n::Int)
    upper = isqrt(n)
    i = 0
    while i <= upper
        ss = i ^ 2 + upper ^ 2
        ss == n && return true
        ss > n ? (upper -= 1) : (i += 1)
    end
    return false
end
## @lc code=end


text:
function getdescriptors(descriptors)

    for i = length(descriptors):-1:1
        if descriptors[i] === nothing
            deleteat!(descriptors, i)
        end
    end

    n = length(descriptors)
    emdescriptors = Array{Any}(nothing, n)
    mldescriptors = nothing;

    n = 0
    for i = 1:length(descriptors)
        if typeof(descriptors[i]) == Potential.PotentialStruct
            n = n + 1
            emdescriptors[n] = descriptors[i]
        else
            mldescriptors = descriptors[i]
        end
    end

    if n>0
        if length(descriptors) > n
            for i = length(descriptors):-1:(n+1)
                deleteat!(emdescriptors, i)
            end
        end
    else
        emdescriptors = []
    end

    return emdescriptors, mldescriptors
end

function setcutoff(emdescriptors, mldescriptors, eta)
    if length(eta)>0
        for i = 1:length(emdescriptors)
            if i==1
                emdescriptors[i].rcut = eta[1] # long-range empirical potentials
            else
                emdescriptors[i].rcut = eta[1] # short-range empirical potentials
            end
        end
        if mldescriptors !== nothing
            if lowercase(mldescriptors.name) == "snap"
                if length(emdescriptors) == 0
                    mldescriptors.rcutfac = eta[1]
                else
                    mldescriptors.rcutfac = eta[1]
                end
            end
        end
    end
    return emdescriptors, mldescriptors
end

function maespace(traindata, validdata, emdescriptors, mldescriptors, etaspace, kappa, pbc=nothing, a=nothing, b=nothing, c=nothing, method=nothing, normalizeenergy=nothing)

    etamin = etaspace[:,1]
    etamax = etaspace[:,2]
    N = Int64.(etaspace[:,3])
    eta = tensornodes(etamin, etamax, N)

    # read data
    trainconfig = Array{Any}(undef, length(traindata))
    for i = 1:length(traindata)
        trainconfig[i], indices = Preprocessing.readconfigdata(traindata[i], pbc, a, b, c)
    end
    validconfig = Array{Any}(undef, length(validdata))
    for i = 1:length(validdata)
        validconfig[i], indices = Preprocessing.readconfigdata(validdata[i], pbc, a, b, c)
    end

    neta = size(eta,1)
    eme = -1.0
    fme = -1.0
    sme = -1.0
    for n = 1:neta
        emdescriptors, mldescriptors = setcutoff(emdescriptors, mldescriptors, eta[n,:])
        ce, cf, cef, cefs = linearfit(trainconfig, emdescriptors, mldescriptors, eta[n,:], kappa, pbc, a, b, c, method, normalizeenergy)
        emae, fmae, smae = validation(validconfig, emdescriptors, mldescriptors, ce, cf, cef, cefs, eta[n,:], kappa, pbc, a, b, c)

        if n == 1
            sz = size(emae)
            eme = -ones(neta, sz[1], sz[2])
            fme = -ones(neta, sz[1], sz[2])
            sme = -ones(neta, sz[1], sz[2])
        end
        eme[n,:,:] = emae
        fme[n,:,:] = fmae
        sme[n,:,:] = smae

        emae = round.(emae; digits=9)
        fmae = round.(fmae; digits=9)
        smae = round.(smae; digits=9)
        print("            Energy Fitting      Force Fitting      Energy+Force Fitting      Energy+Force+Stress Fitting \n");
        print("Energy MAE:   " * string(emae[1]) * "         " * string(emae[2]) * "         " * string(emae[3]) * "           " * string(emae[4]) * "\n");
        print("Force MAE:    " * string(fmae[1]) * "         " * string(fmae[2]) * "         " * string(fmae[3]) * "           " * string(fmae[4]) * "\n");
        print("Stress MAE:   " * string(smae[1]) * "         " * string(smae[2]) * "         " * string(smae[3]) * "           " * string(smae[4]) * "\n");
    end

    return eta, eme, fme, sme
end

function optimize(traindata, testdata, descriptors, etaspace, kappa, weightouter, lossfunc, pbc=nothing, a=nothing, b=nothing, c=nothing, method=nothing, normalizeenergy=nothing)

    # getdescriptors
    emdescriptors, mldescriptors = getdescriptors(descriptors)

    # compute MAE on the eta space grid
    etapts, emae, fmae, smae = maespace(traindata, testdata, emdescriptors, mldescriptors,
                etaspace, kappa, pbc, a, b, c, method, normalizeenergy)

    # compute the loss function on the eta space grid
    if lossfunc == "energy"
        f = weightouter[1]*emae[:,1,end]
    elseif lossfunc == "force"
        f = weightouter[2]*emae[:,2,end]
    elseif lossfunc == "energyforce"
        f = weightouter[1]*emae[:,3,end] + weightouter[2]*fmae[:,3,end]
    elseif lossfunc == "energyforcestress"
        f = weightouter[1]*emae[:,4,end] + weightouter[2]*fmae[:,4,end] + weightouter[3]*smae[:,4,end]
    else
        error("lossfunc is not valid. It must be energy, force, energyforce, or energyforcestress")
    end

    # interpolate the loss function on the eta space grid
    pc = mpolyinterp(etapts, f, etaspace)

    # optimize the loss function on the eta space
    eta, fmin, iter  = gradientdescent(pc, etaspace)

    # set cut-off radius
    emdescriptors, mldescriptors = setcutoff(emdescriptors, mldescriptors, eta)

    # compute linear coefficients
    ce, cf, cef, cefs = linearfit(traindata, emdescriptors, mldescriptors, eta, kappa, pbc, a, b, c, method, normalizeenergy)

    if lossfunc == "energy"
        coeff = ce
    elseif lossfunc == "force"
        coeff = cf
    elseif lossfunc == "energyforce"
        coeff = cef
    elseif lossfunc == "energyforcestress"
        coeff = cefs
    end

    return eta, coeff, fmin, iter, pc, f, etapts, emae, fmae, smae
end

function validate(data, descriptors, coeff, eta, kappa, pbc=nothing, a=nothing, b=nothing, c=nothing)

    # getdescriptors
    emdescriptors, mldescriptors = getdescriptors(descriptors)

    # set cut-off radius
    # if eta !== nothing
    #     emdescriptors, mldescriptors = setcutoff(emdescriptors, mldescriptors, eta)
    # end

    normalizeenergy = 0
    n = length(data)
    emae = -ones(Float64,n)
    fmae = -ones(Float64,n)
    smae = -ones(Float64,n)
    szbe = zeros(Int64, n)
    szbf = zeros(Int64, n)
    szbs = zeros(Int64, n)
    for i = 1:n
        if typeof(data[i]) == Preprocessing.DataStruct
            config, indices = Preprocessing.readconfigdata(data[i], pbc, a, b, c)
        else
            config = data[i]
        end
        valid = Array(1:config.nconfigs)

        Aei, Afi, Asi, bei, bfi, bsi, nfi = Potential.emldescriptors(config, valid, eta, kappa, normalizeenergy, emdescriptors, mldescriptors)
        szbe[i] = length(bei)
        szbf[i] = length(bfi)
        szbs[i] = length(bsi)

        emae[i] = calmae(Aei, bei, coeff)
        fmae[i] = calmae(Afi, bfi, coeff)
        smae[i] = calmae(Asi, bsi, coeff)
    end

    eme = sum(emae.*szbe)/sum(szbe)
    fme = sum(fmae.*szbf)/sum(szbf)
    sme = sum(smae.*szbs)/sum(szbs)

    return eme, fme, sme, emae, fmae, smae, szbe, szbf, szbs
end

function setcutoff2(descriptors, eta)
    if length(eta)>0
        for i = 1:length(descriptors)
            if typeof(descriptors[i]) == Potential.PotentialStruct
                descriptors[i].rcut = eta[1]
            elseif typeof(descriptors[i]) == Potential.SnaStruct
                descriptors[i].rcutfac = eta[1]
                descriptors[3].rcutmax = eta[1]
            end
        end
    end
    return descriptors
end

function validate2(data, descriptors, potential, optim, coeff, pbc=nothing, a=nothing, b=nothing, c=nothing)

    data = deletenothing(data)
    descriptors = deletenothing(descriptors)
    potential = deletenothing(potential)

    lossfunc, method, normalizeenergy, normalizestress, eta, kappa,
        weightouter, etaspace, kappalist = getoptim(optim)

    n = length(data)
    emae = -ones(Float64,n)
    fmae = -ones(Float64,n)
    smae = -ones(Float64,n)
    ermse = -ones(Float64,n)
    frmse = -ones(Float64,n)
    srmse = -ones(Float64,n)
    ersq = -ones(Float64,n)
    frsq = -ones(Float64,n)
    srsq = -ones(Float64,n)
    essr = -ones(Float64,n)
    fssr = -ones(Float64,n)
    sssr = -ones(Float64,n)
    eavg = -ones(Float64,n)
    favg = -ones(Float64,n)
    savg = -ones(Float64,n)
    szbe = zeros(Int64, n)
    szbf = zeros(Int64, n)
    szbs = zeros(Int64, n)
    be = []
    bf = []
    bs = []
    for i = 1:n
        if typeof(data[i]) == Preprocessing.DataStruct
            config, indices = Preprocessing.readconfigdata(data[i], pbc, a, b, c)
        else
            config = data[i]
        end
        valid = Array(1:config.nconfigs)

        Aei, Afi, Asi, bei, bfi, bsi, nfi = Potential.emldescriptors2(config, valid, eta, kappa,
                        normalizeenergy, normalizestress, descriptors, potential)

        szbe[i] = length(bei)
        szbf[i] = length(bfi)
        szbs[i] = length(bsi)

        emae[i], ermse[i], ersq[i], essr[i], eavg[i] = calerrors(Aei, bei, coeff)
        fmae[i], frmse[i], frsq[i], fssr[i], favg[i] = calerrors(Afi, bfi, coeff)
        smae[i], srmse[i], srsq[i], sssr[i], savg[i] = calerrors(Asi, bsi, coeff)

        if szbe[i] > 0
            be = [be; bei]
        end
        if szbf[i] > 0
            bf = [bf; bfi]
        end
        if szbs[i] > 0
            bs = [bs; bsi]
        end
    end

    energyerrors = -ones((n+1),3)
    forceerrors = -ones((n+1),3)
    stresserrors = -ones((n+1),3)

    if sum(szbe) > 0
        ssr = sum(essr)
        avg = sum(eavg.*szbe)/sum(szbe)
        rsq = 1.0 - ssr/sum((be .- avg).^2);
        energyerrors[1,1] = sum(emae.*szbe)/sum(szbe)  # MAE
        energyerrors[1,2] = sqrt(sum((ermse.^2).*szbe)/sum(szbe)) # RMSE
        energyerrors[1,3] = rsq
        energyerrors[2:end,:] = [emae ermse ersq]
    end

    if sum(szbf) > 0
        ssr = sum(fssr)
        avg = sum(favg.*szbf)/sum(szbf)
        rsq = 1.0 - ssr/sum((bf .- avg).^2);
        forceerrors[1,1] = sum(fmae.*szbf)/sum(szbf)
        forceerrors[1,2] =  sqrt(sum((frmse.^2).*szbf)/sum(szbf))
        forceerrors[1,3] = rsq
        forceerrors[2:end,:] = [fmae frmse frsq]
    end

    if sum(szbs) > 0
        ssr = sum(sssr)
        avg = sum(savg.*szbs)/sum(szbs)
        rsq = 1.0 - ssr/sum((bs .- avg).^2);
        stresserrors[1,1] = sum(smae.*szbs)/sum(szbs)
        stresserrors[1,2] = sqrt(sum((srmse.^2).*szbs)/sum(szbs))
        stresserrors[1,3] = rsq
        stresserrors[2:end,:] = [smae srmse srsq]
    end

    return energyerrors, forceerrors, stresserrors
end

function maespace2(traindata, validdata, descriptors, potential, optim, pbc=nothing, a=nothing, b=nothing, c=nothing)

    # generate the eta space grid
    etaspace = optim.etaspace
    etamin = etaspace[:,1]
    etamax = etaspace[:,2]
    N = Int64.(etaspace[:,3])
    eta = tensornodes(etamin, etamax, N)

    # read data
    trainconfig = Array{Any}(undef, length(traindata))
    for i = 1:length(traindata)
        trainconfig[i], indices = Preprocessing.readconfigdata(traindata[i], pbc, a, b, c)
    end
    validconfig = Array{Any}(undef, length(validdata))
    for i = 1:length(validdata)
        validconfig[i], indices = Preprocessing.readconfigdata(validdata[i], pbc, a, b, c)
    end

    neta = size(eta,1)
    emae = -ones(neta)
    fmae = -ones(neta)
    smae = -ones(neta)
    for n = 1:neta
        optim.eta = eta[n,:]
        descriptors = setcutoff2(descriptors, eta[n,:])
        coeff,~ = linearfit2(trainconfig, descriptors, potential, optim, pbc, a, b, c)
        eerr, ferr, serr = validate2(validconfig, descriptors, potential, optim, coeff, pbc, a, b, c)
        emae[n] = eerr[1]
        fmae[n] = ferr[1]
        smae[n] = serr[1]
        a1 = num2string(emae[n], 16);
        a2 = num2string(fmae[n], 16);
        a3 = num2string(smae[n], 16);
        mystr = "$n  "
        for j = 1:length(eta[n,:])
            a = num2string(eta[n,j], 16)
            mystr = mystr * a * "  "
        end
        print(mystr * a1 * "  " * a2 * "  " * a3 * " \n")
    end

    return eta, emae, fmae, smae
end

function optimize2(traindata, testdata, descriptors, potential, optim, pbc=nothing, a=nothing, b=nothing, c=nothing)

    traindata = deletenothing(traindata)
    testdata = deletenothing(testdata)
    descriptors = deletenothing(descriptors)
    potential = deletenothing(potential)

    # compute errors on the eta space grid
    etapts, emae, fmae, smae = maespace2(traindata, testdata, descriptors, potential, optim, pbc, a, b, c)

    lossfunc = optim.lossfunc
    weightouter = optim.weightouter
    etaspace = optim.etaspace

    # compute the loss function on the eta space grid
    if lossfunc == "energy"
        f = weightouter[1]*emae
    elseif lossfunc == "force"
        f = weightouter[2]*emae
    elseif lossfunc == "energyforce"
        f = weightouter[1]*emae + weightouter[2]*fmae
    elseif lossfunc == "energyforcestress"
        f = weightouter[1]*emae + weightouter[2]*fmae + weightouter[3]*smae
    else
        error("lossfunc is not valid. It must be energy, force, energyforce, or energyforcestress")
    end

    # interpolate the loss function on the eta space grid
    pc = mpolyinterp(etapts, f, etaspace)

    # optimize the loss function on the eta space
    eta, fmin, iter  = gradientdescent(pc, etaspace)
    optim.eta = eta

    # set cut-off radius
    descriptors = setcutoff2(descriptors, eta)

    # compute coefficient
    coeff,~ = linearfit2(traindata, descriptors, potential, optim, pbc, a, b, c)

    return eta, coeff, fmin, iter, pc, f, etapts, emae, fmae, smae
end

max_stars_repo_name: cesmix-mit/MDP

max_issues_repo_name: cesmix-mit/MDP

max_forks_repo_name: cesmix-mit/MDP


text:
module LazyCat

using RecursiveArrayTools
using Base: @propagate_inbounds

import Base: getindex, setindex!, size

export lazy_cat

function lazy_cat(arrs::AbstractArray{T}...; dim::Int=1) where T
    if length(arrs) == 1
        arrs[1]
    else
        # Pad arrays with smaller size out to the size of the largest dimension.
        N = maximum(ndims.(arrs))
        LazyCatArray{dim}(
                          (arr -> ndims(arr) < N ?
                           view(arr, (Colon() for i = 1:ndims(arr))...,
                                (1:1 for i = ndims(arr)+1:N)...) :
                           arr).(arrs)...
                         )
    end
end

struct NotImplementedError <: Exception
    msg::String
end

struct LazyCatArray{T, N, D, RANGES, ARR} <: AbstractArray{T, N}
    members::ARR

    function LazyCatArray{D}(members::AbstractArray{T, N}...) where {D, T, N}
        if !(D <= N && D > 0)
            throw(NotImplementedError("LazyCatArray{D}: Not implemented for D > dim of member arrays"))
        end

        members = [members...]
        ranges = UnitRange[]

        for i ∈ 1:length(members)
            for d in 1:N
                if d != D
                    i != length(members) && (size(members[i], d) != size(members[i+1], d)) &&
                        throw(ErrorException("LazyCatArray{$D}: size in dimension $d does not match"))
                end
            end
            if isempty(ranges)
                push!(ranges, 1:size(members[i], D))
            else
                push!(ranges, ranges[end][end]+1:ranges[end][end]+size(members[i], D))
            end
        end

        arr = ArrayPartition(members...)
        new{T, N, D, (ranges...,), typeof(arr)}(arr)
    end
end

# Binary search for the proper array.
@propagate_inbounds function get_arrindex(ranges::NTuple{N, UnitRange},
                                          inds::NTuple{M, Int}, i) where {M, N}
    if M == 1
        i ∈ ranges[inds[1]] ? inds[1] : 0
    elseif M == 0
        0
    else
        half1, half2 = divide_inds(inds)
        if ranges[half2[1]][1] > i
            get_arrindex(ranges, half1, i)
        else
            get_arrindex(ranges, half2, i)
        end
    end
end

@generated function indrange(::Val{N}) where N
    exprs = [:($i) for i in 1:N]
    :(tuple($(exprs...)))
end

@propagate_inbounds function get_arrindex(ranges::NTuple{N, UnitRange}, i) where N
    get_arrindex(ranges, indrange(Val(N)), i)
end

@generated function divide_inds(INDS::NTuple{N, Int}) where N
    exprs1 = [:(INDS[$i]) for i in  1:N÷2]
    exprs2 = [:(INDS[$i]) for i in N÷2+1:N]
    quote
        tuple($(exprs1...)), tuple($(exprs2...))
    end
end

@generated function get_subarr_index(::LazyCatArray{T, N, D, RANGES}, I::NTuple{N, Int},
                                     arrindex) where {T, N, D, RANGES}
    sub_masks = ((
                 ((i == D ? RANGES[j][1] : 0 for i ∈ 1:N)...,)
                 for j ∈ 1:length(RANGES))...,
                )
    add_mask = ( (i == D ? 1 : 0 for i ∈ 1:N)..., )

    quote
        I .- $(sub_masks)[arrindex] .+ $add_mask
    end
end

@propagate_inbounds function getindex(A::LazyCatArray{T, N, D, RANGES},
                                      I::Vararg{Int, N}) where {T, N, D, RANGES}
    arrindex = get_arrindex(RANGES, I[D])

    @boundscheck begin
        if arrindex == 0
            throw(BoundsError(A, I))
        end
    end
    I = get_subarr_index(A, I, arrindex)
    A.members[arrindex, I...]
end

@propagate_inbounds function setindex!(A::LazyCatArray{T, N, D, RANGES}, v,
                                       I::Vararg{Int, N}) where {T, N, D, RANGES}
    arrindex = get_arrindex(RANGES, I[D])
    @boundscheck begin
        if arrindex == 0
            throw(BoundsError(A, I))
        end
    end

    I = get_subarr_index(A, I, arrindex)
    A.members[arrindex, I...] = v
end

@generated function size(A::LazyCatArray{T, N, D, RANGES}) where {T, N, D, RANGES}
    mul_mask = ( (i == D ? 0 : 1 for i ∈ 1:N)..., )
    add_mask = ( (i == D ? sum(length.(RANGES)) : 0 for i ∈ 1:N)..., )
    :(size(A.members.x[1]) .* $mul_mask .+ $add_mask)
end

end #module

max_stars_repo_licenses: MIT Unlicense

max_issues_repo_licenses: MIT Unlicense

max_forks_repo_licenses: MIT Unlicense


text:
"""
    mapgrowth(model; series::AbstractGeoSeries, tspan::AbstractRange)
    mapgrowth(layers...; kwargs...)

Combine growth rates accross layers and subperiods for all required periods.

## Arguments

- `model`: A `Model`, NamedTuple of `Model`,  or `Tuple` of `Layer` components.

`Layer`s can also be passed in as separate arguments. A `NamedTuple` of models
will all be run from the same data source, reducing load time.

## Keyword Arguments

- `series`: any AbstractGeoSeries from [GeoData.jl](http://github.com/rafaqz/GeoData.jl)
- `tspan`: `AbstractRange` for the timespan to run the layers for.
  This will become the `index` values of the output `GeoArray` time-dimension `Ti`.

The output is a GeoArray with the same dimensions as the passed in stack layers, and a Time
dimension with a length of `nperiods`.
"""
mapgrowth(layers::Layer...; kw...) = mapgrowth(layers; kw...)
mapgrowth(layers::Tuple{<:Layer,Vararg}; kw...) = mapgrowth(Model(layers); kw...)
mapgrowth(model::Model; kw...) = mapgrowth((model,); kw...)[1]
mapgrowth(m1::Model, ms::Model...; kw...) = mapgrowth((m1, ms...); kw...)
function mapgrowth(models::Union{Tuple{<:Model,Vararg},NamedTuple{<:Any,<:Tuple{<:Model,Vararg}}};
    series::AbstractGeoSeries, tspan::AbstractRange, arraytype=Array, initval=nothing
)
    models = map(stripparams ∘ parent, models)
    period = step(tspan); nperiods = length(tspan)
    startdate, enddate = first(tspan), last(tspan)
    required_keys = Tuple(union(map(keys ∘ _astuple, models)...))

    # Copy only the required keys to a memory-backed stack
    stack = deepcopy(GeoStack(first(series); keys=required_keys))

    A = first(values(stack))
    initval === nothing && (initval = zero(eltype(A)))
    missingval = eltype(A)(NaN)

    # Replace false with NaN
    mask = map(x -> x ? eltype(A)(x) : missingval, boolmask(A)) |> parent |> arraytype
    stackbuffer = GeoData.modify(arraytype, stack)

    # Make a 3 dimensional GeoArray for output, adding the time dimension
    ti = Ti(tspan; mode=Sampled(Ordered(), Regular(period), Intervals(Start())))
    outdims = (dims(A)..., ti)
    outA = arraytype(fill(initval, size(A)..., nperiods))
    outputs = map(models) do m
        GeoArray(deepcopy(outA), outdims; name=:growthrate, missingval=missingval)
    end

    runperiods!(outputs, stackbuffer, series, mask, models, tspan)

    # Return a GeoArray wrapping a regular Array, not arraytype
    map(o -> GeoData.modify(Array, o), outputs)
end

function runperiods!(outputs, stackbuffer, series, mask, models, tspan)
    period = step(tspan); nperiods = length(tspan)
    println("Running for $(1:nperiods)")
    for p in 1:nperiods
        n = 0
        periodstart = tspan[p]
        periodend = periodstart + period
        println("\n", "Processing period between: $periodstart and $periodend")

        # We don't use `Between` as it might unintentionally cut off the
        # last time if it partially extends beyond the period.
        # So we just work with time as Points using `Where`.
        subseries = series[Ti(Where(t -> t >= periodstart && t < periodend))]
        for t in 1:size(subseries, Ti)
            println("    ", val(dims(subseries, Ti))[t])
            # Copy the arrays we need from disk to the buffer stack
            copy!(stackbuffer, subseries[t])
            map(outputs, models, keys(outputs)) do output, model, key
                # For some reason now this is broken with DD getindex, view is a workaround
                parent(view(output, Ti(p))) .+= combinelayers(model, stackbuffer)
            end
            n += 1
        end
        if n > 0
            map(outputs) do output
                parent(view(output, Ti(p))) .*= mask ./ n
            end
        else
            @warn ("No files found for the $period period starting $periodstart")
            map(outputs) do output
                parent(view(output, Ti(p))) .*= mask
            end
        end
    end
end

@inline combinelayers(layer, stackbuffer) = combinelayers((layer,), stackbuffer)
@inline combinelayers(layers::Tuple, stackbuffer::AbstractGeoStack) =
    conditionalrate.(Ref(first(layers)), parent(stackbuffer[keys(first(layers))])) .+ combinelayers(tail(layers), stackbuffer)
@inline combinelayers(layers::Tuple{T}, stackbuffer::AbstractGeoStack) where T =
    conditionalrate.(Ref(first(layers)), parent(stackbuffer[keys(first(layers))]))
@inline combinelayers(layers::Tuple, vals::NamedTuple) =
    conditionalrate.(Ref(first(layers)), vals[keys(first(layers))]) .+ combinelayers(tail(layers), vals)
@inline combinelayers(layers::Tuple{T}, vals::NamedTuple) where T =
    conditionalrate.(Ref(first(layers)), vals[keys(first(layers))])

# Special-case Celcius: some data sets are in °C, but layers
# Can never use °C, so we convert it.
maybetoK(val) = unit(val) == u"°C" ? val |> K : val

@inline conditionalrate(l::Layer, val) = conditionalrate(model(l), maybetoK(unit(l) * val))
@inline conditionalrate(l::Layer, val::Unitful.AbstractQuantity) =
    conditionalrate(model(l), maybetoK(val))
@inline conditionalrate(model::RateModel, val) =
    condition(model, val) ? rate(model, val) : zero(rate(model, val))

_astuple(x::Tuple) = x
_astuple(x) = (x,)


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

mutable struct IPAddressAvailabilityResult <: SwaggerModel
    available::Any # spec type: Union{ Nothing, Bool } # spec name: available
    availableIPAddresses::Any # spec type: Union{ Nothing, Vector{String} } # spec name: availableIPAddresses

    function IPAddressAvailabilityResult(;available=nothing, availableIPAddresses=nothing)
        o = new()
        validate_property(IPAddressAvailabilityResult, Symbol("available"), available)
        setfield!(o, Symbol("available"), available)
        validate_property(IPAddressAvailabilityResult, Symbol("availableIPAddresses"), availableIPAddresses)
        setfield!(o, Symbol("availableIPAddresses"), availableIPAddresses)
        o
    end
end # type IPAddressAvailabilityResult

const _property_map_IPAddressAvailabilityResult = Dict{Symbol,Symbol}(Symbol("available")=>Symbol("available"), Symbol("availableIPAddresses")=>Symbol("availableIPAddresses"))
const _property_types_IPAddressAvailabilityResult = Dict{Symbol,String}(Symbol("available")=>"Bool", Symbol("availableIPAddresses")=>"Vector{String}")
Base.propertynames(::Type{ IPAddressAvailabilityResult }) = collect(keys(_property_map_IPAddressAvailabilityResult))
Swagger.property_type(::Type{ IPAddressAvailabilityResult }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IPAddressAvailabilityResult[name]))}
Swagger.field_name(::Type{ IPAddressAvailabilityResult }, property_name::Symbol) =  _property_map_IPAddressAvailabilityResult[property_name]

function check_required(o::IPAddressAvailabilityResult)
    true
end

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


text:
module SpringCollab2020InvisibleLightSource

using ..Ahorn, Maple

@mapdef Entity "SpringCollab2020/invisibleLightSource" InvisibleLightSource(x::Integer, y::Integer, alpha::Number=1.0, radius::Number=48.0, startFade::Number=24.0, endFade::Number=48.0, color::String="White")

const colors = sort(collect(keys(Ahorn.XNAColors.colors)))

const placements = Ahorn.PlacementDict(
	"Light Source (Spring Collab 2020)" => Ahorn.EntityPlacement(
		InvisibleLightSource
	)
)

Ahorn.editingOptions(entity::InvisibleLightSource) = Dict{String,Any}(
	"color" => colors
)

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

	return Ahorn.Rectangle(x - 4, y - 4, 7, 8)
end

function Ahorn.renderAbs(ctx::Ahorn.Cairo.CairoContext, entity::InvisibleLightSource, room::Maple.Room)
	x, y = Ahorn.position(entity)
	sprite = Ahorn.getTextureSprite("objects/hanginglamp", "Gameplay")

	Ahorn.drawImage(ctx, sprite, x - 4, y - 4, 0, 16, 7, 8, alpha=0.7)
end

end


text:
@kwargs function kron!(C::GBMatOrTrans,
    A::GBMatOrTrans,
    B::GBMatOrTrans,
    op::MonoidSemiringOrBinary = BinaryOps.TIMES
)
    op = getoperator(op, eltype(C))
    A, desc, B = _handletranspose(A, desc, B)
    if op isa Union{Abstract_GrB_BinaryOp, libgb.GrB_BinaryOp}
        libgb.GxB_kron(C, mask, accum, op, A, B, desc)
    elseif op isa Union{Abstract_GrB_Monoid, libgb.GrB_Monoid}
        libgb.GrB_Matrix_kronecker_Monoid(C, mask, accum, op, A, B, desc)
    elseif op isa Union{Abstract_GrB_Semiring, libgb.GrB_Semiring}
        libgb.GrB_Matrix_kronecker_Semiring(C, mask, accum, op, A, B, desc)
    end
end

@kwargs function kron(
    A::GBMatOrTrans,
    B::GBMatOrTrans,
    op::MonoidSemiringOrBinary = BinaryOps.TIMES
    )
    op = getoperator(op, optype(eltype(A), eltype(B)))
    t = juliatype(ztype(op))
    C = GBMatrix{t}(size(A,1) * size(B, 1), size(A, 2) * size(B, 2))
end


text:
## extended methods for computing covariance and scatter matrix

# auxiliary functions

function _symmetrize!(a::DenseMatrix)
    m, n = size(a)
    m == n || error("a must be a square matrix.")
    for j = 1:n
        @inbounds for i = j+1:n
            vl = a[i,j]
            vr = a[j,i]
            a[i,j] = a[j,i] = middle(vl, vr)
        end
    end
    return a
end

function _scalevars(x::DenseMatrix, s::DenseVector, vardim::Int)
    vardim == 1 ? Diagonal(s) * x :
    vardim == 2 ? x * Diagonal(s) :
    error("vardim should be either 1 or 2.")
end

## scatter matrix

scattermat_zm(x::DenseMatrix, vardim::Int) = Base.unscaled_covzm(x, vardim)

scattermat_zm(x::DenseMatrix, wv::WeightVec, vardim::Int) =
    _symmetrize!(Base.unscaled_covzm(x, _scalevars(x, values(wv), vardim), vardim))

"""
    scattermat(X, [wv::WeightVec]; mean=nothing, vardim=1)

Compute the scatter matrix, which is an unnormalized covariance matrix.
A weighting vector `wv` can be specified to weight
the estimate.

# Arguments
* `mean=nothing`: a known mean value. `nothing` indicates that the mean is
unknown, and the function will compute the mean. Specifying `mean=0` indicates
that the data are centered and hence there's no need to subtract the mean.
* `vardim=1`: the dimension along which the variables are organized.
When `vardim = 1`, the variables are considered columns with observations in rows;
when `vardim = 2`, variables are in rows with observations in columns.
"""
function scattermat end

"""
    cov(X, wv::WeightVec; mean=nothing, vardim=1)

Compute the weighted covariance matrix. By default, the covariance
matrix is normalized by the sum of the weights. That is, `cov(X, wv)`
is equivalent to `scattermat(X, wv) / sum(wv)`.
"""
cov

"""
    mean_and_cov(x, [wv::WeightVec]; vardim=1) -> (mean, cov)

Return the mean and covariance matrix as a tuple. A weighting
vector `wv` can be specified. `vardim` that designates whether
the variables are columns in the matrix (`1`) or rows (`2`).
"""
function mean_and_cov end

@static if VERSION < v"0.5.0-dev+679"
    function scattermat(x::DenseMatrix; mean=nothing, vardim::Int=1)
        mean == 0 ? scattermat_zm(x, vardim) :
        mean == nothing ? scattermat_zm(x .- Base.mean(x, vardim), vardim) :
        scattermat_zm(x .- mean, vardim)
    end

    function scattermat(x::DenseMatrix, wv::WeightVec; mean=nothing, vardim::Int=1)
        mean == 0 ? scattermat_zm(x, wv, vardim) :
        mean == nothing ? scattermat_zm(x .- Base.mean(x, wv, vardim), wv, vardim) :
        scattermat_zm(x .- mean, wv, vardim)
    end

    ## weighted cov
    Base.cov(x::DenseMatrix, wv::WeightVec; mean=nothing, vardim::Int=1) =
        scale!(scattermat(x, wv; mean=mean, vardim=vardim), inv(sum(wv)))

    function mean_and_cov(x::DenseMatrix; vardim::Int=1)
        m = mean(x, vardim)
        return m, Base.covm(x, m; vardim=vardim)
    end
    function mean_and_cov(x::DenseMatrix, wv::WeightVec; vardim::Int=1)
        m = mean(x, wv, vardim)
        return m, Base.cov(x, wv; mean=m, vardim=vardim)
    end
else
    scattermatm(x::DenseMatrix, mean, vardim::Int=1) =
        scattermat_zm(x .- mean, vardim)

    scattermatm(x::DenseMatrix, mean, wv::WeightVec, vardim::Int=1) =
        scattermat_zm(x .- mean, wv, vardim)

    scattermat(x::DenseMatrix, vardim::Int=1) =
        scattermatm(x, Base.mean(x, vardim), vardim)

    scattermat(x::DenseMatrix, wv::WeightVec, vardim::Int=1) =
        scattermatm(x, Base.mean(x, wv, vardim), wv, vardim)

    ## weighted cov
    Base.covm(x::DenseMatrix, mean, wv::WeightVec, vardim::Int=1) =
        scale!(scattermatm(x, mean, wv, vardim), inv(sum(wv)))

    Base.cov(x::DenseMatrix, wv::WeightVec, vardim::Int=1) =
        Base.covm(x, Base.mean(x, wv, vardim), wv, vardim)

    function mean_and_cov(x::DenseMatrix, vardim::Int=1)
        m = mean(x, vardim)
        return m, Base.covm(x, m, vardim)
    end
    function mean_and_cov(x::DenseMatrix, wv::WeightVec, vardim::Int=1)
        m = mean(x, wv, vardim)
        return m, Base.cov(x, wv, vardim)
    end
end