text:
module KMLPolynetTools

using Meshes

#types
export Polynet, Region, triangulate

#methods
export extract_polynet_from_kml, load, save, scaled_svg
export get_or_cache_polynet, inRegions

# dependencies
using Serialization
using LightXML
using SVG

struct Region{T}
    meta::Dict
    areas::Vector{T}
    Region{PolyArea}(m::Dict) = new{PolyArea}(m, Vector{PolyArea}())
    Region{PolyArea}(m::Dict, pa::Vector{PolyArea}) = new{PolyArea}(m, pa)
    Region{SimpleMesh}(m::Dict) = new{SimpleMesh}(m, Vector{SimpleMesh}())
    Region{SimpleMesh}(m::Dict, sm::Vector{SimpleMesh}) = new{SimpleMesh}(m, sm)
end

const Polynet{T} = Vector{Region{T}}

Base.copy(r::Region) = Region(copy(r.meta), copy(r.areas))

function split_at_intersections(points)
    match(n) = t->t==n

    polys = Vector{typeof(points)}()
    poly = Vector{eltype(points)}()
    s = 1
    i = 1
    while i < length(points)
        push!(poly, points[i])
        n = findlast(match(points[i]), points[i+1:end-1])
        if n !== nothing
            append!(polys, split_at_intersections(points[i:i+n]))
            i = i + n
        end
        i += 1
    end
    push!(poly, points[end])
    push!(polys, poly)
    polys
end

function add_perimeters!(r::Region, points)
    for poly in split_at_intersections(remove_repeats(points))
        if length(poly) > 3
            push!(r.areas, PolyArea(poly))
        end
    end
    r
end

function remove_repeats(src)
    # src = [1,2,3,4,4,4,5,6,6,7,8,8,1]
    # tgt = [1,2,3,4,5,6,7,8,1]

    tgt = Vector{eltype(src)}(undef, length(src))
    srci = 1
    tgti = 0
    lastn = 0
    while srci < length(src)
        if src[srci] != lastn
            tgti += 1
            lastn = tgt[tgti] = src[srci]
        end
        srci += 1
    end
    tgti += 1
    tgt[tgti] = src[end]
    tgt[1:tgti]
end

txtPoint2(txt; digits=5) = Point2(map(t->round(parse(Float64, t); digits), split(txt, ",")))

function load(fn)::Union{Polynet, Nothing}
    pm = nothing
    if filesize(fn) > 0
        open(fn, "r") do io
            pm = deserialize(fn)
        end
    end
    pm
end

function save(fn, pm::Polynet)::Polynet
    open(fn, "w+") do io
        serialize(io, pm)
    end
    pm
end

import Meshes.boundingbox

function boundingbox(pnet::Polynet)
    bxs = convert(Vector{Box}, filter(bx->bx !== nothing, map(boundingbox, pnet)))
    length(bxs) > 0 ? boundingbox(bxs) : nothing
end
boundingbox(sms::Vector{SimpleMesh}) =  boundingbox(map(boundingbox, sms))
boundingbox(region::Region) = length(region.areas) > 0 ? boundingbox(map(boundingbox, region.areas)) : nothing

import Base.in

function in(p::Point2, r::Region)
    for a in r.areas
        if p in a
            return true
        end
    end
    false
end

inRegions(p::Point, pnet::Polynet) = pnet[map(r->p in r, pnet)]

function scaled_svg(pnet, width, height, filename; inhtml=true, digits=3, colorfn=nothing)
    bbx = boundingbox(pnet)
    xmin, ymin = coordinates(bbx.min)
    xmax, ymax = coordinates(bbx.max)
    xmx = xmax - xmin
    ymx = ymax - ymin
    scale = min(width, height) / min(xmx, ymx)
    xmx *= scale
    ymx *= scale
    fx = x -> round(scale * (x - xmin); digits)
    fy = y -> round(ymx - scale * (y - ymin); digits)

    asSvg(pnet, width, height, filename; fx, fy, xmx, ymx, inhtml, colorfn)
end

function asSvg(pnet::Polynet{PolyArea}, width, height, filename; fx=identity, fy=identity, xmx=0, ymx=0, colorfn=nothing, inhtml=true)
    if xmx == 0
        xmx = width
    end
    if ymx == 0
        ymx = height
    end
    if colorfn === nothing
        colorfn = (m)->"none"
    end

    pline(meta, polyarea) = Polyline(coordinates.(polyarea.outer.vertices), fx, fy; style=Style(;fill=colorfn(meta)))
    w = (io, svg) -> foreach(reg->foreach(polya->write(io, pline(reg.meta, polya)), reg.areas), pnet)
    SVG.write(filename, SVG.Svg(), width, height ; viewbox="0 0 $xmx $ymx", inhtml, objwrite_fn=w)
end

function asSvg(mnet::Polynet{SimpleMesh}, width, height, filename; fx=identity, fy=identity, xmx=0, ymx=0, colorfn=nothing, inhtml=true)
    if xmx == 0
        xmx = width
    end
    if ymx == 0
        ymx = height
    end
    if colorfn === nothing
        colorfn = (m)->"none"
    end

    function pline(meta, tri)
        coords = coordinates.(vertices(tri))
        Polyline(vcat(coords, [coords[1]]), fx, fy; style=Style(;fill=colorfn(meta)))
    end
    w = (io, svg) -> foreach(reg->foreach(smesh->foreach(tri->write(io, pline(reg.meta, tri)), smesh), reg.areas), mnet)
    SVG.write(filename, SVG.Svg(), width, height ; viewbox="0 0 $xmx $ymx", inhtml, objwrite_fn=w)
end

function polynet_from_kml(xdoc; digits=5)
    pnet = Polynet{PolyArea}()
    for fr in get_elements_by_tagname(get_elements_by_tagname(root(xdoc), "Document")[1], "Folder")
        for pk in get_elements_by_tagname(fr, "Placemark")
            meta = Dict{String, Union{String, Float64}}()
            for ed in get_elements_by_tagname(pk, "ExtendedData")
                for scd in get_elements_by_tagname(ed, "SchemaData")
                    for sd in get_elements_by_tagname(scd, "SimpleData")
                        for a in attributes(sd)
                            if LightXML.name(a) == "name"
                                meta[value(a)] = content(sd)
                            end
                        end
                    end
                end
            end
            region = Region{PolyArea}(meta)
            for mg in get_elements_by_tagname(pk, "MultiGeometry")
                for pol in get_elements_by_tagname(mg, "Polygon")
                    for bound in get_elements_by_tagname(pol, "outerBoundaryIs")
                        for lr in get_elements_by_tagname(bound, "LinearRing")
                            for cords in get_elements_by_tagname(lr, "coordinates")
                                add_perimeters!(region, map(t->txtPoint2(t; digits), split(content(cords), " ")))
                            end
                        end
                    end
                end
            end
            push!(pnet, region)
        end
    end
    pnet
end

function get_or_cache_polynet(kml, cachefn; digits=5, force=false)
    pnet::Union{Polynet, Nothing} = nothing
    if !force
        pnet = load(cachefn)
    end
    if pnet === nothing
        pnet = polynet_from_kml(parse_file(kml); digits)
        save(cachefn, pnet)
    end
    pnet
end

function triangulate(pnet::Polynet)
    meshnet = Polynet{SimpleMesh}()
    for reg in pnet
        meshes = Vector{SimpleMesh}()
        for pa in reg.areas
            try
                push!(meshes, discretize(pa, Dehn1899()))
            catch
                nothing
            end
        end
        push!(meshnet, Region{SimpleMesh}(reg.meta, meshes))
    end
    meshnet
end

#==
using KMLPolynetTools
kml = "Local_Planning_Authorities_May_2021_UK_BFC.kml"
geodir = "/home/matt/wren/UkGeoData"
psj = "polynet_2dp.sj"
pnet = get_or_cache_polynet(joinpath(geodir, kml), joinpath(geodir, psj));

==#

###
end


text:
export PoissonSparseCGILU,solve!
using IterativeSolvers
using IncompleteLU

struct PoissonSparseCGILU
    spLxy::SparseMatrixCSC{Float64,Int64}
    iluLxy::Any

    function PoissonSparseCGILU(sp)
        tilu=@elapsed lu = ilu(sp, τ = 0.2)
        println("ILU factorization time:",tilu)
        new(sp,lu)
    end
end

function solve!(Xxy,Bxy,sp::PoissonSparseCGILU)
    p = sp.iluLxy
    b=view(Bxy,1:length(Bxy))
    x=view(Xxy,1:length(Xxy))
    a,cgilulog=IterativeSolvers.cg!(x,sp.spLxy,b,Pl = p;log=true,reltol = 1e-16,maxiter=1000)
    @show cgilulog
    cgilulog
end

ext: jl


text:
"""

```
eulermat(beta, gamma, theta, deg::String)
```

This routine was added because it defines
the rotation in the "human" way, an is thus used
to set the position of the fixed molecules. `deg` can only be `"degree"`, in which
case the angles with be considered in degrees. If no `deg` argument
is provided, radians are used.

That means: `beta` is a counterclockwise rotation around `x` axis.
            `gamma` is a counterclockwise rotation around `y` axis.
            `theta` is a counterclockwise rotation around `z` axis.

"""
function eulermat(beta, gamma, theta, deg::String)

   if deg != "degree"
     error("ERROR: to use radians just omit the last parameter")
   end

   beta = beta*π/180
   gamma = gamma*π/180
   theta = theta*π/180

   return eulermat( beta, gamma, theta )
end

function eulermat(beta::T, gamma::T, theta::T) where T
  c1 = cos(beta)
  s1 = sin(beta)
  c2 = cos(gamma)
  s2 = sin(gamma)
  c3 = cos(theta)
  s3 = sin(theta)
  @SMatrix [   c2*c3                -c2*s3                s2
              (c1*s3 + c3*s1*s2)    (c1*c3 - s1*s2*s3)   -c2*s1
              (s1*s3 - c1*c3*s2)    (c1*s2*s3 + c3*s1)    c1*c2  ]
end


text:
using ClusterTrees
using StaticArrays
using CompScienceMeshes

const P = SVector{3,Float64}

tree = ClusterTrees.Octrees.Octree(Int[])
mesh = meshsphere(1.0, 0.15)
# points = [rand(P) for i in 1:800]
points = vertices(mesh)

struct OTRouter
    smallest_box_size::Float64
    target_point::P
end

ClusterTrees.route!(t::ClusterTrees.Octrees.Octree, state, router) = ClusterTrees.Octrees.route!(t, state, router)

smallest_box_size = 0.1
root_sector = 0
root_center = P(0,0,0)
root_size = 1.0
for i in 1:length(points)
    router = OTRouter(smallest_box_size, points[i])
    root_state = root(tree), root_center, root_size, 1
    ClusterTrees.update!(tree, root_state, i, router) do tree, node, data
        push!(ClusterTrees.data(tree,node).values, data)
    end
end

ClusterTrees.print_tree(tree)

num_bins = 8
bins = [P[] for i in 1:num_bins]

num_printed = 0
num_points = 0
num_nodes = length(tree.pbtree.nodes)
for (i,node) in enumerate(ClusterTrees.DepthFirstIterator(tree, root(tree)))
    println(node, ", ", ClusterTrees.data(tree,node))
    b = div((i-1)*num_bins, num_nodes) + 1
    append!(bins[b], points[ClusterTrees.data(tree,node).values])
    global num_printed += 1
    global num_points += length(data(tree,node).values)
end

ordered_points = reduce(append!, bins)
ordered_points = [p[i] for p in ordered_points, i = 1:3]
@assert num_printed == length(tree.pbtree.nodes)
@show num_points
@show length(points)
@assert num_points == length(points)

# using Plots
# plot()
# for i in 1:length(bins)
#     scatter!(bins[i])
# end
# plot!(legend=false)


text:
# Copyright (c) 2021, Joseph M. Rutherford
# All rights reserved.
#
# Code provided under the license contained in the LICENSE file.

include("../src/sudoku.jl") # Defines Sudoku module

# Form a random solution and puzzle pair
using Random
solution, puzzle = Sudoku.random_puzzle(3,MersenneTwister(12345),10000,45)

# Solve puzzle and report findings
print("\nUnsolved puzzle\n")
print(Sudoku.as_text_grid(puzzle))
iterations, uncertainty = Sudoku.iterative_solve!(puzzle,2)
print("\n\nSolution after $iterations step(s)\n")
print(Sudoku.as_text_grid(puzzle))

if uncertainty != 0
    print("Failure: puzzle solution is incomplete\n")
else
    if Sudoku.as_values(puzzle) == Sudoku.as_values(solution)
        print("Success: solved puzzle matches expected solution\n")
    else
        print("Failure: solved puzzle does not match expected solution\n")
        print("Expected\n")
        print(Sudoku.as_text_grid(solution))
    end
end

max_stars_repo_name: joseph-m-rutherford/sudoku-julia

max_issues_repo_name: joseph-m-rutherford/sudoku-julia

max_forks_repo_name: joseph-m-rutherford/sudoku-julia


text:
immutable PolyUnion
    vars::Vector{Symbol}
    mat::Array{Int, 2}
    eltype::Type
end

function PolyUnion(ps)
    newvars = isempty(ps) ? [] : union([vars(p) for p in ps]...)
    mat = Array(Int, length(newvars), length(ps))
    for pindex = 1:length(ps)
        for i = 1:length(newvars)
            mat[i, pindex] = findfirst(vars(ps[pindex]), newvars[i])
        end
    end
    PolyUnion(newvars, mat, promote_type([eltype(p) for p in ps]...))
end

eltype(pu::PolyUnion) = pu.eltype

vars(pu::PolyUnion) = pu.vars

nvars(pu::PolyUnion) = length(vars(pu))

newexps(pu::PolyUnion, exps::Vector{Int}, polindex::Int) =
    [pu.mat[i, polindex] == 0 ? 0 : exps[pu.mat[i, polindex]] for i = 1:nvars(pu)]

function oldexps(pu::PolyUnion, exps::Vector{Int}, polindex::Integer)
    v = zeros(Int, length(exps))
    l = 0
    for i = 1:length(exps)
        t = pu.mat[i, polindex]
        if t != 0
            v[t] = exps[i]
            l += 1
        end
    end
    v[1:l]
end

zero(pu::PolyUnion) = zero(MPoly{eltype(pu)}, vars=vars(pu))

one(pu::PolyUnion) = one(MPoly{eltype(pu)}, vars=vars(pu))


text:
using Cassette

Cassette.@context CounterCtx;

const binops = (
    (:add32, Core.Intrinsics.add_float, Float32),
    (:sub32, Core.Intrinsics.sub_float, Float32),
    (:mul32, Core.Intrinsics.mul_float, Float32),
    (:div32, Core.Intrinsics.div_float, Float32),
    (:add64, Core.Intrinsics.add_float, Float64),
    (:sub64, Core.Intrinsics.sub_float, Float64),
    (:mul64, Core.Intrinsics.mul_float, Float64),
    (:div64, Core.Intrinsics.div_float, Float64),
)

const unops = (
    (:sqrt32, Core.Intrinsics.sqrt_llvm, Float32),
    (:sqrt64, Core.Intrinsics.sqrt_llvm, Float64),
)

const ops = Iterators.flatten((binops, unops)) |> collect

@eval mutable struct Counter
    $((:($(op[1]) ::Int) for op in ops)...)
    Counter() = new($((0 for _ in 1:length(ops))...))
end

for typ1 in (Float32, Float64)
    @eval function Cassette.prehook(ctx::CounterCtx,
                                    op::Core.IntrinsicFunction,
                                    ::$typ1,
                                    ::$typ1)
        $(Expr(:block,
               (map(binops) do (name, op, typ2)
                  typ1 == typ2 || return :nothing
                  quote
                    if op == $op
                       ctx.metadata.$name += 1
                       return
                    end
                  end
                end)...))
    end

    @eval function Cassette.prehook(ctx::CounterCtx,
                                    op::Core.IntrinsicFunction,
                                    ::$typ1)
        $(Expr(:block,
               (map(unops) do (name, op, typ2)
                  typ1 == typ2 || return :nothing
                  quote
                    if op == $op
                       ctx.metadata.$name += 1
                       return
                    end
                  end
                end)...))
    end
end

# Relatively inefficient, but there should be no need for performance here...

flop(c::Counter) = sum(getfield(c, field) for field in fieldnames(Counter))

import Base: ==, *, show

function Base.show(io::IO, c::Counter)
    println(io, "Flop Counter:")
    for field in fieldnames(Counter)
        println(io, " $field: $(getfield(c, field))")
    end
end

function ==(c1::Counter, c2::Counter)
    all(getfield(c1, field)==getfield(c2, field) for field in fieldnames(Counter))
end

function *(n::Int, c::Counter)
    ret = Counter()
    for field in fieldnames(Counter)
        setfield!(ret, field, n*getfield(c, field))
    end
    ret
end


text:
using DataFrames, DataFramesMeta, CSV
using FileIO, FCSFiles
using StatsBase

"""
    clip(X;p=[0.001,0.999])

Clips the data outside of range specified by `quantile(X,p)`.

    clip(X,Y;p=[0.001,0.999])

Clips the data outside of ranges specified by `quantile(X,p)` and `quantile(Y,p)`.

"""
function clip(X;p=[0.001,0.999])
    q = quantile(X,p)
    keep = (X .> q[1]) .& (X .< q[2])
    X[keep]
end
function clip(X,Y;p=[0.001,0.999])
    q1 = quantile(X,p)
    q2 = quantile(Y,p)
    keep = (X .> q1[1]) .& (X .< q1[2]) .& (Y .> q2[1]) .& (Y .< q2[2])
    X[keep],Y[keep]
end

# Create CSV containing all fcs data
files = readdir("Data/FCS")[2:end]
df1   = Array{Any,1}(undef,length(files))   # Store all data
df2   = Array{Any,1}(undef,length(files))   # Store summary (geometric means)
for (i,file) = enumerate(files)

    fcs  = load("Data/FCS/$file")
    df   = DataFrame(fcs.data)
    gain = Dict([fcs.params["\$P$(i)N"] => parse(Float64,fcs.params["\$P$(i)G"]) for i = 1:18]...)
    info = split(file,"_")  # Metadata from filename

    # Normalise BDP and compensatedd Cy5 signal
    signal_BDP = df[:,"488nm (520_40) LinH"] ./ gain["488nm (520_40) LinH"]
    signal_Cy5 = df[:,"FJComp-642nm (676_29) LinH"] ./ gain["642nm (676_29) LinH"]

    # Filter cells with no signal
    keep = (signal_BDP .> 0.0) .& (signal_Cy5 .> 0.0);
    signal_BDP = signal_BDP[keep]
    signal_Cy5 = signal_Cy5[keep]

    # Clip top and bottom 0.1%
    keep = (quantile(signal_BDP,0.001) .< signal_BDP .< quantile(signal_BDP,0.999)) .&
           (quantile(signal_Cy5,0.001) .< signal_Cy5 .< quantile(signal_Cy5,0.999))
    signal_BDP = signal_BDP[keep]
    signal_Cy5 = signal_Cy5[keep]

    # Convert to data frame (get metadata from filename)
    df1[i] = DataFrame(
        Exp      = fill("DualMarker",                       length(signal_BDP)),
        FCSCond  = fill(info[1],                            length(signal_BDP)),
        Temp     = fill(parse(Float64,info[2][1:end-6]),    length(signal_BDP)),
        CellLine = fill(info[4],                            length(signal_BDP)),
        Quenched = fill(info[3] == "Q",                     length(signal_BDP)),
        Time     = fill(parse(Float64,info[5][1:end-3]),    length(signal_BDP)),
        Repeat   = fill(parse(Int64,info[6][1:end-4]),      length(signal_BDP)),
        Signal_BDP   = signal_BDP,
        Signal_Cy5   = signal_Cy5
    )

    # Summarise data (using MFI/geomean)
    df2[i] = DataFrame(
        Exp      = "DualMarker",
        FCSCond  = info[1],
        Temp     = parse(Float64,info[2][1:end-6]),
        CellLine = info[4],
        Quenched = info[3] == "Q",
        Time     = parse(Float64,info[5][1:end-3]),
        Repeat   = parse(Int64,info[6][1:end-4]),
        CellCountNZ  = length(signal_Cy5),
        CellCount    = nrow(df),
        Signal_BDP   = geomean(signal_BDP),
        Signal_Cy5   = geomean(signal_Cy5),
    )

end

mkdir("Data/CSV")
CSV.write("Data/CSV/DualMarker.csv",vcat(df1...))
CSV.write("Data/CSV/DualMarker_Summary.csv",vcat(df2...))

max_stars_repo_name: ap-browning/internalisation

max_issues_repo_name: ap-browning/internalisation

max_forks_repo_name: ap-browning/internalisation


text:
module SSLTestWorldclim
using SimpleSDMLayers
using Test

wc1and2 = SimpleSDMPredictor(WorldClim, BioClim, [1,2])
@test typeof(first(wc1and2)) <: SimpleSDMPredictor

temp = SimpleSDMPredictor(WorldClim, BioClim, 1)
@test size(temp) == (1080, 2160)
@test round(first(stride(temp)); digits=2) ≈ round(last(stride(temp)); digits=2)
@test length(longitudes(temp)) == 2160
@test length(latitudes(temp)) == 1080

wc3 = SimpleSDMPredictor(WorldClim, BioClim, 3)
@test typeof(wc3) <: SimpleSDMPredictor

wcrange = SimpleSDMPredictor(WorldClim, BioClim, 1:5)
@test eltype(wcrange) <: SimpleSDMPredictor

future = SimpleSDMPredictor(WorldClim, BioClim, MIROC6, SSP126, 1)
@test typeof(future) <: SimpleSDMPredictor

future = SimpleSDMPredictor(WorldClim, BioClim, CanESM5, SSP126, 1:3)
@test eltype(future) <: SimpleSDMPredictor

future = SimpleSDMPredictor(WorldClim, BioClim, CanESM5, SSP370, [1,2,5]; year="2061-2080")
@test eltype(future) <: SimpleSDMPredictor

end


text:
struct InteractivePlay end

hear(s::Scene{InteractivePlay}, ::Genesis!) = for a in stage_ref(s).state.actors
    a.state isa Person && put!(a.inbox, Genesis!())
    break
end

struct Person end

"Start a minimal actor system which you can interact with from the REPL"
interact!() = interact!(InteractivePlay())

"""Like [`play!`](@ref), but allows you to send and receive messages from the REPL

Calling this will asynchronously start a `play` and return a 2-tuple of
[`Scene`](@ref) and `Task`. The `Scene` can be used to call various methods in
the actor system and the task can be checked to see if the actor system has
failed.

The function [`local_addresses`](@ref) can be used to find out what
[`Id`](@ref) a given actor has.

!!! warning

    This is only intended for interactive use in the Julia REPL or similar
    scenario.

### Example

```jldoctest
julia> using Actors

julia> struct Ponger end

julia> struct Echo
       who::Id
       msg::String
       end

julia> Actors.hear(s::Scene{Ponger}, msg::Echo) = say(s, msg.who, msg.msg)

julia> (s, t) = interact!()
(Scene{Actors.Person}, Task (runnable) @0x00007fafc7dad600)

julia> me(s)
Actors.Person@1

julia> ponger = enter!(s, Ponger())
Ponger@5

julia> ask(s, ponger, Echo(me(s), "Echo, echo, ..."), String)
"Echo, echo, ..."

julia> local_addresses(s)
5-element Array{Pair{Int64,DataType},1}:
 1 => Actors.Person
 2 => Actors.Logger{Base.TTY}
 3 => Actors.PassiveMinder{Id{Actors.Logger{Base.TTY},Union{Leave!, Actors.LogDied!, LogInfo!}}}
 4 => Actors.InteractivePlay
 5 => Ponger

```
"""
function interact!(play)
    inb = Channel{Any}(512)
    env = capture_environment(id)
    st_id = Id{Stage}(UInt32(0))
    local a_id

    task = @async begin
        st = Actor(Channel{Any}(1024), Stage(), st_id, current_task(), st_id)
        a_id = push!(st.state.actors, Actor(inb, Person(), st_id))

        put!(st.inbox, PreGenesis!(play))
        put!(inb, st)

        play!(Scene(st, st), env)
    end

    st = take!(inb)
    a = Actor(inb, Person(), st_id, current_task(), a_id)
    st.state.actors[a_id] = a

    yield()

    Scene(a, st), task
end

listen!(s::Scene{Person}) = take!(inbox(s))

"""Return a list of all the [`Id`](@ref)s and what type of [`Actor`](@ref) they point to

!!! warning

    This is only intended for interactive use in the REPL or for diagnostics.
"""
local_addresses(s::Scene) = [i => typeof(a.state) for (i, a) in enumerate(stage_ref(s).state.actors)]


text:
# Read the spectrum data file format associated with handheld Bruker XRF units.

function readbrukerpdz(fn::AbstractString)
    open(fn, "r") do f
        res = readbrukerpdz(f)
        res[:Name] = splitpath(fn)[end]
        return res
    end
end

"""
    readbrukerpdz(io::IO|fn::AbstractString)

Read the default data file type for the Bruker handheld XRF unit. (Reverse engineered.)
"""
function readbrukerpdz(io::IO)
    mm = Mmap.mmap(io)
    nCh = reinterpret(UInt16, mm[7:8])[1]
    eVpCh = reinterpret(Float64, mm[51:58])[1]
    yr, mon, _, day, hour, min, sec = reinterpret(UInt16, mm[147:160])
    e0, pc = reinterpret(Float32, mm[163:170])
    rt, dt, _, lt = reinterpret(Float32, mm[343:358])
    filt = reinterpret(UInt16, mm[115:130]) # Filter / thickness
    filters = [
        (PeriodicTable.elements[filt[2*i-1]], filt[2*i] * 1.0e-6) for
        i in filter(i -> filt[2*i-1] > 0, 1:4)
    ]
    dt = DateTime(yr, mon, day, hour, min, sec)
    props = Dict{Symbol,Any}(
        :RealTime => convert(Float64, rt),
        :LiveTime => convert(Float64, lt),
        :AcquisitionTime => dt,
        :BeamEnergy => convert(Float64, e0) * 1000.0,
        :ProbeCurrent => convert(Float64, pc) * 1000.0,
        :XRFFilter => filters,
    )
    counts = reinterpret(UInt32, mm[length(mm)-4*nCh+1:length(mm)])
    return Spectrum(LinearEnergyScale(0.0, eVpCh), collect(counts), props)
end

function detectBrukerPDZ(ios::IO)
    try
        return read(io, UInt32) == 0x00170101
    catch
        return false
    end
end

max_stars_repo_name: NicholasWMRitchie/NeXLSpectrum

max_issues_repo_name: NicholasWMRitchie/NeXLSpectrum

max_forks_repo_name: NicholasWMRitchie/NeXLSpectrum


text:
using EffectiveWaves, Test

# large volume fraction scatterers,  small size amd on strong scatterer. This is a problamatic case.
@testset "large volume fraction and low frequency" begin
    # should pass for the list of angular frequencies below
    # ωs = [0.001, 2.0, 9.0]
    ω = 0.001

    medium = Medium(1.0,1.0+0.0im)
    species = [
        Specie(ρ=5.,r=0.004, c=1.2, volfrac=0.4),
        Specie(ρ=0.3, r=0.002, c=0.4, volfrac=0.3)
    ]

    eff_medium = effective_medium(medium, species)

    tol = 1e-7
    k_eff_low = ω/eff_medium.c
    k_effs = wavenumbers(ω, medium, species; tol = tol, num_wavenumbers=1)
    i = findmin(abs.(k_effs .- k_eff_low))[2]
    k_eff = k_effs[i]

    @test abs(k_eff - k_eff_low)/norm(k_eff_low) < 10*tol

    R = begin
        wave = EffectiveWave(ω, k_eff, medium, species)
        reflection_coefficient(ω, wave, medium, species)
    end
    R_low2 = begin
        wave = EffectiveWave(ω, k_eff_low, medium, species)
        reflection_coefficient(ω, wave, medium, species)
    end
    R_low = reflection_coefficient_halfspace(medium, eff_medium)

    @test norm(R_low - R) < tol
    @test norm(R_low2 - R) < tol
end


text:
using Regression
using SVM
using Optim

require("preprocess.jl")

require("svm.jl")

X = randn(100, 5000)

I = randbool(100, 5000)
X[I] = 0.0

X[1,:] = 1.0

theta_star = randn(100)

y = 1.0*((X' * theta_star + 0.05 * randn(5000)) .>= 0.0)

data = [X; y']

num_features = size(X, 1)

theta = randn(num_features)/100

lambda = 0.001

function func(x::Vector)
  return loss(data, x, lambda)
end

function gradient!(x::Vector, storage::Vector)
  copy!(storage, full_grad(grad!, data, x, lambda))
end

res = optimize(func, gradient!, theta, method = :l_bfgs; grtol=1e-13)

dataset_info = ["name" => "synthetic", "theta_star" => res.minimum]using Regression

max_stars_repo_name: pcmoritz/slbfgs

max_issues_repo_name: pcmoritz/slbfgs

max_forks_repo_name: pcmoritz/slbfgs


text:
#=
"""
```
hessizero{T<:AbstractFloat}(fcn::Function, x::Vector{T};
                            check_neg_diag::Bool=false,
                            verbose::Symbol=:none,
                            distr::Bool=true)
```

Compute Hessian of function `fcn` evaluated at `x`.

### Arguments
- `check_neg_diag`: Throw an error if any negative diagonal elements are detected.
- `verbose`: Print verbose output
- `distr`: Use available parallel workers to increase performance.
"""
function hessizero{T<:AbstractFloat}(fcn::Function,
                                    x::Vector{T};
                                    check_neg_diag::Bool=false,
                                    verbose::Symbol=:none,
                                    distr::Bool=true)
    n_para = length(x)
    hessian  = zeros(n_para, n_para)

    # Compute diagonal elements first
    if distr && nworkers() > 1
        diag_elements = @sync @distributed (hcat) for i = 1:n_para
            hess_diag_element(fcn, x, i; check_neg_diag=check_neg_diag, verbose=verbose)
        end
        for i = 1:n_para
            hessian[i, i] = diag_elements[i]
        end
    else
        for i=1:n_para
            hessian[i,i] = hess_diag_element(fcn, x, i; check_neg_diag=check_neg_diag,
                                             verbose=verbose)
        end
    end

    # Now compute off-diagonal elements
    # Make sure that correlations are between -1 and 1
    # invalid_corr indexes elements that are invalid
    invalid_corr = Dict{Tuple{Int,Int}, Float64}()

    # Build indices to iterate over
    n_off_diag_els = Int(n_para*(n_para-1)/2)

    #this allows for base case where there is only 1 free parameter
    if n_off_diag_els > 0

        off_diag_inds = Vector{Tuple{Int,Int}}(n_off_diag_els)
        k=1
        for i=1:(n_para-1), j=(i+1):n_para
            off_diag_inds[k] = (i,j)
            k = k+1
        end

        # Iterate over off diag elements
        if distr
            off_diag_out = @sync @distributed (hcat) for (i,j) in off_diag_inds
                σ_xσ_y = sqrt(abs(hessian[i, i]*hessian[j, j]))
                hess_offdiag_element(fcn, x, i, j, σ_xσ_y; verbose=verbose)
            end
            # Ensure off_diag_out is array
            off_diag_out = hcat(off_diag_out)
        else
            off_diag_out = Array{Tuple{T, T},1}(n_off_diag_els)
            for (k,(i,j)) in enumerate(off_diag_inds)
                σ_xσ_y = sqrt(abs(hessian[i, i]*hessian[j, j]))
                off_diag_out[k] = hess_offdiag_element(fcn, x, i, j, σ_xσ_y; verbose=verbose)
            end
        end

        # Fill in values
        for k=1:n_off_diag_els
            (i,j) = off_diag_inds[k]
            (value, ρ_xy) = off_diag_out[k]

            hessian[i,j] = value
            hessian[j,i] = value

            if ρ_xy < -1 || 1 < ρ_xy
                invalid_corr[(i, j)] = ρ_xy
            end
        end
    end

    has_errors = false
    if !isempty(invalid_corr)
        println("Errors: $invalid_corr")
        has_errors = true
    end

    return hessian, has_errors
end
=#


text:
using Test

if !@isdefined RandomNumbers
    include("../common.jl")
end

@testset "MersenneTwisters" begin

    using RandomNumbers.MersenneTwisters

    @info "Testing MersenneTwisters"
    stdout_ = stdout
    pwd_ = pwd()
    cd(dirname(@__FILE__))
    rm("./actual"; force=true, recursive=true)
    mkpath("./actual")

    @test seed_type(MT19937) == NTuple{RandomNumbers.MersenneTwisters.N, UInt32}

    for mt_name in (:MT19937, )

        outfile = open(string(
            "./actual/check-$(lowercase("$mt_name")).out"
        ), "w")
        redirect_stdout(outfile)

        @eval $mt_name()
        x = @eval $mt_name(123)
        @test copyto!(copy(x), x) == x

        for i in 1:100
            @printf "%.9f\n" rand(x)
        end

        close(outfile)
    end
    redirect_stdout(stdout_)

    compare_dirs("expected", "actual")

    # Package content should not be modified.
    rm("./actual"; force=true, recursive=true)

    cd(pwd_)

end


text:
using WaterWaves
using Test

@testset "WaterWaves.jl" begin
    # Write your own tests here.
    T = 6.;
    ω = 2π/T;
    g=9.81;
    μ0 = ω^2/g;
    h=1000.;
    kappa = μ0;
    @test waveNumber(0,μ0,h)==kappa
    @test isapprox(waveNumber(0,μ0,h), waveNumberHunt(μ0,h), atol = 1e-12)
    #println(waveNumber(0,μ0,h))
    #println(waveNumberHunt(μ0,h))
    h=0.5
    kappa = sqrt(ω^2/(g*h))
    @test isapprox(waveNumber(0,μ0,h), waveNumberHunt(μ0,h), atol = 1e-11)
    #println(waveNumber(0,μ0,h))
    #println(waveNumberHunt(μ0,h))
end;


text:
using CSV
using DataFrames
using DelimitedFiles
using Glob
using Plots; pgfplotsx()

function create_subplot(nthreads, kernels, systems)

	# Create plot
	plt = plot(ylabel="GFlops", xaxis=:log10, xlabel="Matrix size",
			   legend=:topleft, framestyle=:box)

	# Linestyles
	linestyles = Dict("sys1, par, 32" => :solid, "sys1, mkl, 32" => :dash,
					  "sys2, par, 32" => :dash,   "sys1, par, 16" => :dash,
					  "sys1, par, 8" => :dashdot,   "sys1, par, 4" => :dashdotdot,
					  "sys1, par, 1" => :dot)

	# Colors
	colors = Dict("sys1, par, 32" => :blue, "sys1, mkl, 32" => :green,
				  "sys2, par, 32" => :darkorange3,  "sys1, par, 16" => :purple,
				  "sys1, par, 8" => :red,  "sys1, par, 4" => :fuchsia,
				  "sys1, par, 1" => :grey)

	# Labels
	syslabel = Dict("sys1" => "sys. I", "sys2" => "sys. II")
	kernellabel = Dict("par" => "MtSpMV.jl", "mkl" => "MKLSparse.jl")

	# Count
	counter = 0

	# Iterate over number of threads
	for (it, t) in enumerate(nthreads)

		# Iterate over kernel
		for (ikernel, kernel) in enumerate(kernels)

			# Iterate over systems
			for (isys, sys) in enumerate(systems)

				# Count up
				counter += 1

				# Add number of threads to label
				if t == "1"
				    label = "SparseArrays, 1 thread, "
				    label = string(label, syslabel[sys])
				else
				    label = ""
				    label = string(label, kernellabel[kernel])
				    label = string(label, ", ")
				    label = string(label, t, " threads")
				    label = string(label, ", ")
				    label = string(label, syslabel[sys])
			    end

				# Get folder with benchmarks
				folder = joinpath(sys, kernel, t)
				folder = joinpath(folder, readdir(folder)[1], "csv")

				# Get all csv files for metric in sorted order
				file = glob("gflops*", folder)[1]

				# Get benchmark results starting from 10^2
			    df = DataFrame!(CSV.File(file, header=true))#[9:65,:]
			    x = convert(Vector, df[!, 2])
			    y = convert(Vector, df[!, 5])

				# Add to plot
				key = string(sys, ", ", kernel, ", ", t)
				plot!(plt, x, y, label=label, linestyle=linestyles[key],
					  color=colors[key])

			end #for

		end # for

	end # for

	return plt

end # function

# Create plots
function plot_spmv_size()

	# Plot 1
	plt1 = create_subplot(["32"], ["par", "mkl"], ["sys1"])

	# Plot 2
	plt2 = create_subplot(["32"], ["par"], ["sys1", "sys2"])

	# Plot 3
	plt3 = create_subplot(["32", "16", "8", "4", "1"], ["par"], ["sys1"])
	# plt3 = plot(plt3, xlabel="Matrix size")

	# Plot side by side and save
	width = 520
	height = 700
	plt = plot(plt1, plt2, plt3, layout=(3,1), size=(width, height))
	filename = string("bench_spmv_size")
	savefig(string(filename, ".png"))
	savefig(string(filename, ".tex"))

end # function

max_stars_repo_name: jkrch/JuliaParGMRES-Plots

max_issues_repo_name: jkrch/JuliaParGMRES-Plots

max_forks_repo_name: jkrch/JuliaParGMRES-Plots


text:
#level 2 matrix vector multiplication

################################################################################
## GBMV - general band matrix vector multiply

function BLAS.gbmv!{T<:Posit}(tA::Char, m::Integer, kl::Integer, ku::Integer, alpha::T, A::Matrix{T}, x::Vector{T}, beta::T, y::Vector{T})
    if tA == 'N' || tA == 'n'
        ngbmv!(m, kl, ku, alpha, A, x, beta, y)
    elseif tA == 'T' || tA == 't' || tA == "C" || tA == "c"
        tgbmv!(m, kl, ku, alpha, A, x, beta, y)
    else
        throw(ArgumentError("unknown tA specification for gbmv"))
    end
end

function BLAS.gbmv!{N,ES}(tA::Char,
        m::Integer,
        kl::Integer,
        ku::Integer,
        alpha::PositOrComplex{N,ES},
        A::Matrix{Complex{Posit{N,ES}}},
        x::Vector{Complex{Posit{N,ES}}},
        beta::PositOrComplex{N,ES},
        y::Vector{Complex{Posit{N,ES}}})

    println("using special blas!")

    if tA == 'N' || tA == 'n'
        ngbmv!(m, kl, ku, alpha, A, x, beta, y)
    elseif tA == 'T' || tA == 't'
        tgbmv!(m, kl, ku, alpha, A, x, beta, y)
    elseif tA == "C" || tA == "c"
        cgbmv!(m, kl, ku, alpha, A, x, beta, y)
    else
        throw(ArgumentError("unknown tA specification for gbmv"))
    end
end

function ngbmv!{N, ES}(
    m::Integer,
    kl::Integer,
    ku::Integer,
    alpha::PositOrComplex{N,ES},
    A::Matrix{<:PositOrComplex{N,ES}},
    x::Vector{<:PositOrComplex{N,ES}},
    beta::PositOrComplex{N,ES},
    y::Vector{<:PositOrComplex{N,ES}})

    #dimension checking.
    if (size(A, 2) != length(x))
        throw(DimensionMismatch("matrix A has dimension 2 $(size(A,2)), but x has length $(length(x))"))
    elseif (size(A, 1) != length(y))
        throw(DimensionMismatch("matrix A has dimension 1 $(size(A,1)), but y has length $(length(y))"))
    end

    for row = 1:size(A,1)
        res = naked_dot(A[start:stop, row], x[start:stop], alpha, beta, y[row])
        y[row] = res
    end
end

function tgbmv!{N,ES}(
    m::Integer,
    kl::Integer,
    ku::Integer,
    alpha::PositOrComplex{N,ES},
    A::Matrix{<:PositOrComplex{N,ES}},
    x::Vector{<:PositOrComplex{N,ES}},
    beta::PositOrComplex{N,ES},
    y::Vector{<:PositOrComplex{N,ES}})

    #dimension checking.
    if (size(A, 1) != length(x))
        throw(DimensionMismatch("matrix A' has dimension 2 $(size(A,2)), but x has length $(length(x))"))
    elseif (size(A, 1) != m)
        throw(DimensionMismatch("matrix A' has dimension 1 $(size(A,1)), but passed value m is $(m)"))
    elseif (size(A, 2) != length(y))
        throw(DimensionMismatch("matrix A' has dimension 1 $(size(A,1)), but y has length $(length(y))"))
    end

    for row = 1:size(A,1)
        start = clamp(row - kl, 1, m)
        stop =  clamp(row + kl, 1, m)
        res = naked_dot(A[row, start:stop], x[start:stop], alpha, beta, y[row])
        y[row] = res
    end
end

function cgbmv!{N,ES}(
    m::Integer,
    kl::Integer,
    ku::Integer,
    alpha::PositOrComplex{N,ES},
    A::Matrix{<:PositOrComplex{N,ES}},
    x::Vector{<:PositOrComplex{N,ES}},
    beta::PositOrComplex{N,ES},
    y::Vector{<:PositOrComplex{N,ES}})

    #dimension checking.
    if (size(A, 1) != length(x))
        throw(DimensionMismatch("matrix A' has dimension 2 $(size(A,2)), but x has length $(length(x))"))
    elseif (size(A, 2) != length(y))
        throw(DimensionMismatch("matrix A' has dimension 1 $(size(A,1)), but y has length $(length(y))"))
    end

    for row = 1:size(A,1)
        start = clamp(row - kl, 1, m)
        stop =  clamp(row + kl, 1, m)
        res = naked_dot(A[row, start:stop], x[start:stop], alpha, beta, y[row], Val{:conj})
        y[row] = res
    end
end

################################################################################
## GEMV - general matrix vector multiply

function BLAS.gemv!{T<:Posit}(tA::Char, alpha::T, A::Matrix{T}, x::Vector{T}, beta::T, y::Vector{T})

    if tA == 'N' || tA == 'n'
        ngemv!(alpha, A, x, beta, y)
    elseif tA == 'T' || tA == 't' || tA == "C" || tA == "c"
        tgemv!(alpha, A, x, beta, y)
    else
        throw(ArgumentError("unknown tA specification for gemv"))
    end
end

function BLAS.gemv!{N,ES}(tA::Char,
        alpha::PositOrComplex{N,ES},
        A::Matrix{Complex{Posit{N,ES}}},
        x::Vector{Complex{Posit{N,ES}}},
        beta::PositOrComplex{N,ES},
        y::Vector{Complex{Posit{N,ES}}})

    if tA == 'N' || tA == 'n'
        ngemv!(alpha, A, x, beta, y)
    elseif tA == 'T' || tA == 't'
        tgemv!(alpha, A, x, beta, y)
    elseif tA == "C" || tA == "c"
        cgemv!(alpha, A, x, beta, y)
    else
        throw(ArgumentError("unknown tA specification for gemv"))
    end
end

################################################################################

function ngemv!{N,ES}(
    alpha::PositOrComplex{N,ES},
    A::Matrix{<:PositOrComplex{N,ES}},
    x::Vector{<:PositOrComplex{N,ES}},
    beta::PositOrComplex{N,ES},
    y::Vector{<:PositOrComplex{N,ES}})

    #dimension checking.
    if (size(A, 2) != length(x))
        throw(DimensionMismatch("matrix A has dimension 2 $(size(A,2)), but x has length $(length(x))"))
    elseif (size(A, 1) != length(y))
        throw(DimensionMismatch("matrix A has dimension 1 $(size(A,1)), but y has length $(length(y))"))
    end

    for row = 1:size(A,1)
        res = naked_dot(A[:, row], x, alpha, beta, y[row])
        y[row] = res
    end
end

function tgemv!{N,ES}(
    alpha::PositOrComplex{N,ES},
    A::Matrix{<:PositOrComplex{N,ES}},
    x::Vector{<:PositOrComplex{N,ES}},
    beta::PositOrComplex{N,ES},
    y::Vector{<:PositOrComplex{N,ES}})

    #dimension checking.
    if (size(A, 1) != length(x))
        throw(DimensionMismatch("matrix A' has dimension 2 $(size(A,2)), but x has length $(length(x))"))
    elseif (size(A, 2) != length(y))
        throw(DimensionMismatch("matrix A' has dimension 1 $(size(A,1)), but y has length $(length(y))"))
    end

    for row = 1:size(A,1)
        res = naked_dot(A[row, :], x, alpha, beta, y[row])
        y[row] = res
    end
end

function cgemv!{N,ES}(
    alpha::PositOrComplex{N,ES},
    A::Matrix{<:PositOrComplex{N,ES}},
    x::Vector{<:PositOrComplex{N,ES}},
    beta::PositOrComplex{N,ES},
    y::Vector{<:PositOrComplex{N,ES}})

    #dimension checking.
    if (size(A, 1) != length(x))
        throw(DimensionMismatch("matrix A' has dimension 2 $(size(A,2)), but x has length $(length(x))"))
    elseif (size(A, 1) != length(y))
        throw(DimensionMismatch("matrix A has dimension 1 $(size(A,1)), but y has length $(length(y))"))
    end

    for row = 1:size(A,1)
        res = naked_dot(A[row, :], x, alpha, beta, y[row], Val{:conj})
        y[row] = res
    end
end

#hermitian matrices
#hbmv
#hemv
#hpmv

#symmetric band matrices
#sbmv
#spmv

#symmetric matrices
#symv

#triangular matrices
#tbmv
#tpmv
#trmv

max_stars_repo_name: ityonemo/SigmoidNumbers

max_issues_repo_name: ityonemo/SigmoidNumbers

max_forks_repo_name: ityonemo/SigmoidNumbers


text:
using LogDensityProblems
import LogDensityProblems: logdensity, logdensity_and_gradient

# ### snippet 9.6

function HMC(model, grad, epsilon, L, current_q)
  q = current_q
  p = rand(Normal(0, 1), length(q)) # random flick - p is momentum
  #p = [0.96484367, -0.06740435]
  current_p = p
  # Make a half step for momentum at the beginning
  v, g = logdensity_and_gradient(grad, q)
  value_U = -v; grad_U = -g
  p = p - epsilon .* grad_U / 2.0
  # Initialize bookkeeping
  qtraj = zeros(L+1, length(q)+3)
  qtraj[1, :] = [q[1], q[2], value_U, grad_U[1], grad_U[2]]

  ptraj = zeros(L+1, length(q))
  ptraj[1, :] = p

  # ### snippet 9.7

  # Alternate full steps for position and momentum

  for i in 1:L
    # Full position step
    q = q + epsilon .* p
    # Full step for momentum,, except for last step
    if i !== L
      v, g = logdensity_and_gradient(grad, q)
      value_U = -v; grad_U = -g
      p = p - epsilon .* grad_U
      ptraj[i+1, :] = p
    end
    qtraj[i+1, :] = [q[1], q[2], value_U, grad_U[1], grad_U[2]]
  end

  # ### snippet 9.8

  # Make a halfstep for momentum at the end

  v, g = logdensity_and_gradient(grad, q)
  value_U = -v; grad_U = -g
  p = p - epsilon .* grad_U / 2.0
  ptraj[L+1, :] = p
  # Negate momentum to make proposal symmatric
  p = -p
  # Evaluate potential and kinetic energies ate beginning and end
  current_U = -logdensity(grad, [current_q[1], current_q[2]])
  current_K = sum(current_p .^ 2) / 2
  proposed_U = -logdensity(grad, [q[1], q[2]])
  proposed_K = sum(p .^ 2) / 2
  dH = proposed_U + proposed_K - current_U - current_K
  # Accept or reject the state at the end of trajectory
  # Return either position at the end or initial position
  local accept = 0
  local new_q
  if rand(Uniform(0, 1), 1)[1] < exp(dH)
    new_q = q # Accept
    accept = 1
  else
    new_q = current_q # Reject
  end

  (new_q, ptraj, qtraj, accept, dH)

end


text:
using ..NEPCore
using ..NEPTypes
using ..LinSolvers
export mslp

"""
     mslp([eltype],nep::NEP;[errmeasure,][tol,][maxit,][λ,][logger,][eigsolvertype::Type])

Runs the method of successive linear problems. The  method requires the solution of a
generalized eigenvalue problem in every iteration. The method used for the eigenvalue
computation is specified in eigsolvertype.
See [`augnewton`](@ref) for other parameters.

# Example
Create a rational NEP using a [`SPMF_NEP`](@ref).
```julia-repl
julia> eye=Matrix{Float64}(I,3,3);
julia> Av=[ones(3,3),eye,triu(ones(3,3))];
julia> fv=[S-> S, S -> S^2, S->inv(S-one(S)*10)];
julia> nep=SPMF_NEP(Av,fv);
julia> (λ,v)=mslp(nep);
julia> compute_Mlincomb(nep,λ,v)
3-element Array{Complex{Float64},1}:
 -1.38778e-17+1.65715e-18im
 -5.55112e-17+1.30633e-17im
 -4.16334e-17-1.54436e-17im
```

# References
* A. Ruhe, Algorithms for the nonlinear eigenvalue problem, SIAM J. Numer. Anal. 10 (1973) 674-689

"""
mslp(nep::NEP;params...)=mslp(ComplexF64,nep;params...)
function mslp(::Type{T},
                 nep::NEP;
                 errmeasure::ErrmeasureType = DefaultErrmeasure(nep),
                 tol::Real=eps(real(T))*100,
                 maxit::Integer=100,
                 λ::Number=zero(T),
                 logger=0,
                 eigsolvertype::Type=DefaultEigSolver) where {T<:Number}

    @parse_logger_param!(logger)

    # Ensure types λ is of type T
    λ::T = T(λ)

    # Allocate memory for the eigenvector approximation
    v = zeros(T,size(nep,1));
    d = zeros(T,1);

    err=Inf;

    # Main loop
    for k=1:maxit

        # solve generalized eigenvalue problem
        solver=eigsolvertype(compute_Mder(nep,λ,0),compute_Mder(nep,λ,1));

        # This will throw an error if the eigenvector is not of correct type
        d[:],v[:] = eig_solve(solver,target=0,nev=1);

        # update eigenvalue
        λ += -d[1] # This will throw an error if the eigval update d is not of correct type

        # Normalize
        normalize!(v)

        # Checck for convergence
        err=estimate_error(errmeasure,λ,v)
        push_iteration_info!(logger,k,err=err,λ=λ);

        if (err< tol)
            return (λ,v)
        end

    end

    msg="Number of iterations exceeded. maxit=$(maxit)."
    throw(NoConvergenceException(λ,v,err,msg))
end


text:
using Documenter, FlexPlan

Documenter.makedocs(
    modules = FlexPlan,
    format = Documenter.HTML(),
    sitename = "FlexPlan",
    authors = "Hakan Ergun, Matteo Rossini, Damien Lapage, Iver Bakken Sperstad, Espen Flo Bødal, Marco Rossi, Merkebu Zenebe Degefa, Reinhilde D'hulst",
    pages = [
        "Home" => "index.md"
        "Manual" => [
            "Getting started" => "quickguide.md"
            "Example scripts" => "example_scripts.md"
        ]
        "Library" => [
            "Problem types" => "problem_types.md"
            "Network formulations" => "network_formulations.md"
            "Multi - period, multi-stage modelling" => [
                "Modelling assumptions" => "modeling_assumptions.md"
                "Model dimensions" => "dimensions.md"
            ]
            "Data model" => "data_model.md"
        ]
    ]
)

Documenter.deploydocs(
    repo = "github.com/Electa-Git/FlexPlan.jl.git"
)


text:
export projector, print_blocktree
using InteractiveUtils

"""
    projector(x)

Return projector on `0` or projector on `1`.
"""
projector(x) = code == 0 ? mat(P0) : mat(P1)

"""
    print_subtypetree(::Type[, level=1, indent=4])

Print subtype tree, `level` specify the depth of the tree.
"""
function print_subtypetree(t::Type, level = 1, indent = 4)
    level == 1 && println(t)
    for s in subtypes(t)
        println(join(fill(" ", level * indent)) * string(s))
        print_subtypetree(s, level + 1, indent)
    end
end


text:
# ---
# title: 1470. Shuffle the Array
# id: problem1470
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Array
# link: <https://leetcode.com/problems/shuffle-the-array/description/>
# hidden: true
# ---
#
# Given the array `nums` consisting of `2n` elements in the form
# `[x1,x2,...,xn,y1,y2,...,yn]`.
#
# _Return the array in the form_ `[x1,y1,x2,y2,...,xn,yn]`.
#
#
#
# **Example 1:**
#
#
#
#     Input: nums = [2,5,1,3,4,7], n = 3
#     Output: [2,3,5,4,1,7]
#     Explanation: Since x1=2, x2=5, x3=1, y1=3, y2=4, y3=7 then the answer is [2,3,5,4,1,7].
#
#
# **Example 2:**
#
#
#
#     Input: nums = [1,2,3,4,4,3,2,1], n = 4
#     Output: [1,4,2,3,3,2,4,1]
#
#
# **Example 3:**
#
#
#
#     Input: nums = [1,1,2,2], n = 2
#     Output: [1,2,1,2]
#
#
#
#
# **Constraints:**
#
#   * `1 <= n <= 500`
#   * `nums.length == 2n`
#   * `1 <= nums[i] <= 10^3`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
"""
    assemble_model(;
        aerodynamic_model = nothing,
        structural_model = nothing,
        dynamics_model = nothing,
        control_surface_model = nothing,
        controller_model = nothing,
        coupling_model = nothing)

Assemble a coupled model from a collection of models.
"""
function assemble_model(;
    aerodynamic_model = nothing,
    structural_model = nothing,
    dynamics_model = nothing,
    control_surface_model = nothing,
    controller_model = nothing,
    coupling_model = nothing)

    # construct tuple of provided models
    models = ()
    for model in (aerodynamic_model, structural_model, dynamics_model,
        control_surface_model, controller_model)

        if !isnothing(model)
            models = (models..., model)
        end

    end

    # define submodels using provided models
    submodels = Submodel.(models)

    # use default coupling if a coupling model is not provided
    coupling = isnothing(coupling_model) ? Coupling(models...) : coupling_model

    # construct coupled model
    return CoupledModel(submodels, coupling)
end

"""
    assemble_states(model;
        aerodynamic_states = nothing,
        structural_states = nothing,
        dynamics_states = nothing,
        control_surface_states = nothing,
        controller_states = nothing)

Assemble the state vector for a coupled model using the specified state variables.  State
variables are initialized with zeros.
"""
function assemble_states(model::CoupledModel; kwargs...)

    x = zeros(number_of_states(model))

    return assemble_states!(x, model; kwargs...)
end

"""
    assemble_states!(x, model;
        aerodynamic_states = nothing,
        structural_states = nothing,
        dynamics_states = nothing,
        control_surface_states = nothing,
        controller_states = nothing)

In-place version of [`assemble_states`](@ref).  Only replaces the values corresponding to
the specified state variables.
"""
function assemble_states!(x0, model::CoupledModel;
    aerodynamic_states = nothing,
    structural_states = nothing,
    dynamics_states = nothing,
    control_surface_states = nothing,
    controller_states = nothing)

    # construct tuple of provided state variables
    coupled_states = ()
    for states in (aerodynamic_states, structural_states, dynamics_states,
        control_surface_states, controller_states)

        if !isnothing(states)
            coupled_states = (coupled_states..., states)
        end

    end

    # populate state vector
    for i = 1:length(coupled_states)
        set_states!(x0, model, i; coupled_states[i]...)
    end

    return x0
end

"""
    assemble_parameters(model;
        aerodynamic_parameters = nothing,
        structural_parameters = nothing,
        dynamics_parameters = nothing,
        control_surface_parameters = nothing,
        controller_parameters = nothing,
        additional_parameters = nothing)

Assemble the parameter vector for a coupled model using the specified parameter variables.
"""
function assemble_parameters(model::CoupledModel; kwargs...)

    x = zeros(number_of_parameters(model))

    return assemble_parameters!(x, model; kwargs...)
end

"""
    assemble_parameters!(p, model;
        aerodynamic_parameters = nothing,
        structural_parameters = nothing,
        dynamics_parameters = nothing,
        control_surface_parameters = nothing,
        controller_parameters = nothing,
        additional_parameters = nothing)

In-place version of [`assemble_parameters`](@ref). Only replaces the values corresponding
to the specified parameters.
"""
function assemble_parameters!(p, model::CoupledModel;
    aerodynamic_parameters = nothing,
    structural_parameters = nothing,
    dynamics_parameters = nothing,
    control_surface_parameters = nothing,
    controller_parameters = nothing,
    additional_parameters = nothing)

    # define coupled parameters
    coupled_parameters = ()
    for parameters in (aerodynamic_parameters, structural_parameters, dynamics_parameters,
        control_surface_parameters, controller_parameters)

        if !isnothing(parameters)
            coupled_parameters = (coupled_parameters..., parameters)
        end

    end

    # populate parameter vector
    for i = 1:length(coupled_parameters)
        set_parameters!(p, model, i; coupled_parameters[i]...)
    end
    set_additional_parameters!(p, model; additional_parameters...)

    return p
end

"""
    linearize(model, x, p; kwargs...)

Calculate the jacobian of the residual function for coupled model `model` with respect to
the state variables and their rates.
"""
function linearize(model::CoupledModel, x, p; dx = FillArrays.Zeros(x), t = 0,
    y = get_coupling_inputs(model, dx, x, p, t))

    return _linearize(inplaceness(model), model, dx, x, y, p, t)
end

function _linearize(::OutOfPlace, model, args...)
    K = Array(get_state_jacobian(model, args...))
    M = Array(get_rate_jacobian(model, args...))
    return K, M
end

function _linearize(::InPlace, model, args...)
    K = get_state_jacobian(model, args...)
    M = get_rate_jacobian(model, args...)
    return K, M
end

"""
    get_eigen(model::TM, K, M; kwargs...)

Return the eigenvalues, left eigenvector matrix, and right eigenvector matrix
corresponding to the model.

For in-place models, the number of eigenvalues to compute may be specified using
the `nev` keyword argument.
"""
function get_eigen(model, K, M; kwargs...)

    return _get_eigen(inplaceness(model), K, M; kwargs...)
end

function _get_eigen(::OutOfPlace, K, M)
    A = -K
    B = M
    E = eigen(A, B) # eigenvalue decomposition
    λ = eigvals(E) # eigenvalues
    V = eigvecs(E) # right eigenvector matrix
    U = I/(A*V) # left eigenvector matrix
    return λ, U, V
end

function _get_eigen(::InPlace, K, M; nev=min(20, number_of_states(model)))

    # calculate state and rate jacobians
    A = -K
    B = M

    # construct linear map
    nx = size(A, 1)
    TF = promote_type(eltype(A), eltype(B))
    Afact = lu(A)
    f! = (b, x) -> ldiv!(b, Afact, B*x)
    fc! = (b, x) -> mul!(b, B', Afact'\x)
    Abar = LinearMap{TF}(f!, fc!, nx, nx; ismutating=true)

    # compute eigenvalues and eigenvectors
    λ, V = partialeigen(partialschur(Abar; nev=min(nx,nev), which=LM())[1])

    # sort eigenvalues by magnitude
    perm = sortperm(λ, by=(λ)->(abs(λ),imag(λ)), rev=true)
    λ = λ[perm[1:nev]]
    V = V[:,perm[1:nev]]

    # eigenvalues are actually 1/λ, no modification necessary for eigenvectors
    λ .= 1 ./ λ

    # also return left eigenvectors
    U = GXBeam.left_eigenvectors(A, -B, λ, V)

    return λ, U, V
end

"""
    ODEFunction(model::Aeroelasticity.CoupledModel, p=nothing)

Construct an `ODEFunction` for a coupled model which may be used with
`DifferentialEquations`.
"""
function SciMLBase.ODEFunction(model::CoupledModel, p=nothing)

    # always construct in-place models, always compile
    iip = true
    compile = true

    # construct coupling function
    fy = ode_coupling_function(model)

    # construct rate function
    f = ode_rate_function(model, fy)

    # construct mass matrix
    mass_matrix = ode_mass_matrix(model, fy)

    # construct time gradient function
    tgrad = ode_time_gradient(model, fy)

    # construct jacobian function
    jac = ode_state_jacobian(model, fy)

    # construct parameter jacobian function
    paramjac = ode_parameter_jacobian(model, fy)

    # return ODE function
    return ODEFunction{iip, compile}(f; mass_matrix, tgrad, jac, paramjac)
end

function ode_coupling_function(model::CoupledModel)

    # check whether cached variables should be used
    if isinplace(model)
        # model is inplace so use cached variables if possible

        # problem dimensions
        Nx = number_of_states(model)
        Ny = number_of_inputs(model)
        Np = number_of_parameters(model)

        # cached variables
        xcache = fill(NaN, Nx)
        ycache = fill(NaN, Ny)
        pcache = fill(NaN, Np)
        tcache = fill(NaN)

        # state rate vector
        dx = FillArrays.Zeros(x)

        # coupling function
        fy = (x, p, t) -> begin
            # check if we can use the cache variables (no custom types)
            if promote_type(eltype(x), eltype(p), typeof(t)) <: Float64
                # update the cache variables (if necessary)
                if (x != xcache) && (p != pcache) && (t != tcache[])
                    xcache .= x
                    ycache .= get_coupling_inputs!(ycache, model, dx, x, p, t)
                    pcache .= p
                    tcache .= t
                end
                # set inputs to the cached inputs
                y = ycache
            else
                # calculate inputs (out-of-place to accomodate custom type)
                y = get_coupling_inputs(model, dx, x, p, t)
            end
            # return cached or computed inputs
            return y
        end
    else
        # model is out of place so calculate coupling inputs directly
        fy = (x, p, t) -> get_coupling_inputs(model, FillArrays.Zeros(x), x, p, t)
    end

    # return coupling function
    return fy
end

function ode_rate_function(model::CoupledModel, fy)

    f = (dx, x, p, t) -> begin
        get_residual!(dx, model, FillArrays.Zeros(x), x, fy(x, p, t), p, t)
        return dx .*= -1
    end

    return f
end

function ode_mass_matrix(model::CoupledModel, fy, p=nothing)

    submodels = model.submodels
    coupling = model.coupling

    if all(isidentity.(getproperty.(submodels, :ratejac))) && iszero(coupling.ratejac)
        # mass matrix is the identity matrix
        mass_matrix = I
    elseif all(isinvariant.(getproperty.(submodels, :ratejac))) && isinvariant(coupling.ratejac)
        # mass matrix is independent of all inputs
        mass_matrix = get_rate_jacobian(model)
    elseif !isnothing(p) && all(isconstant.(getproperty.(submodels, :ratejac))) &&
        isconstant(coupling.ratejac)
        # mass matrix only depends on the parameter vector
        mass_matrix = get_rate_jacobian(model, p)
    else
        # mass matrix defined as a function of state variables, parameters, and time

        # problem dimensions
        Nx = number_of_states(model)
        Ny = number_of_inputs(model)

        # state rate vector
        dx = FillArrays.Zeros(Nx)

        # initialize mass matrix
        M = zeros(Nx, Nx)

        # construct update function

        # cached variables
        drdy_cache = fill(NaN, Nx, Ny)
        dyddx_cache = fill(NaN, Ny, Nx)

        # update function
        update_func = (M, x, p, t) -> begin
            # check if we can use the cache variables (no custom types)
            if promote_type(eltype(x), eltype(p), typeof(t)) <: Float64
                # use cache variables
                get_rate_jacobian!(M, model, dx, x, fy(x, p, t), p, t;
                    drdy_cache = drdy_cache, dyddx_cache = dyddx_cache)
            else
                # don't use cache variables (to accomodate custom types)
                get_rate_jacobian!(M, model, dx, x, fy(x, p, t), p, t)
            end
        end

        mass_matrix = DiffEqArrayOperator(M; update_func)
    end

    return mass_matrix
end

function ode_time_gradient(model::CoupledModel, fy, p=nothing)

    submodels = model.submodels
    coupling = model.coupling

    if all(iszero.(getproperty.(submodels, :tgrad))) && iszero(coupling.tgrad)
        # time gradient is a zero vector
        tgrad = (dT, x, p, t) -> dT .= 0
    elseif all(isinvariant.(getproperty.(submodels, :tgrad))) && isinvariant(coupling.tgrad)
        # time gradient is independent of all inputs
        dT0 = get_time_gradient(model)
        tgrad = (dT, x, p, t) -> dT .= dT0
    elseif !isnothing(p) && all(isconstant.(getproperty.(submodels, :tgrad))) &&
        isconstant(coupling.tgrad)
        # time gradient only depends on the parameter vector
        dT0 = get_rate_jacobian(model, p)
        tgrad = (dT, x, p, t) -> dT .= dT0
    else
        # time gradient defined as a function of state variables, parameters, and time

        # problem dimensions
        Nx = number_of_states(model)
        Ny = number_of_inputs(model)

        # cached variables
        drdy_cache = fill(NaN, Nx, Ny),
        dyddx_cache = fill(NaN, Ny, Nx)

        # state rate vector
        dx = FillArrays.Zeros(x)

        # update function
        tgrad = (dT, x, p, t) -> begin
            # check if we can use the cache variables (no custom types)
            if promote_type(eltype(x), eltype(p), typeof(t)) <: Float64
                # use cache variables
                get_time_gradient!(dT, model, dx, x, fy(x, p, t), p, t;
                    drdy_cache = drdy_cache, dyddx_cache = dyddx_cache)
            else
                # don't use cache variables (to accomodate custom types)
                get_time_gradient!(dT, model, dx, x, fy(x, p, t), p, t)
            end
        end
    end

    return tgrad
end

function ode_state_jacobian(model::CoupledModel, fy, p=nothing)

    submodels = model.submodels
    coupling = model.coupling

    if all(iszero.(getproperty.(submodels, :statejac))) && iszero(coupling.statejac)
        # state jacobian is a zero vector
        jac = (J, x, p, t) -> J .= 0
    elseif all(isinvariant.(getproperty.(submodels, :statejac))) && isinvariant(coupling.statejac)
        # state jacobian is independent of all inputs
        J0 = get_state_jacobian(model)
        jac = (J, x, p, t) -> J .= J0
    elseif !isnothing(p) && all(isconstant.(getproperty.(submodels, :statejac))) &&
        isconstant(coupling.statejac)
        # state jacobian only depends on the parameter vector
        J0 = get_state_jacobian(model, p)
        jac = (J, x, p, t) -> J .= J0
    else
        # state jacobian defined as a function of state variables, parameters, and time

        # problem dimensions
        Nx = number_of_states(model)
        Ny = number_of_inputs(model)

        # state rate vector
        dx = FillArrays.Zeros(Nx)

        # construct state jacobian function

        # cached variables
        drdy_cache = fill(NaN, Nx, Ny)
        dyddx_cache = fill(NaN, Ny, Nx)

        # update function
        jac = (J, x, p, t) -> begin
            # check if we can use the cache variables (no custom types)
            if promote_type(eltype(x), eltype(p), typeof(t)) <: Float64
                # use cache variables
                get_state_jacobian!(J, model, dx, x, fy(x, p, t), p, t;
                    drdy_cache = drdy_cache, dyddx_cache = dyddx_cache)
            else
                # don't use cache variables (to accomodate custom types)
                get_state_jacobian!(J, model, dx, x, fy(x, p, t), p, t)
            end
        end
    end

    return jac
end

function ode_parameter_jacobian(model::CoupledModel, fy, p=nothing)

    submodels = model.submodels
    coupling = model.coupling

    if all(iszero.(getproperty.(submodels, :paramjac))) && isidentity(coupling.paramjac)
        # parameter jacobian is a zero matrix
        jac = (pJ, x, p, t) -> pJ .= 0
    elseif all(isinvariant.(getproperty.(submodels, :paramjac))) && isinvariant(coupling.paramjac)
        # parameter jacobian is independent of all inputs
        pJ0 = get_parameter_jacobian(model)
        paramjac = (pJ, x, p, t) -> pJ .= pJ0
    elseif !isnothing(p) && all(isconstant.(getproperty.(submodels, :paramjac))) &&
        isconstant(coupling.paramjac)
        # time gradient only depends on the parameter vector
        pJ0 = get_parameter_jacobian(model, p)
        paramjac = (pJ, x, p, t) -> pJ .= pJ0
    else
        # construct time gradient function

        # problem dimensions
        Nx = number_of_states(model)
        Ny = number_of_inputs(model)

        # cached variables
        drdy_cache = fill(NaN, Nx, Ny)
        dyddx_cache = fill(NaN, Ny, Nx)

        # state rate vector
        dx = FillArrays.Zeros(Nx)

        # update function
        pjac = (dT, x, p, t) -> begin
            # check if we can use the cache variables (no custom types)
            if promote_type(eltype(x), eltype(p), typeof(t)) <: Float64
                # use cache variables
                get_parameter_jacobian!(pJ, model, dx, x, fy(x, p, t), p, t;
                    drdy_cache = drdy_cache, dyddx_cache = dyddx_cache)
            else
                # don't use cache variables (to accomodate custom types)
                get_parameter_jacobian!(pJ, model, dx, x, fy(x, p, t), p, t)
            end
        end

    end

    return pjac
end

"""
    DAEFunction(model::Aeroelasticity.CoupledModel, p=nothing)

Construct an `DAEFunction` for a coupled model which may be used with
`DifferentialEquations`.
"""
function SciMLBase.DAEFunction(model::CoupledModel, p=nothing)

    # always construct in-place models, always compile
    iip = true
    compile = true

    # construct coupling function
    fy = dae_coupling_function(model)

    # construct rate function
    f = dae_residual_function(model, fy)

    # return ODE function
    return DAEFunction{iip, compile}(f)
end

function dae_coupling_function(model::CoupledModel)

    # check whether cached variables should be used
    if isinplace(model)
        # model is inplace so use cached variables if possible

        # problem dimensions
        Nx = number_of_states(model)
        Ny = number_of_inputs(model)
        Np = number_of_parameters(model)

        # cached variables
        dxcache = fill(NaN, Nx)
        xcache = fill(NaN, Nx)
        ycache = fill(NaN, Ny)
        pcache = fill(NaN, Np)
        tcache = fill(NaN)

        # coupling function
        fy = (dx, x, p, t) -> begin
            # check if we can use the cache variables (no custom types)
            if promote_type(eltype(dx), eltype(x), eltype(p), typeof(t)) <: Float64
                # update the cache variables (if necessary)
                if (dx != dxcache) && (x != xcache) && (p != pcache) && (t != tcache[])
                    dxcache .= dx
                    xcache .= x
                    ycache .= get_coupling_inputs!(ycache, model, dx, x, p, t)
                    pcache .= p
                    tcache .= t
                end
                # set inputs to the cached inputs
                y = ycache
            else
                # calculate inputs (out-of-place to accomodate custom type)
                y = get_coupling_inputs(model, dx, x, p, t)
            end
            # return cached or computed inputs
            return y
        end
    else
        # model is out of place so calculate coupling inputs directly
        fy = (dx, x, p, t) -> get_coupling_inputs(model, dx, x, p, t)
    end

    # return coupling function
    return fy
end

function dae_residual_function(model::CoupledModel, fy)

    return (resid, dx, x, p, t) -> get_residual!(resid, model, dx, x, fy(dx, x, p, t), p, t)
end


text:
struct Reaction{S <: Variable, T <: Number}
    rate
    substrates::Vector{Operation}
    products::Vector{Operation}
    substoich::Vector{T}
    prodstoich::Vector{T}
    netstoich::Vector{Pair{S,T}}
    only_use_rate::Bool
end

function Reaction(rate, subs, prods, substoich, prodstoich;
                  netstoich=nothing, only_use_rate=false, kwargs...)

    subsv  = isnothing(subs) ? Vector{Operation}() : subs
    prodsv = isnothing(prods) ? Vector{Operation}() : prods
    ns = isnothing(netstoich) ? get_netstoich(subsv, prodsv, substoich, prodstoich) : netstoich
    Reaction(rate, subsv, prodsv, substoich, prodstoich, ns, only_use_rate)
end

# three argument constructor assumes stoichiometric coefs are one and integers
function Reaction(rate, subs, prods; kwargs...)

    sstoich = isnothing(subs) ? Int[] : ones(Int,length(subs))
    pstoich = isnothing(prods) ? Int[] : ones(Int,length(prods))
    Reaction(rate, subs, prods, sstoich, pstoich; kwargs...)
end

# calculates the net stoichiometry of a reaction as a vector of pairs (sub,substoich)
function get_netstoich(subs, prods, sstoich, pstoich)
    # stoichiometry as a Dictionary
    nsdict = Dict{Variable,eltype(sstoich)}(sub.op => -sstoich[i] for (i,sub) in enumerate(subs))
    for (i,p) in enumerate(prods)
        coef = pstoich[i]
        prod = p.op
        @inbounds nsdict[prod] = haskey(nsdict, prod) ? nsdict[prod] + coef : coef
    end

    # stoichiometry as a vector
    ns = [el for el in nsdict if el[2] != zero(el[2])]

    ns
end

struct ReactionSystem <: AbstractSystem
    eqs::Vector{Reaction}
    iv::Variable
    states::Vector{Variable}
    ps::Vector{Variable}
    name::Symbol
    systems::Vector{ReactionSystem}
end

function ReactionSystem(eqs, iv, species, params; systems = ReactionSystem[],
                                                  name = gensym(:ReactionSystem))

    ReactionSystem(eqs, iv, convert.(Variable,species), convert.(Variable,params),
                   name, systems)
end

# Calculate the ODE rate law
function oderatelaw(rx)
    @unpack rate, substrates, substoich, only_use_rate = rx
    rl = rate
    if !only_use_rate
        coef = one(eltype(substoich))
        for (i,stoich) in enumerate(substoich)
            coef *= factorial(stoich)
            rl   *= isone(stoich) ? substrates[i] : substrates[i]^stoich
        end
        (!isone(coef)) && (rl /= coef)
    end
    rl
end

function assemble_drift(rs)
    D   = Differential(rs.iv())
    eqs = [D(x(rs.iv())) ~ 0 for x in rs.states]
    species_to_idx = Dict((x => i for (i,x) in enumerate(rs.states)))

    for rx in rs.eqs
        rl = oderatelaw(rx)
        for (spec,stoich) in rx.netstoich
            i = species_to_idx[spec]
            if iszero(eqs[i].rhs)
                signedrl = (stoich > zero(stoich)) ? rl : -rl
                rhs      = isone(abs(stoich)) ? signedrl : stoich * rl
            else
                Δspec = isone(abs(stoich)) ? rl : abs(stoich) * rl
                rhs   = (stoich > zero(stoich)) ? (eqs[i].rhs + Δspec) : (eqs[i].rhs - Δspec)
            end
            eqs[i] = Equation(eqs[i].lhs, rhs)
        end
    end
    eqs
end

function assemble_diffusion(rs)
    eqs = Expression[Constant(0) for x in rs.states, y in rs.eqs]
    species_to_idx = Dict((x => i for (i,x) in enumerate(rs.states)))

    for (j,rx) in enumerate(rs.eqs)
        rlsqrt = sqrt(oderatelaw(rx))
        for (spec,stoich) in rx.netstoich
            i            = species_to_idx[spec]
            signedrlsqrt = (stoich > zero(stoich)) ? rlsqrt : -rlsqrt
            eqs[i,j]     = isone(abs(stoich)) ? signedrlsqrt : stoich * rlsqrt
        end
    end
    eqs
end

function Base.convert(::Type{<:ODESystem},rs::ReactionSystem)
    eqs = assemble_drift(rs)
    ODESystem(eqs,rs.iv,rs.states,rs.ps,name=rs.name,
              systems=convert.(ODESystem,rs.systems))
end

function Base.convert(::Type{<:SDESystem},rs::ReactionSystem)
    eqs = assemble_drift(rs)
    noiseeqs = assemble_diffusion(rs)
    SDESystem(eqs,noiseeqs,rs.iv,rs.states,rs.ps,
              name=rs.name,systems=convert.(SDESystem,rs.systems))
end


text:
################################################################################

"""

    shifter(inputAr)

# Description
Reshape signals for autoencoder.

See also: [`reshifter`](@ref)
"""
function shifter(inputAr)
  outAr = [
    Array{Float64}(
      undef,
      prod(size(inputAr)[1:2])
    )
    for _ ∈ 1:size(inputAr, 3)
  ]
  for ι ∈ 1:size(inputAr, 3)
    outAr[ι] = vec(inputAr[:, :, ι])
  end
  return outAr
end

################################################################################

"""

    reshifter(inputAr;
    outSize = length(inputAr[1]))

# Description
Reshape signals from autoencoder.

See also: [`shifter`](@ref)
"""
function reshifter(inputAr; outSize = length(inputAr[1]))
  ψ = length(inputAr[1]) / outSize |> π -> convert(Int64, π)
  outAr  = Array{Float64}(
    undef,
    ψ,
    outSize,
    length(inputAr)
  )

  for ι ∈ 1:length(inputAr)
    outAr[:, :, ι] = inputAr[ι] |> π -> reshape(π, ψ, outSize)
  end

  return outAr
end

################################################################################

max_stars_repo_name: DanielRivasMD/MindReader

max_issues_repo_name: DanielRivasMD/MindReader

max_forks_repo_name: DanielRivasMD/MindReader