text:
####
#### Output writer utilities
####

ext(fw::AbstractOutputWriter) = throw("Extension for $(typeof(fw)) is not implemented.")

# When saving stuff to disk like a JLD2 file, `saveproperty!` is used, which
# converts Julia objects to language-agnostic objects.
saveproperty!(file, location, p::Number)        = file[location] = p
saveproperty!(file, location, p::AbstractRange) = file[location] = collect(p)
saveproperty!(file, location, p::AbstractArray) = file[location] = Array(p)
saveproperty!(file, location, p::AbstractField)         = file[location] = Array(p.data.parent)
saveproperty!(file, location, p::Function) = @warn "Cannot save Function property into $location"
saveproperty!(file, location, p) = [saveproperty!(file, location * "/$subp", getproperty(p, subp))
                                        for subp in propertynames(p)]

# Special saveproperty! so boundary conditions are easily readable outside julia.
function saveproperty!(file, location, cbcs::CoordinateBoundaryConditions)
    for endpoint in propertynames(cbcs)
        endpoint_bc = getproperty(cbcs, endpoint)
        if isa(endpoint_bc.condition, Function)
            @warn "$field.$coord.$endpoint boundary is of type Function and cannot be saved to disk!"
            file["boundary_conditions/$field/$coord/$endpoint/type"] = string(bctype(endpoint_bc))
            file["boundary_conditions/$field/$coord/$endpoint/condition"] = missing
        else
            file["boundary_conditions/$field/$coord/$endpoint/type"] = string(bctype(endpoint_bc))
            file["boundary_conditions/$field/$coord/$endpoint/condition"] = endpoint_bc.condition
        end
    end
end

# Special savepropety! for AB2 time stepper struct used by the checkpointer so
# it only saves the fields and not the tendency BCs or χ value (as they can be
# constructed by the `Model` constructor).
function saveproperty!(file, location, ts::AdamsBashforthTimestepper)
    saveproperty!(file, location * "/Gⁿ", ts.Gⁿ)
    saveproperty!(file, location * "/G⁻", ts.G⁻)
end

saveproperties!(file, structure, ps) = [saveproperty!(file, "$p", getproperty(structure, p)) for p in ps]

# When checkpointing, `serializeproperty!` is used, which serializes objects
# unless they need to be converted (basically CuArrays only).
serializeproperty!(file, location, p) = file[location] = p
serializeproperty!(file, location, p::Union{AbstractArray, AbstractField}) = saveproperty!(file, location, p)
serializeproperty!(file, location, p::Function) = @warn "Cannot serialize Function property into $location"

serializeproperties!(file, structure, ps) = [serializeproperty!(file, "$p", getproperty(structure, p)) for p in ps]

# Don't check arrays because we don't need that noise.
has_reference(T, ::AbstractArray{<:Number}) = false

# These two conditions are true, but should not necessary.
has_reference(::Type{Function}, ::AbstractField) = false
has_reference(::Type{T}, ::NTuple{N, <:T}) where {N, T} = true

"""
    has_reference(has_type, obj)

Check (or attempt to check) if `obj` contains, somewhere among its
subfields and subfields of fields, a reference to an object of type
`has_type`. This function doesn't always work.
"""
function has_reference(has_type, obj)
    if typeof(obj) <: has_type
        return true
    elseif applicable(iterate, obj) && length(obj) > 1
        return any([has_reference(has_type, elem) for elem in obj])
    elseif applicable(propertynames, obj) && length(propertynames(obj)) > 0
        return any([has_reference(has_type, getproperty(obj, p)) for p in propertynames(obj)])
    else
        return typeof(obj) <: has_type
    end
end

####
####  JLD2 output writer
####

"""
    JLD2OutputWriter{F, I, O, IF, IN, KW} <: AbstractOutputWriter

An output writer for writing to JLD2 files.
"""
mutable struct JLD2OutputWriter{F, I, O, IF, IN, KW} <: AbstractOutputWriter
        filepath :: String
         outputs :: O
        interval :: I
       frequency :: F
            init :: IF
       including :: IN
        previous :: Float64
            part :: Int
    max_filesize :: Float64
           async :: Bool
           force :: Bool
         verbose :: Bool
         jld2_kw :: KW
end

noinit(args...) = nothing

"""
    JLD2OutputWriter(model, outputs; interval=nothing, frequency=nothing, dir=".",
                     prefix="", init=noinit, including=[:grid, :coriolis, :buoyancy, :closure],
                     part=1, max_filesize=Inf, force=false, async=false, verbose=false)

Construct a `JLD2OutputWriter` that writes `label, func` pairs in `outputs` (which can be a `Dict` or `NamedTuple`)
to a JLD2 file, where `label` is a symbol that labels the output and `func` is a function of the form `func(model)`
that returns the data to be saved.

Keyword arguments
=================

    - `frequency::Int`   : Save output every `n` model iterations.
    - `interval::Int`    : Save output every `t` units of model clock time.
    - `dir::String`      : Directory to save output to. Default: "." (current working
                           directory).
    - `prefix::String`   : Descriptive filename prefixed to all output files. Default: "".
    - `init::Function`   : A function of the form `init(file, model)` that runs when a JLD2
                           output file is initialized. Default: `noinit(args...) = nothing`.
    - `including::Array` : List of model properties to save with every file. By default, the
                           grid, equation of state, coriolis parameters, buoyancy parameters,
                           and turbulence closure parameters are saved.
    - `part::Int`        : The starting part number used if `max_filesize` is finite.
                           Default: 1.
    - `max_filesize::Int`: The writer will stop writing to the output file once the file size
                           exceeds `max_filesize`, and write to a new one with a consistent
                           naming scheme ending in `part1`, `part2`, etc. Defaults to `Inf`.
    - `force::Bool`      : Remove existing files if their filenames conflict. Default: `false`.
    - `async::Bool`      : Write output asynchronously. Default: `false`.
    - `verbose::Bool`    : Log what the output writer is doing with statistics on compute/write
                           times and file sizes. Default: `false`.
    - `jld2_kw::Dict`    : Dict of kwargs to be passed to `jldopen` when data is written.
"""
function JLD2OutputWriter(model, outputs; interval=nothing, frequency=nothing, dir=".", prefix="",
                          init=noinit, including=[:grid, :coriolis, :buoyancy, :closure],
                          part=1, max_filesize=Inf, force=false, async=false, verbose=false,
                          jld2_kw=Dict{Symbol, Any}())

    validate_interval(interval, frequency)

    mkpath(dir)
    filepath = joinpath(dir, prefix * ".jld2")
    force && isfile(filepath) && rm(filepath, force=true)

    jldopen(filepath, "a+"; jld2_kw...) do file
        init(file, model)
        saveproperties!(file, model, including)
    end

    return JLD2OutputWriter(filepath, outputs, interval, frequency, init, including,
                            0.0, part, max_filesize, async, force, verbose, jld2_kw)
end

function write_output(model, fw::JLD2OutputWriter)
    verbose = fw.verbose
    verbose && @info @sprintf("Calculating JLD2 output %s...", keys(fw.outputs))
    tc = Base.@elapsed data = Dict((name, f(model)) for (name, f)
                                    in zip(keys(fw.outputs), values(fw.outputs)))
    verbose && @info "Calculation time: $(prettytime(tc))"

    iter = model.clock.iteration
    time = model.clock.time

    filesize(fw.filepath) >= fw.max_filesize && start_next_file(model, fw)

    path = fw.filepath
    verbose && @info "Writing JLD2 output $(keys(fw.outputs)) to $path..."
    t0, sz = time_ns(), filesize(path)

    if fw.async
        @async remotecall(jld2output!, 2, path, iter, time, data, fw.jld2_kw)
    else
        jld2output!(path, iter, time, data, fw.jld2_kw)
    end

    t1, newsz = time_ns(), filesize(path)

    verbose && @info @sprintf("Writing done: time=%s, size=%s, Δsize=%s",
                              prettytime((t1-t0)/1e9), pretty_filesize(newsz), pretty_filesize(newsz-sz))

    return
end

"""
    jld2output!(path, iter, time, data, kwargs)

Write the (name, value) pairs in `data`, including the simulation
`time`, to the JLD2 file at `path` in the `timeseries` group,
stamping them with `iter` and using `kwargs` when opening
the JLD2 file.
"""
function jld2output!(path, iter, time, data, kwargs)
    jldopen(path, "r+"; kwargs...) do file
        file["timeseries/t/$iter"] = time
        for (name, datum) in data
            file["timeseries/$name/$iter"] = datum
        end
    end
    return
end

function start_next_file(model::Model, fw::JLD2OutputWriter)
    verbose = fw.verbose
    sz = filesize(fw.filepath)
    verbose && @info begin
        "Filesize $(pretty_filesize(sz)) has exceeded maximum file size $(pretty_filesize(fw.max_filesize))."
    end

    if fw.part == 1
        part1_path = replace(fw.filepath, r".jld2$" => "_part1.jld2")
        verbose && @info "Renaming first part: $(fw.filepath) -> $part1_path"
        mv(fw.filepath, part1_path, force=fw.force)
        fw.filepath = part1_path
    end

    fw.part += 1
    fw.filepath = replace(fw.filepath, r"part\d+.jld2$" => "part" * string(fw.part) * ".jld2")
    fw.force && isfile(fw.filepath) && rm(fw.filepath, force=true)
    verbose && @info "Now writing to: $(fw.filepath)"

    jldopen(fw.filepath, "a+"; fw.jld2_kw...) do file
        fw.init(file, model)
        saveproperties!(file, model, fw.including)
    end
end

"""
    FieldOutput([return_type=Array], field)

Returns a `FieldOutput` type intended for use with the `JLD2OutputWriter`.
Calling `FieldOutput(model)` returns `return_type(field.data.parent)`.
"""
struct FieldOutput{O, F}
    return_type :: O
          field :: F
end

FieldOutput(field) = FieldOutput(Array, field) # default
(fo::FieldOutput)(model) = fo.return_type(fo.field.data.parent)

"""
    FieldOutputs(fields)

Returns a dictionary of `FieldOutput` objects with key, value
pairs corresponding to each name and value in the `NamedTuple` `fields`.
Intended for use with `JLD2OutputWriter`.

Examples
========

```julia
julia> output_writer = JLD2OutputWriter(model, FieldOutputs(model.velocities), frequency=1);

julia> keys(output_writer.outputs)
Base.KeySet for a Dict{Symbol,FieldOutput{UnionAll,F} where F} with 3 entries. Keys:
  :w
  :v
  :u
```
"""
function FieldOutputs(fields)
    names = propertynames(fields)
    nfields = length(fields)
    return Dict((names[i], FieldOutput(fields[i])) for i in 1:nfields)
end

####
#### NetCDF output writer
####

"A type for writing NetCDF output."
Base.@kwdef mutable struct NetCDFOutputWriter <: AbstractOutputWriter
             dir  :: AbstractString = "."
          prefix  :: AbstractString = ""
        frequency :: Int    = 1
          padding :: Int    = 9
    naming_scheme :: Symbol = :iteration
      compression :: Int    = 3
            async :: Bool   = false
end

ext(fw::NetCDFOutputWriter) = ".nc"

function filename(fw, name, iteration)
    if fw.naming_scheme == :iteration
        fw.prefix * name * lpad(iteration, fw.padding, "0") * ext(fw)
    elseif fw.naming_scheme == :file_number
        file_num = Int(iteration / fw.frequency)
        fw.prefix * name * lpad(file_num, fw.padding, "0") * ext(fw)
    else
        throw(ArgumentError("Invalid naming scheme: $(fw.naming_scheme)"))
    end
end

function write_output(model::Model, fw::NetCDFOutputWriter)
    fields = Dict(
        "xC" => collect(model.grid.xC),
        "yC" => collect(model.grid.yC),
        "zC" => collect(model.grid.zC),
        "xF" => collect(model.grid.xF),
        "yF" => collect(model.grid.yF),
        "zF" => collect(model.grid.zF),
        "u" => Array(parentdata(model.velocities.u)),
        "v" => Array(parentdata(model.velocities.v)),
        "w" => Array(parentdata(model.velocities.w)),
        "T" => Array(parentdata(model.tracers.T)),
        "S" => Array(parentdata(model.tracers.S))
    )

    if fw.async
        # Execute asynchronously on worker 2.
        i = model.clock.iteration
        @async remotecall(write_output_netcdf, 2, fw, fields, i)
    else
        write_output_netcdf(fw, fields, model.clock.iteration)
    end

    return
end

function write_output_netcdf(fw::NetCDFOutputWriter, fields, iteration)
    xC, yC, zC = fields["xC"], fields["yC"], fields["zC"]
    xF, yF, zF = fields["xF"], fields["yF"], fields["zF"]

    u, v, w = fields["u"], fields["v"], fields["w"]
    T, S    = fields["T"], fields["S"]

    xC_attr = Dict("longname" => "Locations of the cell centers in the x-direction.", "units" => "m")
    yC_attr = Dict("longname" => "Locations of the cell centers in the y-direction.", "units" => "m")
    zC_attr = Dict("longname" => "Locations of the cell centers in the z-direction.", "units" => "m")

    xF_attr = Dict("longname" => "Locations of the cell faces in the x-direction.", "units" => "m")
    yF_attr = Dict("longname" => "Locations of the cell faces in the y-direction.", "units" => "m")
    zF_attr = Dict("longname" => "Locations of the cell faces in the z-direction.", "units" => "m")

    u_attr = Dict("longname" => "Velocity in the x-direction", "units" => "m/s")
    v_attr = Dict("longname" => "Velocity in the y-direction", "units" => "m/s")
    w_attr = Dict("longname" => "Velocity in the z-direction", "units" => "m/s")
    T_attr = Dict("longname" => "Temperature", "units" => "K")
    S_attr = Dict("longname" => "Salinity", "units" => "g/kg")

    filepath = joinpath(fw.dir, filename(fw, "", iteration))

    if fw.async
        println("[Worker $(Distributed.myid()): NetCDFOutputWriter] Writing fields to disk: $filepath")
    end

    isfile(filepath) && rm(filepath)

    nccreate(filepath, "u", "xF", xC, xC_attr,
                            "yC", yC, yC_attr,
                            "zC", zC, zC_attr,
                            atts=u_attr, compress=fw.compression)

    nccreate(filepath, "v", "xC", xC, xC_attr,
                            "yF", yC, yC_attr,
                            "zC", zC, zC_attr,
                            atts=v_attr, compress=fw.compression)

    nccreate(filepath, "w", "xC", xC, xC_attr,
                            "yC", yC, yC_attr,
                            "zF", zC, zC_attr,
                            atts=w_attr, compress=fw.compression)

    nccreate(filepath, "T", "xC", xC, xC_attr,
                            "yC", yC, yC_attr,
                            "zC", zC, zC_attr,
                            atts=T_attr, compress=fw.compression)

    nccreate(filepath, "S", "xC", xC, xC_attr,
                            "yC", yC, yC_attr,
                            "zC", zC, zC_attr,
                            atts=S_attr, compress=fw.compression)

    ncwrite(u, filepath, "u")
    ncwrite(v, filepath, "v")
    ncwrite(w, filepath, "w")
    ncwrite(T, filepath, "T")
    ncwrite(S, filepath, "S")

    ncclose(filepath)

    return
end

function read_output(fw::NetCDFOutputWriter, field_name, iter)
    filepath = joinpath(fw.dir, filename(fw, "", iter))
    field_data = ncread(filepath, field_name)
    ncclose(filepath)
    return field_data
end

####
#### Checkpointer
####

"""
    Checkpointer{I, T, P, A} <: AbstractOutputWriter

An output writer for checkpointing models to a JLD2 file from which models can be restored.
"""
mutable struct Checkpointer{I, T, P, A} <: AbstractOutputWriter
         frequency :: I
          interval :: T
          previous :: Float64
               dir :: String
            prefix :: String
        properties :: P
    has_array_refs :: A
             force :: Bool
           verbose :: Bool
end

"""
    Checkpointer(model; frequency=nothing, interval=nothing, dir=".", prefix="checkpoint",
                        force=false, verbose=false,
                        properties = [:architecture, :boundary_conditions, :grid, :clock,
                                      :eos, :constants, :closure, :velocities, :tracers,
                                      :timestepper])

Construct a `Checkpointer` that checkpoints the model to a JLD2 file every so often as
specified by `frequency` or `interval`. The `model.clock.iteration` is included in the
filename to distinguish between multiple checkpoint files.

Note that extra model `properties` can be safely specified, but removing crucial properties
such as `:velocities` will make restoring from the checkpoint impossible.

The checkpoint file is generated by serializing model properties to JLD2. However,
functions cannot be serialized to disk (at least not with JLD2). So if a model property
contains a reference somewhere in its hierarchy it will not be included in the checkpoint
file (and you will have to manually restore them).

Keyword arguments
=================

    - `frequency::Int`   : Save output every `n` model iterations.
    - `interval::Int`    : Save output every `t` units of model clock time.
    - `dir::String`      : Directory to save output to. Default: "." (current working
                           directory).
    - `prefix::String`   : Descriptive filename prefixed to all output files.
                           Default: "checkpoint".
    - `force::Bool`      : Remove existing files if their filenames conflict.
                           Default: `false`.
    - `verbose::Bool`    : Log what the output writer is doing with statistics on
                           compute/write times and file sizes. Default: `false`.
    - `properties::Array`: List of model properties to checkpoint. By default,
                           properties = [:architecture, :boundary_conditions, :grid,
                                         :clock, :eos, :constants, :closure, :velocities,
                                         :tracers, :timestepper]
"""
function Checkpointer(model; frequency=nothing, interval=nothing, dir=".", prefix="checkpoint", force=false,
                      verbose=false, properties = [:architecture, :boundary_conditions, :grid, :clock, :coriolis,
                                                   :buoyancy, :closure, :velocities, :tracers, :timestepper])

    validate_interval(frequency, interval)

    # Grid needs to be checkpointed for restoring to work.
    :grid ∉ properties && push!(properties, :grid)

    has_array_refs = Symbol[]

    for p in properties
        isa(p, Symbol) || @error "Property $p to be checkpointed must be a Symbol."
        p ∉ propertynames(model) && @error "Cannot checkpoint $p, it is not a model property!"

        if has_reference(Function, getproperty(model, p))
            @warn "model.$p contains a function somewhere in its hierarchy and will not be checkpointed."
            filter!(e -> e != p, properties)
        end

        has_reference(AbstractField, getproperty(model, p)) && push!(has_array_refs, p)
    end

    mkpath(dir)
    return Checkpointer(frequency, interval, 0.0, dir, prefix, properties, has_array_refs, force, verbose)
end

function write_output(model, c::Checkpointer)
    filepath = joinpath(c.dir, c.prefix * string(model.clock.iteration) * ".jld2")
    c.verbose && @info "Checkpointing to file $filepath..."

    t0 = time_ns()
    jldopen(filepath, "w") do file
        file["checkpointed_properties"] = c.properties
        file["has_array_refs"] = c.has_array_refs

        serializeproperties!(file, model, filter(e -> !(e ∈ c.has_array_refs), c.properties))
        saveproperties!(file, model, filter(e -> e ∈ c.has_array_refs, c.properties))
    end

    t1, sz = time_ns(), filesize(filepath)
    c.verbose && @info "Checkpointing done: time=$(prettytime((t1-t0)/1e9)), size=$(pretty_filesize(sz))"
end

defaultname(::Checkpointer, nelems) = :checkpointer
_arr(::CPU, a) = a
_arr(::GPU, a) = CuArray(a)

function restore_fields!(model, file, arch, fieldset; location="$fieldset")
    if fieldset == :timestepper
        restore_fields!(model.timestepper, file, arch, :Gⁿ; location="timestepper/Gⁿ")
        restore_fields!(model.timestepper, file, arch, :G⁻; location="timestepper/G⁻")
    else
        for p in propertynames(getproperty(model, fieldset))
            getproperty(getproperty(model, fieldset), p).data.parent .= _arr(arch, file[location * "/$p"])
        end
    end
end

"""
    restore_from_checkpoint(filepath; kwargs=Dict())

Restore a model from the checkpoint file stored at `filepath`. `kwargs` can be passed
to the `Model` constructor, which can be especially useful if you need to manually
restore forcing functions or boundary conditions that rely on functions.
"""
function restore_from_checkpoint(filepath; kwargs=Dict())
    file = jldopen(filepath, "r")

    cps = file["checkpointed_properties"]
    has_array_refs = file["has_array_refs"]

    # Restore model properties that were just serialized.
    # We skip fields that contain structs and restore them after model creation.
    for p in cps
        if p ∉ has_array_refs
            kwargs[p] = file["$p"]
        end
    end

    model = Model(; kwargs...)

    # Now restore fields.
    for p in cps
        if p in has_array_refs
            restore_fields!(model, file, model.architecture, p)
        end
    end

    close(file)

    return model
end


text:
using RunTests
using Base.Test

@testmodule DataTransformTest begin
using MLLibrary

  function test_proxyColumns()
    gender = ["M" "F" "M" "M" "F"]
    month=[1 3 2 1 8]
    A=[gender' month']
    proxycolumns = [Dict("column"=>1,"proxy"=>["M" "F"]) Dict("column"=>2,"proxy"=>collect(1:12))]
    @test MLLibrary.DataManipulation.DataTransform.addProxyColumns(A,proxycolumns)==
    [1 0 1 0 0 0 0 0 0 0 0 0 0 0;
     0 1 0 0 1 0 0 0 0 0 0 0 0 0;
     1 0 0 1 0 0 0 0 0 0 0 0 0 0;
     1 0 1 0 0 0 0 0 0 0 0 0 0 0;
     0 1 0 0 0 0 0 0 0 1 0 0 0 0]
    # @test 1==0
  end

  function test_zScale()
    A=[1 2 3 4;5 6 7 8;
       9 10 11 12; 13 14 15 16]
    @test MLLibrary.DataManipulation.DataTransform.zScale(MLLibrary.DataManipulation.DataTransform.convertToFloat(A))==[-1.161895003862225 -1.161895003862225 -1.161895003862225 -1.161895003862225
 -0.3872983346207417 -0.3872983346207417 -0.3872983346207417 -0.3872983346207417
 0.3872983346207417 0.3872983346207417 0.3872983346207417 0.3872983346207417
 1.161895003862225 1.161895003862225 1.161895003862225 1.161895003862225]

 @test MLLibrary.DataManipulation.DataTransform.zScale(A, 1, 1)==[0 1 2 3; 4 5 6 7; 8 9 10 11; 12 13 14 15]
  end

  function test_floatConversion()
    A=["1" 2 3.1 "3.2"]
    @test round(MLLibrary.DataManipulation.DataTransform.convertToFloat(A),1)==[1. 2. 3.1 3.2]
  end

  function test_zScaleColumn()
    A=MLLibrary.DataManipulation.DataTransform.convertToFloat([1 2 3 4;5 6 7 8;
       9 10 11 12; 13 14 15 16])
    Expected = [-1.161895003862225;
        -0.3872983346207417;
         0.3872983346207417;
         1.161895003862225]
    @test MLLibrary.DataManipulation.DataTransform.zScaleColumn(A,1)==[Expected A[:,2:end]]
         @test MLLibrary.DataManipulation.DataTransform.zScaleColumn(A,2)==[Expected Expected A[:,3:end]]
         @test MLLibrary.DataManipulation.DataTransform.zScaleColumn(A,3)==[Expected Expected Expected A[:,4]]
         @test MLLibrary.DataManipulation.DataTransform.zScaleColumn(A,4)==[Expected Expected Expected Expected]

  end

end

ext: jl

max_stars_repo_name: ajermaky/CustomMLLibrary

max_issues_repo_name: ajermaky/CustomMLLibrary

max_forks_repo_name: ajermaky/CustomMLLibrary


text:
@testset "693.binary-number-with-alternating-bits.jl" begin
    @test !any(has_alternating_bits.(6:9))
    @test has_alternating_bits(5)
    @test has_alternating_bits(10)
end


text:
module cheta

import PyCall: PyObject, pyimport, PyNULL, pycall

fetch = PyNULL()

function __init__()
    copy!(fetch, pyimport("cheta.fetch_sci"))
end

# MSID

struct MSID
    msid::PyObject
    function MSID(msid::PyObject)
        new(msid)
    end
end

function MSID(msid::String, start, stop)
    m = fetch.MSID(msid, start, stop)
    MSID(m)
end

function Base.getproperty(m::MSID, v::Symbol)
    if v == :msid
        return getfield(m, :msid)
    else
        return getproperty(m.msid, v)
    end
end

Base.show(io::IO, msid::MSID) = print(io, msid.msid."__repr__"())

function interpolate(msid::MSID; dt=nothing, start=nothing, stop=nothing, times=nothing)
    msid.msid.interpolate(dt, start, stop, times)
    nothing
end

# MSIDset

struct MSIDset
    msidset::PyObject
    function MSIDset(msidset::PyObject)
        new(msidset)
    end
end

function MSIDset(msids::Array{String}, start, stop)
    m = pycall(fetch.MSIDset, PyObject, msids, start, stop)
    MSIDset(m)
end

function Base.getindex(m::MSIDset, i::String)
    MSID(m.msidset."__getitem__"(i))
end

Base.show(io::IO, m::MSIDset) = print(io, m.msidset."__repr__"())

end


text:
module WolframLinearAlgebra

import Base: OneTo, HasShape, promote_rule, convert, zero, one, size, axes, ndims, length
import Base: firstindex, lastindex, IteratorSize, keys, first, last, copy
import Base: +, -, *, /, ^, inv, sqrt, conj, adjoint, transpose
import Base: sin, cos, tan, sinh, cosh, tanh, asin, acos, atan, asinh, acosh, atanh
import Base: log, log2, log10, exp
import Base: string, show

import LinearAlgebra: Eigen, eigen, SVD, svd, Schur, schur
import LinearAlgebra: qr, lq
import LinearAlgebra: pinv, nullspace, diag

import MathLink: WSymbol, WExpr, WInteger, WReal, weval, parseexpr, @W_str, @W_cmd

export WElem, WSolver, @WE_str, @WE_cmd, @WS_str, @W_str, @W_cmd
export weval, simplify, n
export eigen, svd, qr, lq, schur, pinv, nullspace

iswlist(x) = false
iswlist(x::WExpr) = x.head == W"List"

include("welem.jl")
include("wsolver.jl")
include("linargfunctions.jl")

end # module


text:
module ElasticGD

export elastic_desc
export println1
export get_mse
export get_sens
export get_spec
export get_sens_spec
export get_path_frac
export train_val_test
export sweep_alpha
export make_folds

using Random
using Statistics
using LinearAlgebra
using Lasso
using StatsBase
using Combinatorics

function elastic_desc(X, y, alpha, step_size, alg="norm_ns", break_early=false)
  beta = zeros(size(X)[2],1)
  grad = nothing
  betas = [beta]
  grads = Array{Float64,2}[]
  old_norm_grad = Inf
  XtXn=X'*X/size(X)[1]
  Xtyn=X'*y/size(X)[1]
  c_alpha=(alpha+sqrt(alpha^2+4)-2)/(2*alpha)
  for i in 1:round(Int,100/step_size)
    grad=XtXn*beta-Xtyn
    I01=(abs.(grad)/maximum(abs.(grad)).>=alpha)
    elastic_grad=I01.*grad

    delta_beta=Nothing
    if alg=="unnorm"
      delta_beta=(alpha*sign.(elastic_grad)+(1-alpha)*elastic_grad)*step_size
    elseif alg[1:4]=="norm"
      delta_beta=(alpha*elastic_grad/norm(elastic_grad,1)+(1-alpha)*elastic_grad/norm(elastic_grad,2))*step_size
      if alg=="norm_s"
        q1=(norm(elastic_grad,1)/norm(elastic_grad,2))^2
        c=(sqrt(q1*(alpha^2*q1 + 4*(1-alpha)))-alpha*q1)/(2*(1-alpha)*(sqrt(q1)*(1-alpha)+alpha))
        delta_beta*=c
      end
    elseif alg[1:2]=="gs"
      c=1
      if alg=="gs_s"
        q1=(norm(elastic_grad,1)/norm(elastic_grad,2))^2
        c=((sqrt(sqrt(q1)*2*alpha*sqrt(alpha^2*q1 + 4*step_size*(1-alpha)) + q1*((1-alpha)^3-2*alpha^2)) - sqrt(q1)*(1 - alpha)*sqrt(1-alpha))/(2*sqrt(step_size)*alpha*sqrt(1- alpha)))^2
      end
      delta_beta=c*step_size*alpha*elastic_grad/norm(elastic_grad,1)+np.sqrt(c*step_size)*(1-alpha)*elastic_grad/norm(elastic_grad,2)
    end
    beta-=delta_beta

    push!(betas,beta)
    push!(grads,grad)

    if break_early
      norm_grad=norm(grad)
      if norm_grad<0.1*step_size
        break
      end
      if mod(i,round(Int, 1/step_size))==0
        if round(norm_grad/old_norm_grad; digits=4)==1
          break
        end
        old_norm_grad = norm_grad
      end
    end
  end
  return hcat(betas...)', hcat(grads...)'
end

function println1(LOG_NAME, print_str)
  if "SLURM_SUBMIT_DIR" in keys(ENV)
    f = open(string(ENV["SLURM_SUBMIT_DIR"],"/logs/",LOG_NAME,".txt"), "a")
    println(f, print_str)
    close(f)
  else
    println(print_str)
    #f = open("temp.txt", "a")
    #println(f, print_str)
    #close(f)
  end
end

function get_mse(y,y_hat)
  return mean((y.-y_hat).^2)
end

function get_sens(beta, beta_star)
  tp = sum((beta.!=0) .& (beta_star.!=0))
  p = sum(beta_star.!=0)
  return tp/p
end

function get_spec(beta, beta_star)
  tn = sum((beta.==0) .& (beta_star.==0))
  n = sum(beta_star.==0)
  return tn/n
end

function get_sens_spec(beta_eg, beta_en, beta_star)
  return get_sens(beta_eg, beta_star), get_sens(beta_en, beta_star), get_spec(beta_eg, beta_star), get_spec(beta_en, beta_star)
end

function get_path_frac(beta_path, beta_star)
  norms=sum(abs.(beta_path),dims=2)
  correct_betas=Bool[]
  for norm in 0:0.1:norms[end]
    idx=findlast(norm.>=norms)[1]
    beta_bool=beta_path[idx,:].!=0
    #if !all(beta_bool) #TODO: 0?
    push!(correct_betas, sum(.!((beta_bool.==0) .⊻ (beta_star.==0)))==length(beta_star))
    #end
  end
  return mean(correct_betas)
end

function get_idxs(n, splits)
  @assert sum(splits)==1
  r1 = round(Int, splits[1]*n)
  r2 = r1+round(Int, splits[2]*n)
  train_idxs=1:r1
  val_idxs=r1+1:r2
  if length(splits)==3
    test_idxs=r2+1:n
  else
    test_idxs=1:0
  end
  return train_idxs, val_idxs, test_idxs
end

function make_folds(n,n_folds; test=true, val_frac=0.15)
  all_folds=1:n_folds
  n_val_test=max(1,round(Int,val_frac*n_folds))
  fold_list=Array{Array{Int64,1},1}[]
  if test
    for test_folds in combinations(all_folds,n_val_test)
      train_val_folds=setdiff(all_folds, test_folds)
      for val_folds in combinations(train_val_folds,n_val_test)
        train_folds=setdiff(train_val_folds, val_folds)
        push!(fold_list, [train_folds, val_folds, test_folds])
      end
    end
  else
    for val_folds in combinations(all_folds,n_val_test)
      train_folds=setdiff(all_folds, val_folds)
      push!(fold_list, [train_folds, val_folds])
    end
  end
  n_idx_fold=n÷n_folds
  n_folds_extra=mod(n,n_folds)
  fold_ranges=UnitRange{Int64}[]
  for fold in 1:n_folds
    if fold<=n_folds_extra
      fold_range=(fold-1)*(n_idx_fold+1)+1:fold*(n_idx_fold+1)
    else
      fold_range=n_folds_extra+(fold-1)*n_idx_fold+1:n_folds_extra+fold*n_idx_fold
    end
    push!(fold_ranges,fold_range)
  end
  return fold_list, fold_ranges
end

function normalize_X(X)
  return (X.-mean(X, dims=1))./std(X, dims=1)
end

function normalize_y(y)
  return (y.-mean(y))
end

function train_val_test(X,y; splits=nothing, fold_list=nothing, fold_ranges=nothing, comb_idx=nothing, normalize=true, seed=0)
  Random.seed!(seed)
  p=randperm(size(X)[1])
  Xp=X[p,:]
  yp=y[p]
  if fold_list==nothing
    train_idxs, val_idxs, test_idxs = get_idxs(size(X)[1], splits)
  else
    flc=fold_list[comb_idx]
    if length(flc)==2
      push!(flc,[])
    end
    train_idxs, val_idxs, test_idxs = map(x->collect(Iterators.flatten(fold_ranges[x])), flc)
  end
  X_train=Xp[train_idxs,:]
  X_val=Xp[val_idxs,:]
  X_test=Xp[test_idxs,:]
  y_train=yp[train_idxs]
  y_val=yp[val_idxs]
  y_test=yp[test_idxs]
  if normalize
    X_train = normalize_X(X_train)
    X_val = normalize_X(X_val)
    X_test = normalize_X(X_test)
    y_train = normalize_y(y_train)
    y_val = normalize_y(y_val)
    y_test = normalize_y(y_test)
  end
  return X_train, y_train, X_val, y_val, X_test, y_test
end

function most_common_beta(beta_path, incl_ridge)
  norms=sum(abs.(beta_path),dims=2)
  beta_bools=Array{Bool,1}[]
  for norm in 0:0.1:norms[end]
    idx=findlast(norm.>=norms)[1]
    push!(beta_bools,beta_path[idx,:].!=0)
  end
  beta_counts=countmap(beta_bools)
  best_count=0
  best_beta=ones(size(beta_path)[2])
  for (beta, count) in beta_counts
    if count>best_count && any(beta) && (!all(beta) || incl_ridge)
#      more common        not only 0s      only 1s     include only ones
      best_count=count
      best_beta=beta
    end
  end
  return best_beta
end

function eval_alpha(beta_path, X_val, y_val)
  best_mse=Inf
  best_r=nothing
  for r in 1:size(beta_path)[1]
    beta=beta_path[r,:]
    mse=get_mse(y_val, X_val*beta)
    if mse<best_mse
      best_r=r
      best_mse=mse
    end
  end
  return best_r, best_mse
end

function select_alpha(mses, best_rs, beta_paths; alphas, mse_frac=0)
  alpha_idx=findlast(mses/minimum(mses).<=1+mse_frac)
  best_beta_path=beta_paths[alpha_idx]
  best_r=best_rs[alpha_idx]
  best_beta=best_beta_path[best_r,:]
  return best_beta, best_beta_path, alphas[alpha_idx]
end

function sweep_alpha(X_train, y_train, X_val, y_val, alg; mse_frac=0, n_folds=10, alphas=0.01:0.01:1, step_size=1e-2, nλ=100, λminratio=ifelse(size(X_train, 1) < size(X_train, 2), 0.01, 1e-4))
  best_rs_eg= Int64[]
  best_rs_en= Int64[]
  beta_paths_eg= Array{Float64,2}[]
  beta_paths_en= Array{Float64,2}[]
  mses_eg= Float64[]
  mses_en= Float64[]
  for alpha in alphas
    beta_path_eg = elastic_desc(X_train, y_train, alpha, step_size, alg)[1]
    best_r_eg, best_mse_eg = eval_alpha(beta_path_eg, X_val, y_val)
    push!(beta_paths_eg, beta_path_eg)
    push!(best_rs_eg, best_r_eg)
    push!(mses_eg, best_mse_eg)

    beta_path_en=fit(LassoPath, X_train, y_train; α=alpha, nλ=nλ, λminratio=λminratio).coefs'
    best_r_en, best_mse_en = eval_alpha(beta_path_en, X_val, y_val)
    push!(beta_paths_en, beta_path_en)
    push!(best_rs_en, best_r_en)
    push!(mses_en, best_mse_en)
  end
  beta_mse_eg, beta_path_eg, alpha_eg = select_alpha(mses_eg, best_rs_eg, beta_paths_eg; alphas=alphas, mse_frac=mse_frac)
  beta_mse_en, beta_path_en, alpha_en = select_alpha(mses_en, best_rs_en, beta_paths_en; alphas=alphas, mse_frac=mse_frac)
  beta_cv_eg, beta_cv_en= select_glmnet(X_train, y_train, alpha_eg, alpha_en, alg; n_folds=n_folds, step_size=step_size, nλ=nλ, λminratio=λminratio)
  beta_mc_eg = most_common_beta(beta_path_eg, true)
  beta_mc_en = most_common_beta(beta_path_en, true)
  beta_mc1_eg = most_common_beta(beta_path_eg, false)
  beta_mc1_en = most_common_beta(beta_path_en, false)
  return beta_mse_eg, beta_mse_en, beta_cv_eg, beta_cv_en, beta_path_eg, beta_path_en, alpha_eg, alpha_en, beta_mc_eg, beta_mc_en, beta_mc1_eg, beta_mc1_en
end

function select_glmnet(X_train, y_train, alpha_eg, alpha_en, alg; n_folds=10, step_size=1e-2, nλ=100, λminratio=ifelse(size(X_train, 1) < size(X_train, 2), 0.01, 1e-4))
  fold_list, fold_ranges = make_folds(size(X_train)[1], n_folds; test=false, val_frac=1/n_folds)
  λmax=maximum(abs.(X_train'*y_train/length(y_train)))/alpha_en
  logλmax = log(λmax)
  lbdas=exp.(range(logλmax, stop=logλmax + log(λminratio), length=nλ))
  mses_eg= Array{Float64,1}[]
  mses_en= Array{Float64,1}[]
  for k in 1:n_folds
    X_train_cv, y_train_cv, X_val_cv, y_val_cv, _, _ = train_val_test(X_train, y_train; fold_list=fold_list, fold_ranges=fold_ranges, comb_idx=k, seed=4, normalize=false)

    beta_path_eg = elastic_desc(X_train_cv, y_train_cv, alpha_eg, step_size, alg)[1]
    mse_eg=map(beta->get_mse(y_val_cv, X_val_cv*beta), map(r->beta_path_eg[r,:],1:size(beta_path_eg)[1]))
    push!(mses_eg, mse_eg)

    beta_path_en=fit(LassoPath, X_train_cv, y_train_cv; α=alpha_en, λ=lbdas).coefs'
    mse_en=map(beta->get_mse(y_val_cv, X_val_cv*beta), map(r->beta_path_en[r,:],1:size(beta_path_en)[1]))
    push!(mses_en, mse_en)
  end
  r_eg=select_pen(mses_eg)

  beta_eg = elastic_desc(X_train, y_train, alpha_eg, step_size, alg)[1][r_eg,:]
  r_en=select_pen(mses_en)
  beta_en=collect(fit(LassoPath, X_train, y_train; α=alpha_en, nλ=nλ, λminratio=λminratio).coefs'[r_en,:])
  return beta_eg, beta_en
end

function select_pen(mses)
  min_len=minimum(map(length,mses))
  mses=hcat(map(mses->mses[1:min_len],mses)...)
  min_mean=minimum(mean(mses, dims=2))
  return findfirst(mean(mses, dims=2).-std(mses, dims=2)/sqrt(size(mses)[2]).<min_mean)[1]
end

function to_norm(beta_path)
  betas=Array{Float64,1}[]
  norms=sum(abs.(beta_path),dims=2)
  for norm_i in 0:1:norms[end]
    idx=findfirst(norm_i.<=norms)[1]
    push!(betas,beta_path[idx,:])
  end
  return hcat(betas...)'
end
end


text:
using Documenter, PCM

makedocs(
    modules = [PCM],
    format = :html,
    checkdocs = :exports,
    sitename = "PCM.jl",
    pages = Any["index.md"]
)

deploydocs(
    repo = "github.com/hildebrandmw/PCM.jl.git",
)


text:
# If successful, returns X such that X * A * X^-1 == B
function _issimilar_husert_generic(A, B)
  @assert issquarefree(minpoly(A))
  m = nrows(A)
  f = charpoly(A)
  fac = factor(f)
  mus = typeof(f)[]
  ns = Int[]
  ds = Int[]
  for (p, n) in fac
    push!(mus, p)
    push!(ns, n)
    push!(ds, degree(p))
  end

  s = length(mus)
  Ks = AnticNumberField[]
  vecsA = Matrix{nf_elem}(undef, m, sum(ns))
  vecsB = Matrix{nf_elem}(undef, m, sum(ns))
  k = 1
  for i in 1:s
    K, _ = NumberField(mus[i], "a", cached = false)
    EA = eigenspace(change_base_ring(K, A), gen(K))
    EB = eigenspace(change_base_ring(K, B), gen(K))
    push!(Ks, K)
    if length(EA) != length(EB)
      return false, zero_matrix(FlintQQ, 0, 0)
    end
    @assert length(EA) == length(EB)
    for j in 1:length(EA)
      for l in 1:m
        vecsA[l, k] = EA[j][l]
        vecsB[l, k] = EB[j][l]
      end
      k += 1
    end
  end

  @show Ks, ns
#  for i in 1:m
#    @show view(vecsA, i, :)
#    @show _to_relative_basis(_to_absolute_basis(view(vecsA, i, :), m, ns, Ks), m, ns, Ks)
#  end

  absolute_basis_vec_A = Vector{Vector{fmpq}}()
  absolute_basis_vec_B = Vector{Vector{fmpq}}()
  absolute_basis_A = zero_matrix(FlintQQ, m, m)
  absolute_basis_B = zero_matrix(FlintQQ, m, m)
  for i in 1:m
    vA = _to_absolute_basis(view(vecsA, i, :), m, ns, Ks)
    vB = _to_absolute_basis(view(vecsB, i, :), m, ns, Ks)
    push!(absolute_basis_vec_A, vA)
    push!(absolute_basis_vec_B, vB)
    for j in 1:m
      absolute_basis_A[i, j] = vA[j]
      absolute_basis_B[i, j] = vB[j]
    end
  end

  @show "asds"

  # Now construct the colon ideal
  # First the Q-basis of \prod Mat(n_i, K_i)
  actions = Vector{fmpq_mat}()
  another_basis_of_actions = Vector{Vector{dense_matrix_type(nf_elem)}}()

  for i in 1:s
    ni = ns[i]
    for i1 in 1:ni
      for i2 in 1:ni
        for j in 1:degree(Ks[i])
          V = Vector{dense_matrix_type(nf_elem)}(undef, s)
          M = representation_matrix(gen(Ks[i])^(j - 1))
          z = zero_matrix(Ks[i], ni, ni)
          z[i1, i2] = gen(Ks[i])^(j-1)
          V[i] = z
          for j in 1:s
            if j != i
              V[j] = zero_matrix(Ks[j], ns[j], ns[j])
            end
          end
          push!(another_basis_of_actions, V)
          # I need to construct the block matrix, where M is at block (i1, i2)
          # and all the other ones are zero
          # The overall matrix must have size ni * degree(Ks[i])
          Z1 = zero_matrix(FlintQQ, (i1 - 1) * ds[i], ni * ds[i])
          Z2 = zero_matrix(FlintQQ, ds[i], (i2 - 1) * ds[i])
          Z3 = hcat(Z2, M)
          Z3 = hcat(Z3, zero_matrix(FlintQQ, ds[i], (ni - i2) * ds[i]))
          Z4 = zero_matrix(FlintQQ, (ni - i1) * ds[i], ni * ds[i])
          MM = vcat(Z1, Z3, Z4)
          blocklengthbefore = sum(Int[ns[o] * ds[o] for o in 1:(i - 1)])
          blocklengthafter = sum(Int[ns[o] * ds[o] for o in (i + 1):s])
          D = diagonal_matrix(zero_matrix(FlintQQ, blocklengthbefore, blocklengthbefore), MM, zero_matrix(FlintQQ, blocklengthafter, blocklengthafter))
          push!(actions, D)
        end
      end
    end
  end

  absolute_basis_B_inv = inv(absolute_basis_B)
  absolute_basis_A_inv = inv(absolute_basis_A)

  # So D is a Q-basis of Mat_n(K)

  # Now determine the colon ideal

  la = length(actions)

  M = zero_matrix(FlintQQ, la, 0)

  for a in absolute_basis_vec_A
    _M = zero_matrix(FlintQQ, la, m)
    for i in 1:length(actions)
      D = actions[i]
      v = matrix(FlintQQ, 1, m, a) * D * absolute_basis_B_inv
      for j in 1:m
        _M[i, j] = v[1, j]
      end
    end
    M = hcat(M, _M)
  end

  M = sub(hnf(FakeFmpqMat(transpose(M)), :upperright), 1:la, 1:la)
  N = inv(M)

  SS = N

  M = zero_matrix(FlintQQ, length(actions), 0)

  for a in absolute_basis_vec_B
    _M = zero_matrix(FlintQQ, length(actions), m)
    for i in 1:length(actions)
      D = actions[i]
      v = matrix(FlintQQ, 1, m, a) * D * absolute_basis_B_inv
      for j in 1:m
        _M[i, j] = v[j]
      end
    end
    M = hcat(M, _M)
  end
  M = sub(hnf(FakeFmpqMat(transpose(M)), :upperright), 1:la, 1:la)
  N = inv(M)

  bcolonb = N

  basis_of_colon_ideal = fmpq_mat[sum(SS[i, j] * actions[i] for i in 1:length(actions)) for j in 1:nrows(SS)]

  basis_of_colon_idealAA = fmpq_mat[sum(bcolonb[i, j] * actions[i] for i in 1:length(actions)) for j in 1:nrows(bcolonb)]

  other_basis_of_colon_ideal = Vector{dense_matrix_type(nf_elem)}[]
  other_basis_of_colon_idealAA = Vector{dense_matrix_type(nf_elem)}[]

  for j in 1:nrows(SS)
    push!(other_basis_of_colon_ideal, foldl((x, y) -> x .+ y, [SS[i, j] .* another_basis_of_actions[i] for i in 1:nrows(SS)]))
  end

  for j in 1:nrows(bcolonb)
    push!(other_basis_of_colon_idealAA, foldl((x, y) -> x .+ y, [bcolonb[i, j] .* another_basis_of_actions[i] for i in 1:nrows(bcolonb)]))
  end

  for a in absolute_basis_vec_A
    for i in 1:length(basis_of_colon_ideal)
      @assert denominator(matrix(FlintQQ, 1, m, a)  * basis_of_colon_ideal[i] * absolute_basis_B_inv) == 1
      #@assert denominator(matrix(FlintQQ, 1, m, a)  * basis_of_colon_idealAA[i] * absolute_basis_A_inv) == 1
    end
  end

  A = _create_algebra_husert(Ks, ns)

  idealAgens = elem_type(A)[]

  ordergens = elem_type(A)[]

  dec = decompose(A)

  for b in other_basis_of_colon_ideal
    z = zero(A)
    @assert length(dec) == length(b)
    for i in 1:length(dec)
      B, mB = dec[i]
      C, BtoC = B.isomorphic_full_matrix_algebra
      z = z + mB(preimage(BtoC, C(b[i]))::elem_type(B))
    end
    push!(idealAgens, z)
  end

  for b in other_basis_of_colon_idealAA
    z = zero(A)
    @assert length(dec) == length(b)
    for i in 1:length(dec)
      B, mB = dec[i]
      local C::AlgMat{nf_elem, Generic.MatSpaceElem{nf_elem}}
      C, BtoC = B.isomorphic_full_matrix_algebra
      z = z + mB(preimage(BtoC, C(b[i]))::elem_type(B))
    end
    push!(ordergens, z)
  end

  ide = ideal_from_lattice_gens(A, idealAgens)
  ideO = Order(A, ordergens)

  @assert ideO == right_order(ide)

  #@assert right_order(ide) == ideO

  fl, y = _isprincipal(ide, ideO, :right)

  if fl
    dec = decompose(A)
    D = Vector{Generic.MatSpaceElem{nf_elem}}(undef, length(dec))
    for i in 1:length(dec)
      B, mB = dec[i]
      local C::AlgMat{nf_elem, Generic.MatSpaceElem{nf_elem}}
      C, BtoC = B.isomorphic_full_matrix_algebra
      z = BtoC(mB\y)::elem_type(C)
      D[i] = matrix(z)
    end
    DD = diagonal_matrix(fmpq_mat[_explode(D[i]) for i in 1:length(D)])
    return fl, transpose(absolute_basis_A * DD * absolute_basis_B_inv)
  else
    return false, zero_matrix(FlintQQ, 0, 0)
  end
end

function _explode(x::Generic.MatSpaceElem{nf_elem})
  K = base_ring(x)
  d = degree(K)
  n = nrows(x)
  m = ncols(x)
  z = zero_matrix(FlintQQ, 0, n * d)
  for i in 1:n
    zz = zero_matrix(FlintQQ, d, 0)
    for j in 1:m
      M = representation_matrix(x[i, j])
      zz = hcat(zz, M)
    end
    z = vcat(z, zz)
  end
  return z
end

Base.:(*)(x::fmpq, y::AbstractAlgebra.Generic.MatSpaceElem{nf_elem}) = base_ring(y)(x) * y

function _to_absolute_basis(v, m, ns, Ks)
  w = Vector{fmpq}(undef, m)
  k = 1
  for i in 1:sum(ns)
    vi = v[i]
    #@show i, vi
    #@assert parent(vi) == Ks[i]
    K = parent(vi)
    for j in 1:degree(K)
      w[k] = coeff(vi, j - 1)
      k += 1
    end
  end
  return w
end

function _to_relative_basis(v, m, ns, Ks)
  w = Vector{nf_elem}(undef, sum(ns))
  k = 1
  l = 1
  for i in 1:length(ns)
    a = gen(Ks[i])
    for j in 1:ns[i]
      w[k] = zero(Ks[i])
      for n in 1:degree(Ks[i])
        w[k] = w[k] + a^(n - 1) * v[l]
        l += 1
      end
      k += 1
    end
  end
  return w
end

function _create_algebra_husert(Ks, ns)
  s = length(Ks)
  algs = AlgAss{fmpq}[]
  for i in 1:s
    A = matrix_algebra(Ks[i], ns[i])
    B, BtoA = AlgAss(A)
    C, CtoB = restrict_scalars(B, FlintQQ)
    C.isomorphic_full_matrix_algebra = (A, CtoB * BtoA)
    push!(algs, C)
  end
  return _my_direct_product(algs)
end

function _issimilar_husert(A::fmpz_mat, B::fmpz_mat)
  QA = change_base_ring(FlintQQ, A)
  QB = change_base_ring(FlintQQ, B)
  fl, QC = _issimilar_husert_generic(QA, QB)
  return fl, change_base_ring(FlintZZ, QC)
end

################################################################################
#
#  General approach
#
################################################################################

mutable struct CommutatorAlgebra
  A# matrix
  K# number fields of the minpoly
  Jblocks # Jordanblocks
  Eig# Eigenbases over K
  EigAbs# absolute bases of Eigenspaces

  CommutatorAlgebra() = new()
end

function _issimilar_new(A::fmpz_mat, B::fmpz_mat)
  AQ = change_base_ring(FlintQQ, A)
  BQ = change_base_ring(FlintQQ, B)
  return _issimilar_new(AQ, BQ)
end

global _debug = []

function _issimilar_new(A, B)
  Z = _create_com_alg(A)
  ns = Int[length(E) for E in Z.Eig]
  AA = _create_algebra_husert(Z.K, ns)
  push!(_debug, (Z, AA))
  O = _basis_of_integral_commutator_algebra(A, A)
  I = _basis_of_integral_commutator_algebra(A, B)
  ordergens = elem_type(AA)[]
  idealgens = elem_type(AA)[]
  dec = decompose(AA)

  CA, TA = rational_canonical_form(A)
  CB, TB = rational_canonical_form(B)

  if CA != CB
    return false
  end

  @show Z.K
  @show ns

  @show O

  for bb in O
    b = _induce_action(bb, Z)
    z = zero(AA)
    @assert length(dec) == length(b)
    for i in 1:length(dec)
      B, mB = dec[i]
      local C::AlgMat{nf_elem, Generic.MatSpaceElem{nf_elem}}
      C, BtoC = B.isomorphic_full_matrix_algebra
      @show b[i]
      z = z + mB(preimage(BtoC, C(b[i]))::elem_type(B))
    end
    @show z
    push!(ordergens, z)
  end

  for bb in I
    b = _induce_action(inv(TA) * TB * bb, Z)
    z = zero(AA)
    @assert length(dec) == length(b)
    for i in 1:length(dec)
      B, mB = dec[i]
      local C::AlgMat{nf_elem, Generic.MatSpaceElem{nf_elem}}
      C, BtoC = B.isomorphic_full_matrix_algebra
      z = z + mB(preimage(BtoC, C(b[i]))::elem_type(B))
    end
    push!(idealgens, z)
  end

  @show dim(AA)

  OO = Order(AA, ordergens)
  OI = ideal_from_lattice_gens(AA, idealgens)
  @assert OO == right_order(OI)
  fl, y = _isprincipal(OI, OO, :right)
  return fl, OI
end

function _create_com_alg(A)
  f = factor(charpoly(A))
  Ks = []
  z = CommutatorAlgebra()
  z.A = A
  z.Jblocks = []
  z.Eig = []
  z.EigAbs = []
  z.K = Ks
  for (p, e) in f
    K, a = NumberField(p, cached = false)
    push!(Ks, K)
    res = _extract_bases_of_jordan_blocks(A, K)
    push!(z.Jblocks, res)
    zz = []
    zzz = []
    for i in 1:length(res)
      v = matrix(K, 1, length(res[i][end]), res[i][end])
      push!(zz, v)
      for j in 0:(degree(K) - 1)
        push!(zzz, a^j * v)
      end
    end
    push!(z.Eig, zz)
    push!(z.EigAbs, zzz)
  end
  return z
end

function _basis_of_commutator_algebra(A)
  linind = transpose(LinearIndices(size(A)))
  cartind = cartesian_product_iterator([1:x for x in size(A)])
  n = nrows(A)
  z = zero_matrix(base_ring(A), n^2, n^2)
  for i in 1:n
    for j in 1:n
      for k in 1:n
        z[linind[i, j], linind[i, k]] = z[linind[i, j], linind[i, k]] + A[k, j]
        z[linind[i, j], linind[k, j]] = z[linind[i, j], linind[k, j]] - A[i, k]
      end
    end
  end
  r, K = right_kernel(z)
  res = typeof(A)[]
  for k in 1:ncols(K)
    M = zero_matrix(base_ring(A), nrows(A), ncols(A))
    for (l, v) in enumerate(cartind)
      M[v[2], v[1]] = K[l, k]
    end
    push!(res, M)
  end
  return res
end

function _basis_of_integral_commutator_algebra(A)
  linind = transpose(LinearIndices(size(A)))
  cartind = cartesian_product_iterator([1:x for x in size(A)])
  n = nrows(A)
  z = zero_matrix(FlintZZ, n^2, n^2)
  for i in 1:n
    for j in 1:n
      for k in 1:n
        z[linind[i, j], linind[i, k]] = FlintZZ(z[linind[i, j], linind[i, k]] + A[k, j])
        z[linind[i, j], linind[k, j]] = FlintZZ(z[linind[i, j], linind[k, j]] - A[i, k])
      end
    end
  end
  r, K = right_kernel(z)
  res = typeof(A)[]
  for k in 1:ncols(K)
    M = zero_matrix(base_ring(A), nrows(A), ncols(A))
    for (l, v) in enumerate(cartind)
      M[v[2], v[1]] = K[l, k]
    end
    push!(res, M)
  end
  return res
end

function _basis_of_integral_commutator_algebra(A, B)
  linind = transpose(LinearIndices(size(A)))
  cartind = cartesian_product_iterator([1:x for x in size(A)])
  n = nrows(A)
  z = zero_matrix(FlintZZ, n^2, n^2)
  for i in 1:n
    for j in 1:n
      for k in 1:n
        z[linind[i, j], linind[i, k]] = FlintZZ(z[linind[i, j], linind[i, k]] + A[k, j])
        z[linind[i, j], linind[k, j]] = FlintZZ(z[linind[i, j], linind[k, j]] - B[i, k])
      end
    end
  end
  r, K = right_kernel(z)
  res = typeof(A)[]
  for k in 1:ncols(K)
    M = zero_matrix(base_ring(A), nrows(A), ncols(A))
    for (l, v) in enumerate(cartind)
      M[v[2], v[1]] = K[l, k]
    end
    @assert M * A == B * M
    push!(res, M)
  end
  return res
end

function _induce_action(M, z::CommutatorAlgebra)
  res = []
  for i in 1:length(z.K)
    K = z.K[i]
    @show z.Eig[i]
    MEig = reduce(vcat, identity.(z.Eig[i]))
    @show MEig, M
    m = MEig * change_base_ring(K, M)
    @show m
    fl, x = can_solve(MEig, m, side = :left)
    @assert fl
    push!(res, x)
  end
  return identity.(res)
end

function _extract_bases_of_jordan_blocks(a, K, c = gen(K))
  aK = change_base_ring(K, a)
  J, B = jordan_normal_form(aK)
  @show J
  blocks = []
  d = degree(K)
  i = 1
  while i <= nrows(a)
    b, ro, i = _get_next_jordan_block(J, i)
    @show b, ro, i
    if b != c
      continue
    end

    if degree(b) == d
      push!(blocks, (b, ro))
    end
  end
  @show blocks
  res = []
  for b in blocks
    if length(b[2]) == 0
      continue
    end
    ro = b[2]
    v = [[B[ro[r], j] for j in 1:ncols(B)] for r in 1:length(ro)]
    @assert matrix(K, 1, ncols(B), v[end]) * aK == b[1] * matrix(K, 1, ncols(B), v[end])
    push!(res, v)
  end

  return res
end

function _get_next_jordan_block(A, i)
  @show A, i
  i0 = i
  bad = false
  while i < nrows(A)
    if !iszero(A[i0:i, (i + 1):nrows(A)]) || !iszero(A[(i + 1):nrows(A), i0:i])
      bad = true
      i += 1
    else
      break
    end
  end

  @show i

  for j in i0:(i-1)
    for k in i0:(j-1)
      if !iszero(A[j, k])
        return A[i0, i0], [], i + 1
      end
    end
  end

  @show i

  for j in i0:i
    if A[i0, i0] != A[j, j]
      return A[i0, i0], [], i + 1
    end
  end

  @show "here2"

  for j in i0:(i-1)
    if !isone(A[j, j + 1])
      return A[i0, i0], [], i + 1
    end
  end

  @show "here"
  for j in i0:i
    for j in (i0+1):(i-1)
      @show i, j, A[i, j]
      if !iszero(A[i, j])
        return A[i0, i0], [], i + 1
      end
    end
  end

  @show i

  i = i0
  res = [i]
  b = A[i, i]
  while i < nrows(A) && A[i, i + 1] == 1
    i += 1
    push!(res, i)
  end
  @show res
  return b, res, i + 1
end

function _get_next_jordan_block2(A)
  # first collect the pivot entries
  pivots = Int[]
  for j in 1:nrows(A)
    k = 1
    while iszero(A[j, k])
      k += 1
    end
    push!(pivots, k)
  end

  k = 1

  @show pivots
end

################################################################################
#
#  Lifting again
#
################################################################################

mutable struct Emat
  R
  n
  i
  j
  a
end

function matrix(e::Emat)
  z = identity_matrix(e.R, e.n)
  z[e.i, e.j] = e.a
  return z
end

inv(e::Emat) = Emat(e.R, e.n, e.i, e.j, -e.a)

function elementary_matrix(R, n, i, j, a)
  M = identity_matrix(R, n)
  M[i, j] = a
  return M
end

function _lift2(MM)
  M = deepcopy(MM)
  R = base_ring(M)
  S = _base_ring(R)
  @assert det(M) == 1
  k = 1
  n = nrows(M)
  left = []
  left2 = []
  while k < ncols(M)
    @assert det(M) == 1
    # scan the first column for a unit
    v = Int[euclid(M[i, k]) for i in k:n]
    @show v
    if isone(minimum(abs.(v)))
      l = findfirst(isone, v) + (k - 1)
      @show k, M
      @show isunit(M[l, k])
      @show M[l, k]
      if l != k
        @assert l isa Int
        @show l
        b = M[l, k]
        @show b
        @assert det(M) == 1
        E1 = elementary_matrix(R, n, k, l, one(R))
        E2 = elementary_matrix(R, n, l, k, -one(R))
        E3 = elementary_matrix(R, n, k, l, one(R))
        M = (E1 * E2 * E3) * M
        @assert det(M) == 1
        # push first E3
        push!(left, elementary_matrix(S, n, k, l, one(S)))
        @assert map_entries(R, left[end]) == E3
        push!(left2, Emat(S, n, k, l, one(S)))
        push!(left, elementary_matrix(S, n, l, k, -one(S)))
        push!(left2, Emat(S, n, l, k, -one(S)))
        @assert map_entries(R, left[end]) == E2
        push!(left, elementary_matrix(S, n, k, l, one(S)))
        push!(left2, Emat(S, n, k, l, one(S)))
        @assert map_entries(R, left[end]) == E1
      end
      @assert det(M) == 1
      @assert isunit(M[k, k])
      for j in (k+1):n
        if !iszero(M[j, n])
          N = deepcopy(M)
          q = divexact(M[j, k], M[k, k])
          E = elementary_matrix(R, n, j, k, -q)
          push!(left, elementary_matrix(S, n, j, k, lift(-q)))
          push!(left2, Emat(S, n, j, k, lift(-q)))
          @assert map_entries(R, left[end]) == E
          for o in 1:n
            M[j, o] = M[j, o] - q * M[k, o]
          end
          @assert E * N == M
        end
      end
      k += 1
      @show M
    else # no unit there # I could do this in one go by putting a unit in the first position
      _, l = findmin(abs.(v))
      l = l + (k - 1)
      @show l
      @show euclid(M[l, k])
      @show v
      for j in k:n
        if j == l
          continue
        end
        fl, _ = divides(M[j, k], M[l, k])
        if !fl
          @show M[j, k], M[l, k]
          N = deepcopy(M)
          @show euclid(M[j, k])
          q, r = divrem(M[j, k], M[l, k])
          @show euclid(r)
          E = elementary_matrix(R, n, j, l, -q)
          for o in 1:n
            M[j, o] = M[j, o] - q * M[l, o]
          end
          push!(left, elementary_matrix(S, n, j, l, lift(-q)))
          push!(left2, Emat(S, n, j, l, lift(-q)))
          @assert map_entries(R, left[end]) == E
          @assert E * N == M
          break
        end
      end
    end
    @assert det(M) == 1
  end
  println("M should now be lower triangular")
  @show M
  @assert det(M) == 1
  # Now M is lower triangular with units on the diagonal
  @show "here"
  for i in n:-1:2
    for j in (i - 1):-1:1
      # set jth row to jth row - M[k, j]/M[k, k] * jth row
      @show isunit(M[i, i])
      @show M[i, i]
      q = -divexact(M[j, i], M[i, i])
      @assert (-q) * M[i, i] == M[j, i]
      E = elementary_matrix(R, n, j, i, q)
      push!(left, elementary_matrix(S, n, j, i, lift(q)))
      push!(left2, Emat(S, n, j, i, lift(q)))
      @assert map_entries(R, left[end]) == E
      N = deepcopy(M)
      M[j, i] = M[j, i] + q * M[i, i]
      @assert E * N == M
    end
  end
  @assert det(M) == 1
  # Now sort the diagonal
  # We use this neat trick of Rosenberg, p. 60
  for i in 1:(n - 1)
    a = inv(M[i, i])
    ainv = M[i, i]
    # we multiply with (1,...1,a,a^-1,1,...1)
    E1 = elementary_matrix(R, n, i, i + 1, a)
    E2 = elementary_matrix(R, n, i + 1, i, -ainv)
    E3 = E1
    #E4 = identity_matrix(R, n)
    #E4[i, i] = 0
    #E4[i, i + 1] = -1
    #E4[i + 1, i + 1] = 0
    #E4[i + 1, i] = 1
    E5 = elementary_matrix(R, n, i, i + 1, -one(R))
    E6 = elementary_matrix(R, n, i + 1, i, one(R))
    E4 = E6

    #E4S = identity_matrix(S, n)
    #E4S[i, i] = 0
    #E4S[i, i + 1] = -1
    #E4S[i + 1, i + 1] = 0
    #E4S[i + 1, i] = 1
    E1S = elementary_matrix(S, n, i, i + 1, lift(a))
    E1SS = Emat(S, n, i, i + 1, lift(a))
    E2S = elementary_matrix(S, n, i + 1, i, lift(-ainv))
    E2SS = Emat(S, n, i + 1, i, lift(-ainv))
    E3S = E1S
    E3SS = E1SS
    E5S = elementary_matrix(S, n, i, i + 1, -one(S))
    E5SS = Emat(S, n, i, i + 1, -one(S))
    E6S = elementary_matrix(S, n, i + 1, i, one(S))
    E6SS = Emat(S, n, i + 1, i, one(S))
    E4S = E6S
    E4SS = E6SS
    push!(left, E6S)
    push!(left, E5S)
    push!(left, E4S)
    push!(left, E3S)
    push!(left, E2S)
    push!(left, E1S)
    push!(left2, E6SS)
    push!(left2, E5SS)
    push!(left2, E4SS)
    push!(left2, E3SS)
    push!(left2, E2SS)
    push!(left2, E1SS)

    M = (E1 * E2 * E3 * E4 * E5 * E6) * M
    @assert det(M) == 1
  end

  @assert isone(M)

  for i in 1:length(left)
    @assert matrix(left2[i]) == left[i]
  end

  return prod(matrix.(inv.(left2)))

  @assert map(R, prod(reverse(left))) * MM == 1

  return M, left
end

euclid(n::nmod) = gcd(n.data, modulus(parent(n)))
euclid(n::Nemo.fmpz_mod) = gcd(n.data, modulus(parent(n)))

function Base.divrem(n::T, m::T) where T <: Union{nmod,Nemo.fmpz_mod}
  @assert !iszero(m)
  R = parent(n)
  e = euclid(m)

  cp = coprime_base(fmpz[n.data, m.data, modulus(R)])::Array{fmpz, 1}

  q = Array{Tuple{fmpz, fmpz}, 1}()
  for i=1:length(cp)
    v = valuation(modulus(R), cp[i])::Int
    if v != 0
      pk = cp[i]^v
      nv = valuation(n.data % pk, cp[i])::Int
      mv = valuation(m.data % pk, cp[i])::Int
      if nv < mv
        push!(q, (pk, 0))
      else
        push!(q, (pk, divexact(n.data % pk, cp[i]^nv)))
      end
    end
  end
  qq =  R(crt([x[2] for x = q], [x[1] for x = q])::fmpz)::T
  rr = n-qq*m
  @assert n == qq*m+rr
  @assert rr == 0 || euclid(rr) < e
  return (qq,rr)::Tuple{T, T}
end

################################################################################
#
#  Generate interesting examples
#
################################################################################

function _jordan_block(R, n::Int, a)
  z = identity_matrix(R, n)
  for i in 1:(n - 1)
    z[i + 1, i] = a
  end
  return z
end

function _random_elementary_operations!(a; type = :rows)
  tr = type == :columns
  @assert issquare(a)
  n = nrows(a)
  d = det(a)
  if n == 1
    return a
  end
  i = rand(1:n)
  j = rand(1:n)
  while j == i
    j = rand(1:n)
  end

  r = rand(base_ring(a), 1:10)
  for k in 1:n
    if tr
      a[k, i] = a[k, i] + r * a[k, j]
    else
      a[i, k] = a[i, k] + r * a[j, k]
    end
  end
  @assert d == det(a)
  return a
end

function _random_sln(R, n; num_op = 10)
  a = identity_matrix(R, n)
  for i in 1:num_op
    _random_elementary_operations!(a; type = rand() < 0.5 ? :rows : :columns)
  end
  return a
end

function _random_matrix(R, block_shape, eigval_range = -10:10)
  matrices = dense_matrix_type(R)[]
  for r in block_shape
    push!(matrices, _jordan_block(R, r, rand(eigval_range)))
  end
  return diagonal_matrix(matrices)
end

function _similarity_test_setup(R, n)
  block_shape = Int[]
  nn = n
  while !iszero(nn)
    r = rand(1:nn)
    push!(block_shape, r)
    nn = nn - r
  end
  b = block_shape
  A = _random_matrix(R, block_shape)
  z = _random_sln(R, n)
  @assert isone(det(z))
  zinv = inv(z)
  B = z * A * zinv
  return A, B
end

################################################################################
#
#  Redoing the Jordan decomposition thing
#
################################################################################

function _get_morphism(A::fmpq_mat)
  Qx = Hecke.Globals.Qx
  x = gen(Qx)
  Ax = x - change_base_ring(Qx, A)
  S, _, T = snf_with_transform(Ax)
  return T, diagonal(S)
end


text:
## Primary testing file for CI, slower and more thorough testing, including two tests
## that are not run in `runtests.jl`:

# (1) We only run tests for the SCC here since they are fairly slow
# (2) We only run tests for the extensions here since they too are slow, and
#     also clean out the output folders

using Test
using Mimi

# include the main model
include("../src/main_model.jl")

@testset "MimiPAGE2020-main" begin
    include("runtests.jl")
end

@testset "MimiPAGE2020-extra" begin
    include("test_scc.jl") # Takes very long to run.
    include("test_extensions.jl") # NB will currently clean out the 'output' folder! - tests number of model extension output files for PAGE2020 update
end


text:
using Documenter, ICDataServer

makedocs()


text:
# UnitCommitment.jl: Optimization Package for Security-Constrained Unit Commitment
# Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
# Released under the modified BSD license. See COPYING.md for more details.

using UnitCommitment, LinearAlgebra, Cbc, JuMP, JSON

@testset "build_model" begin
    instance = UnitCommitment.read_benchmark("test/case14")
    for line in instance.lines, t in 1:4
        line.normal_flow_limit[t] = 10.0
    end
    optimizer = optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 0)
    model = UnitCommitment.build_model(
        instance = instance,
        optimizer = optimizer,
        variable_names = true,
    )
    @test name(model[:is_on]["g1", 1]) == "is_on[g1,1]"

    # Optimize and retrieve solution
    UnitCommitment.optimize!(model)
    solution = UnitCommitment.solution(model)

    # Write solution to a file
    filename = tempname()
    UnitCommitment.write(filename, solution)
    loaded = JSON.parsefile(filename)
    @test length(loaded["Is on"]) == 6

    # Verify solution
    @test UnitCommitment.validate(instance, solution)

    # Reoptimize with fixed solution
    UnitCommitment.fix!(model, solution)
    UnitCommitment.optimize!(model)
    @test UnitCommitment.validate(instance, solution)
end


text:
module Dp

export fibonacci

struct FibCACHE
  INT_CACHE::Vector{Int}
  BIGINT_CACHE::Vector{BigInt}
end

const fib_cache = FibCACHE(Int[1,1], BigInt[])

function (fib::FibCACHE)(n)
  if n <= length(fib.INT_CACHE)
    return fib.INT_CACHE[n]
  elseif n > 92
    bn = n - 92
    if bn <= length(fib.BIGINT_CACHE)
      return fib.BIGINT_CACHE[bn]
    else
      fnm1 = fibonacci(n-1) |> BigInt
      fnm2 = fibonacci(n-2) |> BigInt
      return push!(fib.BIGINT_CACHE, fnm1 + fnm2)[bn]
    end
  else
    fnm1 = fibonacci(n-1)
    fnm2 = fibonacci(n-2)
    return push!(fib.INT_CACHE, fnm1 + fnm2)[n]
  end
end

function fibonacci(n)
  global fib_cache
  n == 0 && return 0
  return fib_cache(n)
end

end


text:
# Meta Parsers
# ============

@testset "Named" begin
  parser = EP.Named{:hi}(EP.anysymbol)
  parsed = parse_expr(parser, :symbol)
  @test parsed isa EP.Named{:hi, Base.Symbol}
  test_closure(parser, :symbol)
end

@testset "Indexed" begin
  indexed_parser = EP.Indexed() do dict
    EP.Expr(quote
      a = $(dict[:a] = EP.Isa(Int))
      b = $(dict[:b] = EP.anysymbol)
    end)
  end
  expr = quote
    a = 4
    b = a
  end
  test_closure(indexed_parser, expr)
  indexed_parsed = parse_expr(indexed_parser, expr)
  @test indexed_parsed isa EP.Indexed{EP.Expr_Parsed}
  @test collect(keys(indexed_parsed)) == [:a; :b]
  @test collect(values(indexed_parsed)) == [4, :a]
  @test indexed_parsed[:a] == 4
  @test indexed_parsed[:b] == :a

  indexed_parsed[:a] = 5
  @test indexed_parsed[:a] == 5
  to_expr(indexed_parsed)  # looks good, but hard to test automatically

  @test_throws EP.ParseError parse_expr(indexed_parser, quote
    a = 4.0
    b = a
  end)
end


text:
# signal_haar.jl
##########################
fs=Int(1000);
Ts=1/fs;
N=Int(12000);
T=N/fs;
#t=(0:Int(N-1))*Ts;

x=zeros(Float64, N)
for m=1:N
    x[m]=sin(2*Ts*(m-1)*(Ts*(m-1)-3)*(Ts*(m-1)-6)*(Ts*(m-1)-9)*(Ts*(m-1)-12))
end

using DSP, Wavelets

y1=dwt(x, wavelet(WT.haar), 1);
y2=dwt(x, wavelet(WT.haar), 2);
y3=dwt(x, wavelet(WT.haar), 3);
y4=dwt(x, wavelet(WT.haar), 4);
y5=dwt(x, wavelet(WT.haar), 5);

y1a=zeros(N);
y2a=zeros(N);
y3a=zeros(N);
y4a=zeros(N);
y5a=zeros(N);
y1a[1:Int(N/2)]=y1[1:Int(N/2)];
y2a[1:Int(N/4)]=y2[1:Int(N/4)];
y3a[1:Int(N/8)]=y3[1:Int(N/8)];
y4a[1:Int(N/16)]=y4[1:Int(N/16)];
y5a[1:Int(N/32)]=y5[1:Int(N/32)];
y1d=y1-y1a;
y2d=y2-y2a;
y3d=y3-y3a;
y4d=y4-y4a;
y5d=y5-y5a;

x1a=idwt(y1a, wavelet(WT.haar), 1);
x1d=idwt(y1d, wavelet(WT.haar), 1);
x2a=idwt(y2a, wavelet(WT.haar), 2);
x2d=idwt(y2d, wavelet(WT.haar), 2);
x3a=idwt(y3a, wavelet(WT.haar), 3);
x3d=idwt(y3d, wavelet(WT.haar), 3);
x4a=idwt(y4a, wavelet(WT.haar), 4);
x4d=idwt(y4d, wavelet(WT.haar), 4);
x5a=idwt(y5a, wavelet(WT.haar), 5);
x5d=idwt(y5d, wavelet(WT.haar), 5);

using Plots

# x1a
p1a=plot(x1a,
    xaxis="n = 0, 1, 2, ..., N-1",
    yaxis="amplitude",
    title="Level 1 Approximation",
    grid=false,
    legend=false)
plot!(xticks = ([0 6000 12000;], [0 6000 12000]))
plot!(yticks = ([-1 0 1;], [-1 0 1]))
# x1d
p1d=plot(x1d,
    xaxis="n = 0, 1, 2, ..., N-1",
    yaxis="amplitude",
    title="Level 1 Detail",
    grid=false,
    legend=false)
plot!(xticks = ([0 6000 12000;], [0 6000 12000]))
plot!(yticks = ([-1 0 1;], [-1 0 1]))
# x1
plot(p1a, p1d, layout = (1,2))
savefig("signal_haar_level_1.png")

# x2a
p2a=plot(x2a,
    xaxis="n = 0, 1, 2, ..., N-1",
    yaxis="amplitude",
    title="Level 2 Approximation",
    grid=false,
    legend=false)
plot!(xticks = ([0 6000 12000;], [0 6000 12000]))
plot!(yticks = ([-1 0 1;], [-1 0 1]))
# x2d
p2d=plot(x2d,
    xaxis="n = 0, 1, 2, ..., N-1",
    yaxis="amplitude",
    title="Level 2 Detail",
    grid=false,
    legend=false)
plot!(xticks = ([0 6000 12000;], [0 6000 12000]))
plot!(yticks = ([-1 0 1;], [-1 0 1]))
# x1
plot(p2a, p2d, layout = (1,2))
savefig("signal_haar_level_2.png")

# x3a
p3a=plot(x3a,
    xaxis="n = 0, 1, 2, ..., N-1",
    yaxis="amplitude",
    title="Level 3 Approximation",
    grid=false,
    legend=false)
plot!(xticks = ([0 6000 12000;], [0 6000 12000]))
plot!(yticks = ([-1 0 1;], [-1 0 1]))
# x3d
p3d=plot(x3d,
    xaxis="n = 0, 1, 2, ..., N-1",
    yaxis="amplitude",
    title="Level 3 Detail",
    grid=false,
    legend=false)
plot!(xticks = ([0 6000 12000;], [0 6000 12000]))
plot!(yticks = ([-1 0 1;], [-1 0 1]))
# x1
plot(p3a, p3d, layout = (1,2))
savefig("signal_haar_level_3.png")

# x4a
p4a=plot(x4a,
    xaxis="n = 0, 1, 2, ..., N-1",
    yaxis="amplitude",
    title="Level 4 Approximation",
    grid=false,
    legend=false)
plot!(xticks = ([0 6000 12000;], [0 6000 12000]))
plot!(yticks = ([-1 0 1;], [-1 0 1]))
# x4d
p4d=plot(x4d,
    xaxis="n = 0, 1, 2, ..., N-1",
    yaxis="amplitude",
    title="Level 4 Detail",
    grid=false,
    legend=false)
plot!(xticks = ([0 6000 12000;], [0 6000 12000]))
plot!(yticks = ([-1 0 1;], [-1 0 1]))
# x1
plot(p4a, p4d, layout = (1,2))
savefig("signal_haar_level_4.png")

# x5a
p5a=plot(x5a,
    xaxis="n = 0, 1, 2, ..., N-1",
    yaxis="amplitude",
    title="Level 5 Approximation",
    grid=false,
    legend=false)
plot!(xticks = ([0 6000 12000;], [0 6000 12000]))
plot!(yticks = ([-1 0 1;], [-1 0 1]))
# x5d
p5d=plot(x5d,
    xaxis="n = 0, 1, 2, ..., N-1",
    yaxis="amplitude",
    title="Level 5 Detail",
    grid=false,
    legend=false)
plot!(xticks = ([0 6000 12000;], [0 6000 12000]))
plot!(yticks = ([-1 0 1;], [-1 0 1]))
# x1
plot(p5a, p5d, layout = (1,2))
savefig("signal_haar_level_5.png")


text:
# test partial pose3 constraints and evaluation

using Statistics
using RoME
# , Distributions
# using IncrementalInference, TransformUtils
using Test

global N=50
global fg = initfg()

global v1 = addVariable!(fg,:x1, Pose3, N=N) # 0.001*randn(6,N)
global f0 = addFactor!(fg, [:x1;], PriorPose3(MvNormal(zeros(6),1e-2*Matrix{Float64}(LinearAlgebra.I, 6,6))))

global sigx = 0.01
global sigy = 0.01
global sigth = 0.0281
global mu1 = [0.0;0.0; -10.0]
global prpz = PartialPriorRollPitchZ(
  MvNormal( mu1[1:2], [sigx 0.0; 0.0 sigy]^2 ),
  Normal( mu1[3], sigth )
)

global mu2 = [20.0,5.0,pi/2]
global xyy = PartialPose3XYYaw(
  MvNormal(
    mu2[1:2],
    [sigx 0.0; 0.0 sigy]^2
  ),
  Normal( mu2[3], sigth )
)

global v2 = addVariable!(fg,:x2, Pose3, N=N) # randn(6,N)

global f1 = addFactor!(fg, [:x2], prpz, autoinit=false)
global f2 = addFactor!(fg, [:x1;:x2], xyy, autoinit=false)

# ls(fg, :x2)

@testset "test PartialPriorRollPitchZ evaluations" begin

# ensure that at least the first pose is already initialized
doautoinit!(fg, :x1)
@test isInitialized(fg, :x1)

X1pts = getVal(fg, :x1)
@test sum(isnan.(X1pts)) == 0

ppts = approxConv(fg, :x1x2f1,:x2)
@test sum(isnan.(ppts)) == 0

ensureAllInitialized!(fg)
@test isInitialized(fg, :x2)

# get values and ensure that a re-evaluation produces consistent results
global X2pts = getVal(fg, :x2)
@test sum(isnan.(X2pts)) == 0

# check that only partial states are updated
global pts = IIF.approxConv(fg, :x2f1, :x2, N=N)

global newdims = collect(getData(f1).fnc.usrfnc!.partial)
global olddims = setdiff(collect(1:6), newdims)

@test size(pts, 1) == 6
@test size(pts, 2) == N

# check that untouched dimensions (not in the partial list) truely remain untouched
@test norm(X2pts[olddims,:] - pts[olddims,:]) < 1e-10

# check that the prior new dims are updated to new and correct values
# @show Statistics.mean(pts,dims=2)[newdims]
@test sum(abs.(Statistics.mean(pts,dims=2)[newdims]-mu1[[3;1;2]]) .< [0.5; 0.1; 0.1]) == 3

# ensure a forced re-evaluatoin
@test norm(X2pts[newdims,:] - pts[newdims,:]) < 1.0

# memcheck that the exact same values are used
@test norm(X2pts - getVal(v2)) < 1e-10

end

@testset "test residual function of PartialPose3XYYaw" begin

  global res = zeros(3)
  global idx = 1
  global meas = getSample(xyy)
  global xi = zeros(6,1)
  global xja = zeros(6,1)
  xyy(res, nothing, idx, meas, xi, xja)
  @test abs(res[1]-mu2[1]) < 0.3
  @test abs(res[2]-mu2[2]) < 0.3
  @test abs(res[3]-mu2[3]) < 0.2

  global xjb = zeros(6,1)
  xjb[collect(xyy.partial),1] = mu2
  global res = zeros(3)
  xyy(res, nothing, idx, meas, xi, xjb)
  @test 0.0 < norm(res) < 0.3

  global meas = getSample(xyy,100)
  @test norm(Statistics.std(meas[1],dims=2) - [0.01;0.01;0.002]) < 0.05

end

@testset "test PartialPose3XYYaw evaluations" begin

# get existing and predict new
global X2pts = getVal(fg, :x2)
global pts = approxConv(fg, :x1x2f1, :x2, N=N)

# find which dimensions are and and are not updated by XYYaw partial
global newdims = collect(getData(f2).fnc.usrfnc!.partial)
global olddims = setdiff(collect(1:6), newdims)

# check the number of points are correct
@test size(pts, 1) == 6
@test size(pts, 2) == N

# ensure the unchanged dimensions actually remain unchanged
@test norm(X2pts[olddims,:] - pts[olddims,:]) < 1e-10

for i in 1:N
    pts[6,i] = wrapRad(pts[6,i])
end

# ensure the newly updated values match what is specified in mu2
@show mu2 # mu2 is used for XYYaw
# trying to compare world and body frame values -- not right!
@show Statistics.mean(pts[newdims,:],dims=2)
@test sum(abs.(Statistics.mean(pts[newdims,:],dims=2)-mu2) .< [0.7;0.7;0.2]) == 3

# ensure a re-evaluation of the partial factor updates the partial variable dimensions correclty
@test norm(X2pts[newdims,:] - pts[newdims,:]) < 1.0

# ensure that memory pointers are working correctly
global memcheck = getVal(v2)
@test norm(X2pts - memcheck) < 1e-10

end

@testset "test predictbelief with two functions" begin
  global val, = predictbelief(fg, :x2, ls(fg, :x2), N=N)

  for i in 1:N
    val[6,i] = wrapRad(val[6,i])
  end

  @test size(val, 1) == 6
  @test size(val, 2) == N

  global estmu1mean = Statistics.mean(val[collect(getData(f1).fnc.usrfnc!.partial),:],dims=2)
  global estmu2mean = Statistics.mean(val[collect(getData(f2).fnc.usrfnc!.partial),:],dims=2)

  @test sum(abs.(estmu1mean - mu1[[3;1;2]]) .< [0.7; 0.1; 0.1]) == 3
  @test sum(abs.(estmu2mean - mu2) .< [0.7; 0.7; 0.15] ) == 3

  global memcheck = getVal(v2)
  @test 1e-10 < norm(val - memcheck)
end

#


text:
data = openBF.loadSimulationFiles("test.yml")

blood = openBF.buildBlood(data["blood"])
vessels, edges = openBF.buildArterialNetwork(data["network"], blood, 100)
openBF.makeResultsFolder(data)
Ccfl = data["solver"]["Ccfl"]
heart = vessels[1].heart

A = vessels[1].A[1]
A0 = vessels[1].A0[1]
beta = vessels[1].beta[1]
gamma = vessels[1].gamma[1]
Pext = vessels[1].Pext

p = openBF.pressure(A, A0, beta, Pext)
@test p == 1e4
@test p == openBF.pressure(sqrt(A/A0), beta, Pext)

c = openBF.waveSpeed(A, gamma)
@test isapprox(c, 6.33, atol=1e-2)

dt = openBF.calculateDeltaT(vessels, Ccfl)
@test isapprox(dt, 1e-4, atol=1e-4)

@test_nowarn openBF.solveModel(vessels, edges, blood, dt, 0.0)

q = vessels[1].Q[1]
p = vessels[1].P[end]
@test_nowarn openBF.solveVessel(vessels[1], blood, dt, 0.0)
@test q != vessels[1].Q[1]
@test p != vessels[1].P[end]

@test_nowarn openBF.solveOutlet(1, vessels[1], blood, vessels, edges, dt)

v = vessels[1]
@test_nowarn openBF.computeLimiter(v, v.vA, v.invDx, v.dU, v.slopesA)
@test_nowarn openBF.computeFlux(v, v.Al, v.Ql, v.Fl)

@test openBF.maxMod(1.0, 0.0) == 1.0
@test openBF.maxMod(0.0, 0.5) == 0.5

@test openBF.minMod(-1.0, 1.0) == 0.0
@test openBF.minMod(1.0, -1.0) == 0.0
@test openBF.minMod(0.5, 1.0) == 0.5
@test openBF.minMod(1.0, 0.5) == 0.5

@test_nowarn openBF.superBee(v, v.dU, v.slopesA)

cd("..")
rm("test_results", recursive=true)


text:
module MaximumEntropyBasisFunctions

# export CreateMaxEntBasis, getBasisDerivative, EvaluateMaxEntBasis

# Maximum entropy basis functions
# ===============================

using  LineSearches, Optim, NLsolve, NearestNeighbors, LinearAlgebra

struct MaxEntBasis
	x::Matrix{Real};
	β::Real; # Constant β. Future version we can incorporate x_i dependent β
end

"""
CreateMaxEntBasis(x::Matrix; nn=3, γ=1)::MaxEntBasis

# DESCRIPTION

	Creates MaxEntBasis object with given points x, and tuning parameters nn, flag.
	This object is used to generate basis and derivatives.

# INPUT

	x::Matrix of size d x N, where d is the dimension of state space, N is the number of points
	nn is an integer to define number of nearest neighbours -- tune local vs global property
	flag tunes how distances from nn neighbours are used to define locality of maximum entropy basis functions
		flag = 0: minimum distance (default)
		flag = 1: maximum distance
		flag = 2: average distance
	a is scaling of distance based β, detault = 1

# OUTPUT

	MaxEntBasis	object with fields x::Matrix{Real}, and β::Vector{Real}, where β defines the locality of the gaussian RBs.
"""
function CreateMaxEntBasis(x::Matrix;nn=3,γ=1)::MaxEntBasis
	# 1. We need to compute β, which is needed for computing the maxent basis functions.

	kdtree = KDTree(x)
	idxs, dists = knn(kdtree, x, nn, true); # Compute knn for all points in the data

	# # Compute β according to flag -- use list comprehension to do this efficiently
	# if flag == 0 # Minimum
	# 	b = [d[2] for d in dists];
	# elseif flag == 1 # Maximum
	# 	b = [d[end] for d in dists];
	# elseif flag == 2 # Average
	# 	b = [sum(d[2:end])/(length(d)-1) for d in dists]; # Skip the first element, which is distance with self.
	# end
	# β = γ ./ b; # See paper

	D = vcat(dists'...);
	D1 = vec(D[:,2:nn]);
	avgD = sum(D1)/length(D1); # Average nn distance.
	β = γ/avgD; # See paper

	return(MaxEntBasis(x,β));
end

function EvaluateMaxEntBasis(objϕ::MaxEntBasis, x::Array{T,1})::Array{T,2} where T <: Real
	xTilde = objϕ.x .- x;
	nBasis = size(xTilde,2);
	λ = getLambda(x,objϕ);
	ϕ = meBasis(λ,x,objϕ)'; # Return row vector
	return ϕ;
end

function priorFunction(x,β)
    f = exp(-β*norm(x)^2);
	return f;
end

function getLambda(x,objϕ)
	function costFunction(λ)
		phi = meBasis(λ,x,objϕ);
		nBasis = length(phi);
		d = length(x);

		J = zeros(d);
		xTilde = x .- objϕ.x;

		for i in 1:nBasis
			J += xTilde[:,i]*phi[i];
		end
		return J'*J
	end

	d = length(x)
	initial_x = zeros(d);
	sol = optimize(costFunction,initial_x,Newton());
	λ = Optim.minimizer(sol);
	# println("λ: $λ")
	return(λ);
end

# Generate maximum entropy-basis functions
# ========================================
function meBasis(λ, x, objϕ)
	xTilde = x .- objϕ.x;
	nBasis = size(xTilde,2);
	Z = zeros(nBasis);

	for i in 1:nBasis
		m = priorFunction(xTilde[:,i],objϕ.β);
		Z[i] = m*exp(-λ'*xTilde[:,i]);
	end
	phi = Z./sum(Z);
	return phi;
end

# Compute gradient of maxent basis function
# =================================================================
# See following paper for technical details
#	Local maximum-entropy approximation schemes:
#	A seamless bridge between finite elements and meshfree methods
#   by M. Arroyo and M. Ortiz
# -----------------------------------------------------------------
function getBasisDerivative(objϕ,xEval,ϕ)
	xTilde = objϕ.x .- xEval;
	d, nBasis = size(xTilde);

	J = zeros(d,d);
	r = zeros(d);

	for i in 1:nBasis
		xTilde = (xEval-objϕ.x[:,i]);
		r .+= ϕ[i]*xTilde;
		J .+= ϕ[i]*xTilde*xTilde';
	end
	J = J - r*r';
	# invJ = inv(J); # replace
	invJ = J\I(d);
	∇ϕ = zeros(d,nBasis);
	for i in 1:nBasis
		∇ϕ[:,i] = -ϕ[i]*invJ*(xEval-objϕ.x[:,i]);
	end
	return ∇ϕ
end
end # module

max_stars_repo_name: isrlab/MaximumEntropyBasisFunctions

max_issues_repo_name: isrlab/MaximumEntropyBasisFunctions

max_forks_repo_name: isrlab/MaximumEntropyBasisFunctions


text:
module TestMeasures

using MLJFair, Test

include(joinpath("..", "data", "data.jl"))

include("calcmetrics.jl")
include("boolmetrics.jl")

end  # module


text:
#
# Factors Join
#

"""
    join(op, ϕ1::Factor, ϕ2::Factor)

Joins two `Factors` with an outer join using `op`

An outer join returns a `Factor` with the union of the two dimensions
The two factors are combined with `Base.broadcast(op, ...)`
"""
function Base.join(op, ϕ1::Factor, ϕ2::Factor)
    # avoid all the broadcast overhead with a larger array (ideally)
    # more useful for edge cases where one is (or both are) singleton
    if length(ϕ1) < length(ϕ2)
        ϕ2, ϕ1 = ϕ1, ϕ2
        op_new = (x, y) -> op(y, x)
    else
        op_new = op
    end

    common_names = intersect(names(ϕ1), names(ϕ2))

    index_common1 = indexin(common_names, names(ϕ1))
    index_common2 = indexin(common_names, names(ϕ2))

    ϕ1.dimensions[index_common1] != ϕ2.dimensions[index_common2] &&
        throw(DimensionMismatch("Common dimensions must have same support"))

    index_unique1 = setdiff(1:ndims(ϕ1), index_common1)
    index_unique2 = setdiff(1:ndims(ϕ2), index_common2)

    # the first dimensions are all from ϕ1
    new_dims = vcat(ϕ1.dimensions, ϕ2.dimensions[index_unique2])

    if ndims(ϕ2) != 0
        # permute the common dimensions of ϕ2 to the beginning
        #  in the order that they appear in ϕ1 (and new_dims)
        perm = vcat(index_common2, index_unique2)
        temp = permutedims(ϕ2.potential, perm)

        # reshape by lining up the common dims in ϕ2 with those in ϕ1
        size_unique2 = size(ϕ2)[index_unique2]
        # set those dims to have dimension 1 for data in ϕ2
        reshape_lengths = vcat(size(ϕ1)..., size_unique2...)
        #new_pot = duplicate(ϕ1.potential, size_unique2)
        new_pot = Array{Float64}(reshape_lengths...)
        reshape_lengths[index_unique1] = 1
        temp = reshape(temp, (reshape_lengths...))
    else
        new_pot = similar(ϕ1.potential)
        temp = ϕ2.potential
    end

    # ndims(ϕ1) == 0 implies ndims(ϕ2) == 0 b/c len(ϕ1) >= len(ϕ2)
    if ndims(ϕ1) == 0
        new_pot = op_new(ϕ1.potential[1], temp[1])
    else
        broadcast!(op_new, new_pot, ϕ1.potential, temp)
    end

    return Factor(new_dims, new_pot)
end

*(ϕ1::Factor, ϕ2::Factor) = join(*, ϕ1, ϕ2)
/(ϕ1::Factor, ϕ2::Factor) = join(/, ϕ1, ϕ2)
+(ϕ1::Factor, ϕ2::Factor) = join(+, ϕ1, ϕ2)
-(ϕ1::Factor, ϕ2::Factor) = join(-, ϕ1, ϕ2)

# code for an inner join, if its ever done
# error("not supported")
# new_dims = getdim(ϕ1, common)
#
# if isempty(common)
#     # weird magic for a zero-dimensional array
#     v_new = squeeze(zero(eltype(ϕ1), 0), 1)
# else
#     if reducehow == nothing
#         throw(ArgumentError("Need reducehow"))
#     end
#
#     inds1 = (findin(ϕ1.dimensions, common)...)
#     inds2 = (findin(ϕ2.dimensions, common)...)
#
#     if v0 != nothing
#         v1_new = squeeze(reducedim(op, ϕ1.v, inds1, v0), inds)
#         v2_new = squeeze(reducedim(op, ϕ2.v, inds2, v0), inds)
#     else
#         v1_new = squeeze(reducedim(op, ϕ1.v, inds1), inds)
#         v2_new = squeeze(reducedim(op, ϕ2.v, inds2), inds)
# new_m = zeros(eltype(ϕ1.v), reshape_lengths)
#     end
#     v_new = op(v1_new, v2_new)
# end


text:
using BlackBoxOptim, Plots, Random

"""
    plot_crossover(xover::CrossoverOperator, parents::Matrix{Float64};
                   n=5000, shuffle_parents=false)

Plot the distribution of offsprings for the crossover operator `xover`
and given `parents`.
"""
function plot_crossover(xover::CrossoverOperator, parents::Matrix{Float64};
                        n=5000, shuffle_parents=false)
    @assert size(parents, 2) == numparents(xover)
    parents_pop = FitPopulation(parents, NaN)
    parent_ixs = collect(1:numparents(xover))
    children_mtx = hcat([hcat(apply!(xover, [fill!(Individual(undef, size(parents, 1)), NaN) for _ in 1:numchildren(xover)],
                              zeros(Int, numchildren(xover)), parents_pop,
                              shuffle_parents ? shuffle(parent_ixs) : parent_ixs)...) for _ in 1:n]...)
    plot(view(children_mtx, 1, :), view(children_mtx, 2, :),
         seriestype=:scatter, title=string(typeof(xover)), label="children",
         markersize=1.5, markerstrokewidth=0, markercolor=colorant"orange")
    scatter!(view(parents_pop.individuals, 1, :), view(parents_pop.individuals, 2, :),
             markersize=5, markercolor=colorant"purple", label="parents")
end

plot_crossover(SimplexCrossover{3}(1.2),
               [[0.5, 0.0] [-0.2, 0.0] [0.0, 1.0]])

plot_crossover(SimulatedBinaryCrossover(0.2, 5.0),
               [[1.0, 0.2] [-0.2, 1.0]])

plot_crossover(UnimodalNormalDistributionCrossover{2}(0.5, 0.05),
               [[0.5, 0.0] [1.0, 3.0]])

plot_crossover(UnimodalNormalDistributionCrossover{3}(0.5, 0.05),
               [[0.5, 0.0] [1.0, 3.0] [3.0, 1.0]])

plot_crossover(ParentCentricCrossover{3}(0.1, 0.2),
               [[0.5, 0.0] [1.0, 3.0] [3.0, 1.0]])

plot_crossover(ParentCentricCrossover{3}(0.1, 0.2),
               [[0.5, 0.0] [1.0, 3.0] [3.0, 1.0]], shuffle_parents=true)


text:
function comparefields(dataobject1, dataobject2)

    pr = propertynames(dataobject1);
    println("Unequal/not defined fields: ")

    println()
    println("Legend:")
    println("do1 = dataobject1, do2 = dataobject2")
    rr =96; gg=168; bb=48 # green
    print("- type in ", "\e[38;2;$rr;$gg;$bb;249m","green", "\e[0m", " => unequal fields are fine as long as their fields are equal or defined  \n")

    #rr =175; gg=120; bb=2 # orange
    #print("- type in ", "\e[38;2;$rr;$gg;$bb;249m","orange", "\e[0m", " => marks dictionaries \n")

    rr =200; gg=75; bb=75 # red
    print("- type in ", "\e[38;2;$rr;$gg;$bb;249m","red", "\e[0m", " => undefined fields = ndf \n")
    println("(if do1 & do2 are both undefined, than everything is fine! ;)")
    println("=========================================")

    for i in pr
        level = 1
        j1, j2 = checkfield(dataobject1, dataobject2, i, level)
        comparefield(i, level, j1, j2)

        if length(propertynames(j1)) != 0 # if subfield
                for j in propertynames(j1)
                    level =2
                    k1, k2 = checkfield(j1, j2, j, level)

                    if typeof(k1) !== Dict{Any,Any} # distinguish dicts and "normal" fields

                        if length(propertynames(k1)) != 0 # if subfield
                            for k in propertynames(k1)
                                level =3
                                l1, l2 = checkfield(k1, k2, k, level)
                                comparefield(k, level, l1, l2)

                                if length(propertynames(l1)) != 0 # if subfield
                                    for l in propertynames(l1)
                                        level =4
                                        m1, m2 = checkfield(l1, l2, l, level)
                                        comparefield(l, level, m1, m2)
                                    end
                                end
                            end
                        end

                    else # handle dicts

                        if length(k1) !== 0
                            for p in keys(k1)
                                println("----",p)
                            end
                        else # compare empty dicts
                            comparefield(k1, level, k1, k2)
                        end
                    end # distinguish dicts and "normal" fields

                end # j: forloop

        end # if subfield
    end # i: forloop

end

function checkfield(dataobject1, dataobject2, field, level)
    j1 = :nofield
    if isdefined(dataobject1, field)
        j1 = getfield(dataobject1, field)

    else
        rr =200; gg=75; bb=75 # red
        if level == 1
            print("\e[1m l$level: ", "\e[0m",field)
        elseif level == 2
            print("\e[1m \t l$level: ", "\e[0m",field)
        elseif level == 3
            print("\e[1m \t \t l$level: ", "\e[0m",field)
        end
        print(" -> \e[38;2;$rr;$gg;$bb;249m ndf do1 \e[0m \n")
    end

    j2 = :nofield
    if isdefined(dataobject2, field)
        j2 = getfield(dataobject2, field)
    else
        rr =200; gg=75; bb=75 # red
        if level == 1
            #print("\e[1m","\t \t l", field_levels, ": ", "\e[0m", j, " = ", ctype)
            print("\e[1m l$level: ", "\e[0m",field)
        elseif level == 2
            print("\e[1m \t l$level: ", "\e[0m",field)
        elseif level == 3
            print("\e[1m \t \t l$level: ", "\e[0m",field)
        end
        print(" -> \e[38;2;$rr;$gg;$bb;249m ndf do2 \e[0m \n")
    end

    return j1, j2
end

function comparefield(field, level, j1, j2)
     if j1 != :nofield || j2 != :nofield
            if j1 != j2

                rr =96; gg=168; bb=48 # green
                if level == 1
                    println("\e[1m l$level: ", "\e[0m",field, " -> \e[38;2;$rr;$gg;$bb;249m unequal \e[0m")
                elseif level ==2
                    println("\e[1m \t l$level: ", "\e[0m",field, " -> \e[38;2;$rr;$gg;$bb;249m unequal \e[0m")
                elseif level ==3
                    println("\e[1m \t \t l$level: ", "\e[0m",field, " -> \e[38;2;$rr;$gg;$bb;249m unequal \e[0m")
                end
            end
        #end
    end
end


text:
using SparseArrays
include("../src/Lar.jl")
LarA = Lar.Arrangement

triforce = Lar.Model([
    0.0  4.0  8.0  4.0  2.0  6.0  4.0 -4.0 12.0  4.0  3.0  5.0
    0.0  6.0  0.0  0.0  3.0  3.0 -6.0  6.0  6.0  1.0  2.5  2.5
])
triforce.T[1] = SparseArrays.sparse(Array{Int8, 2}([
    [1 1 0 0 0 0 0 0 0 0 0 0]
    [0 1 1 0 0 0 0 0 0 0 0 0]
    [1 0 1 0 0 0 0 0 0 0 0 0]
    [0 0 0 1 1 0 0 0 0 0 0 0]
    [0 0 0 0 1 1 0 0 0 0 0 0]
    [0 0 0 1 0 1 0 0 0 0 0 0]
    [0 0 0 0 0 0 1 1 0 0 0 0]
    [0 0 0 0 0 0 0 1 1 0 0 0]
    [0 0 0 0 0 0 1 0 1 0 0 0]
    [0 0 0 0 0 0 0 0 0 1 1 0]
    [0 0 0 0 0 0 0 0 0 0 1 1]
    [0 0 0 0 0 0 0 0 0 1 0 1]
]))
σ = convert(Lar.Chain, SparseArrays.sparse([
    1; 1; 1; 0; 0; 0; 0; 0; 0; 0; 0; 0
]))

model, sigma, edge_map = LarA.planar_arrangement_1(triforce, σ, true, false)
cleanedM, cleanedEM = LarA.cleandecomposition(model, sigma, edge_map)
bicon_comps = Lar.Arrangement.biconnected_components(cleanedM.T[1])
arrangedM = LarA.planar_arrangement_2(cleanedM, bicon_comps)


text:
# test/runtests.jl

using MIRT

test_all_mirt()


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule libfdk_aac_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("libfdk_aac")
JLLWrappers.@generate_main_file("libfdk_aac", UUID("f638f0a6-7fb0-5443-88ba-1cc74229b280"))
end  # module libfdk_aac_jll


text:
using Interact
using LightGraphs
using VisJS

# Little slow to update!
@manipulate for i=1:100
    VisJS.drawGraph(WheelGraph(i))
end


text: export A,B,   C


text:
using BinaryProvider

# Parse some basic command-line arguments
const verbose = "--verbose" in ARGS
const prefix = Prefix(get([a for a in ARGS if a != "--verbose"], 1, joinpath(@__DIR__, "usr")))
products = Product[
    LibraryProduct(prefix, String["libbid"], :libbid),
]

# Download binaries from hosted location
bin_prefix = "https://github.com/quinnj/DecFPBuilder/releases/download/v0.4"

# Listing of files generated by BinaryBuilder:
download_info = Dict(
    BinaryProvider.Linux(:i686, :glibc) => ("$bin_prefix/DecFP.i686-linux-gnu.tar.gz", "7a6f2460ffdf4f07a5b03246400a65f57f7d378490252af6b67f25da2df5b53d"),
    BinaryProvider.Windows(:i686) => ("$bin_prefix/DecFP.i686-w64-mingw32.tar.gz", "f7dbf57fce1450d1eaf891580491a1a982b8c3491b9e2aba7992d62da89d202d"),
    BinaryProvider.MacOS() => ("$bin_prefix/DecFP.x86_64-apple-darwin14.tar.gz", "b4fe9a370eb4b9bbe486c9a93a0aae0c1f1cca51e74d9b39bd7cdb0d9553e221"),
    BinaryProvider.Linux(:x86_64, :glibc) => ("$bin_prefix/DecFP.x86_64-linux-gnu.tar.gz", "ac336480e4c0200ff0d320e68f30feb40b2894105b5eaf49ae7b4ddf73078a67"),
    BinaryProvider.Windows(:x86_64) => ("$bin_prefix/DecFP.x86_64-w64-mingw32.tar.gz", "25eceff1f2b4d39042704f9a44a0295032d481d4240ec22e8000fc38ccf19be7"),
)

# First, check to see if we're all satisfied
if any(!satisfied(p; verbose=verbose) for p in products)
    if haskey(download_info, platform_key())
        # Download and install binaries
        url, tarball_hash = download_info[platform_key()]
        install(url, tarball_hash; prefix=prefix, force=true, verbose=verbose)
    else
        # If we don't have a BinaryProvider-compatible .tar.gz to download, complain.
        # Alternatively, you could attempt to install from a separate provider,
        # build from source or something more even more ambitious here.
        error("Your platform $(Sys.MACHINE) is not supported by this package!")
    end
end

# Write out a deps.jl file that will contain mappings for our products
write_deps_file(joinpath(@__DIR__, "deps.jl"), products)


text:
# WIP

@fastmath function ltc_gj_ca!(du,u,p,t, n_in, n_neurons, w_sens, w_syns, w_gj)
  dv = @view du[1:n_neurons]
  dCa = @view du[n_neurons:end]
  v = @view u[1:n_neurons]
  Ca = @view u[n_neurons:end]
  stim, C, G_L, E_L, G_sens, s_sens, h_sens, E_sens, G_syns, s_syns, h_syns, E_syns, G_gj = ltc_gj_unpack_p(n_in,n_neurons,p)

  I_senss = zeros(eltype(stim), n_in, n_neurons)
  I_synss = zeros(eltype(stim), n_neurons, n_neurons)
  I_gjs = zeros(eltype(stim), n_neurons, n_neurons)
  #I_tot = zeros(eltype(stim), n_neurons)
  @inbounds @simd for n in 1:n_neurons
    for s in 1:n_in
      I_senss[s,n] = w_sens[s,n] * G_sens[s,n] * mysigm((stim[s] - h_sens[s,n]) * s_sens[s,n]) * (v[n] - E_sens[s,n])
    end
    for s in 1:n_neurons
      I_synss[s,n] = w_syns[s,n] * G_syns[s,n] * mysigm((v[s] - h_syns[s,n]) * s_syns[s,n]) * (v[n] - E_syns[s,n])
      I_gjs[s,n] = w_gj[s,n] * G_gj[s,n] * (v[n] - v[s])
    end
  end
  I_sens = reshape(sum(I_senss, dims=1), :)
  I_syns = reshape(sum(I_synss, dims=1), :)
  I_gj = reshape(sum(I_gjs, dims=1), :)

  I_L = (G_L .* (v .- E_L))

  #I_Ca

  I_tot = I_L .+ I_sens .+ I_syns .+ I_gj #.+ I_Ca

  dv .= (-(I_tot) ./ C)
  #dCa .= (-I_Ca / 2*F*vol) - β*Ca

  nothing
end

function ltc_gj_ca_init_p(n_in, n_neurons; T=Float32)
  stim   = fill(T(0), n_in)
  #w_sens = zeros(T, n_in, n_neurons)
  #w_syns = zeros(T, n_neurons, n_neurons)
  C      = rand_uniform(T,  1.0,     1.002,   n_neurons)
  G_L    = rand_uniform(T,  0.0001,  0.1,     n_neurons)
  E_L    = rand_uniform(T, -0.3,     0.3,     n_neurons) #.± 0.1f0
  G_sens = rand_uniform(T,  0.0001,  0.1,     n_in, n_neurons)
  s_sens = rand_uniform(T,  3,       8,       n_in, n_neurons)
  h_sens = rand_uniform(T,  0.3,     0.8,     n_in, n_neurons)
  E_sens = rand_uniform(T, -0.9,     0.9,     n_in, n_neurons)
  G_syns = rand_uniform(T,  0.0001,  0.1,     n_neurons, n_neurons)
  s_syns = rand_uniform(T,  3,       8,       n_neurons, n_neurons)
  h_syns = rand_uniform(T,  0.3,     0.8,     n_neurons, n_neurons)
  E_syns = rand_uniform(T, -0.9,     0.9,     n_neurons, n_neurons)

  G_gj   = rand_uniform(T,  0.0001,  0.1,     n_neurons, n_neurons)

  v = rand_uniform(T, 0.01, 0.1, n_neurons)
  Ca = rand_uniform(T, 0.001, 0.01, n_neurons)

  #u0 = ComponentArray(v=v)
  #p = ComponentArrays(stim=stim, C=C, G_L=G_L, E_L=E_L,
  #                G_sens=G_sens, s_sens=s_sens, h_sens=h_sens, E_sens=E_sens,
  #                G_syns=G_syns, s_syns=s_syns, h_syns=h_syns, E_syns=E_syns)
  u0 = vcat(v, Ca)
  p = vcat(stim, #vec(w_sens), vec(w_syns),
           C, G_L, E_L,
           vec(G_sens), vec(s_sens), vec(h_sens), vec(E_sens),
           vec(G_syns), vec(s_syns), vec(h_syns), vec(E_syns),
           vec(G_gj),)
  return u0, p
end

function ltc_gj_ca_bounds(n_in, n_neurons; T=Float32)
  lb = vcat(
    [1.0 for _ in 1:n_neurons],               # C
    [0 for _ in 1:n_neurons],                 # G_L
    [-1 for _ in 1:n_neurons],                # E_L
    [0 for _ in 1:n_in*n_neurons],            # G_sens
    [2 for _ in 1:n_in*n_neurons],            # s_sens
    [-0.3 for _ in 1:n_in*n_neurons],          # h_sens
    [-1 for _ in 1:n_in*n_neurons],           # E_sens
    [0 for _ in 1:n_neurons*n_neurons],       # G_syns
    [2 for _ in 1:n_neurons*n_neurons],       # s_syns
    [-0.3 for _ in 1:n_neurons*n_neurons],     # h_syns
    [-1 for _ in 1:n_neurons*n_neurons],      # E_syns
    [0 for _ in 1:n_neurons*n_neurons],       # G_gj

    [-0.4 for _ in 1:n_neurons],              # v
    [0 for _ in 1:n_neurons]),             # Ca

  ub = vcat(
    [2.0 for _ in 1:n_neurons],
    [0.4 for _ in 1:n_neurons],
    [1 for _ in 1:n_neurons],
    [1 for _ in 1:n_in*n_neurons],
    [9 for _ in 1:n_in*n_neurons],
    [0.9 for _ in 1:n_in*n_neurons],
    [1 for _ in 1:n_in*n_neurons],
    [1 for _ in 1:n_neurons*n_neurons],
    [9 for _ in 1:n_neurons*n_neurons],
    [0.9 for _ in 1:n_neurons*n_neurons],
    [1 for _ in 1:n_neurons*n_neurons],
    [1 for _ in 1:n_neurons*n_neurons],       # G_gj

    [0.8 for _ in 1:n_neurons], # v
    [100 for _ in 1:n_neurons],) # Ca

  T.(lb), T.(ub)
end

function ltc_gj_ca_unpack_p(n_in,n_neurons,p)
  stiml = n_in
  Cl = n_neurons
  G_Ll = n_neurons
  E_Ll = n_neurons
  sensl = n_in*n_neurons
  synsl = n_neurons*n_neurons
  G_sensl = sensl
  s_sensl = sensl
  h_sensl = sensl
  E_sensl = sensl
  G_synsl = synsl
  s_synsl = synsl
  h_synsl = synsl
  E_synsl = synsl

  G_gjl = synsl

  s = 1
  stim = @view p[s:s+stiml-1]
  s += stiml
  C = @view p[s:s+Cl-1]
  s += Cl
  G_L = @view p[s:s+G_Ll-1]
  s += G_Ll
  E_L = @view p[s:s+E_Ll-1]
  s += E_Ll
  G_sens = @view p[reshape(s:s+G_sensl-1, n_in, n_neurons)]
  s += G_sensl
  s_sens = @view p[reshape(s:s+s_sensl-1, n_in, n_neurons)]
  s += s_sensl
  h_sens = @view p[reshape(s:s+h_sensl-1, n_in, n_neurons)]
  s += h_sensl
  E_sens = @view p[reshape(s:s+E_sensl-1, n_in, n_neurons)]
  s += E_sensl

  G_syns = @view p[reshape(s:s+G_synsl-1, n_neurons, n_neurons)]
  s += G_synsl
  s_syns = @view p[reshape(s:s+s_synsl-1, n_neurons, n_neurons)]
  s += s_synsl
  h_syns = @view p[reshape(s:s+h_synsl-1, n_neurons, n_neurons)]
  s += h_synsl
  E_syns = @view p[reshape(s:s+E_synsl-1, n_neurons, n_neurons)]
  s += E_synsl

  G_gj = @view p[reshape(s:s+G_gjl-1, n_neurons, n_neurons)]
  s += G_gjl

  return stim, C, G_L, E_L, G_sens, s_sens, h_sens, E_sens, G_syns, s_syns, h_syns, E_syns, G_gj
end


text:
"""
    PLists

This is a simple example of how to read an XML file.

```julia-repl
julia> doc = readxml("examples/note.xml");
julia> r = root(doc);
julia> ns = nodes(r);
```

You can use an array or dictionary interface to work with child
nodes or attribute nodes.

To read a NeXT or Apple style property list:

```julia-repl
julia> dict = readplist("test/example.plist")
```

This gives you a regular dictionary you can work with.
"""
module PLists

isalnum(c::Char) = isletter(c) || isnumeric(c)

include("plist_lexer.jl")
include("xml_lexer.jl")
include("parser.jl")
include("plist_parser.jl")
include("plist_writer.jl")
include("xml_dom.jl")
include("xml_parser.jl")
include("xml_show.jl")
include("xml_plist_parser.jl")

end


text:
struct VariableReference
    kind::Symbol
    value
end

mutable struct DebuggerState
    last_exception
    expr_splitter::Union{JuliaInterpreter.ExprSplitter,Nothing}
    frame
    not_yet_set_function_breakpoints::Set{Any}
    debug_mode::Symbol
    compile_mode
    sources::Dict{Int,String}
    next_source_id::Int
    varrefs::Vector{VariableReference}
    next_cmd::Channel{Any}
    not_yet_set_compiled_items::Vector{String}

    function DebuggerState()
        return new(nothing, nothing, nothing, Set{String}(), :unknown, JuliaInterpreter.finish_and_return!, Dict{Int,String}(), 1, VariableReference[], Channel{Any}(Inf), String[])
    end
end

is_toplevel_return(frame) = frame.framecode.scope isa Module && JuliaInterpreter.isexpr(JuliaInterpreter.pc_expr(frame), :return)

function attempt_to_set_f_breakpoints!(bps)
    for bp in bps
        @debug "Trying to set function breakpoint for '$(bp.name)'."
        try
            f = Core.eval(bp.mod, bp.name)

            signat = if bp.signature !== nothing
                Tuple{(Core.eval(Main, i) for i in bp.signature)...}
            else
                nothing
            end

            JuliaInterpreter.breakpoint(f, signat, bp.condition)
            delete!(bps, bp)

            @debug "Setting function breakpoint for '$(bp.name)' succeeded."
        catch err
            @debug "Setting function breakpoint for '$(bp.name)' failed."
        end
    end
end

function get_next_top_level_frame(state)
    state.expr_splitter === nothing && return nothing
    x = iterate(state.expr_splitter)
    x === nothing && return nothing

    (mod, ex), _ = x
    if Meta.isexpr(ex, :global, 1)
        # global assignment can be lowered, but global declaration can't,
        # let's just evaluate and iterate to next
        Core.eval(mod, ex)
        return get_next_top_level_frame(state)
    end
    return JuliaInterpreter.Frame(mod, ex)
end

function our_debug_command(cmd, state)
    while true
        @debug "Running a new frame." state.frame state.compile_mode

        ret = Base.invokelatest(JuliaInterpreter.debug_command, state.compile_mode, state.frame, cmd, true)

        state.not_yet_set_compiled_items = set_compiled_functions_modules!(state.not_yet_set_compiled_items)
        attempt_to_set_f_breakpoints!(state.not_yet_set_function_breakpoints)

        @debug "Finished running frame."

        if ret !== nothing && is_toplevel_return(ret[1])
            ret = nothing
        end

        if ret !== nothing
            state.frame = ret[1]
            return ret[2]
        end

        state.frame = get_next_top_level_frame(state)

        if state.frame === nothing
            return nothing
        end

        ret !== nothing && error("Invalid state.")

        if ret === nothing && (cmd == :n || cmd == :s || cmd == :finish || JuliaInterpreter.shouldbreak(state.frame, state.frame.pc))
            return state.frame.pc
        end
    end
end

function send_stopped_msg(conn, ret_val, state)
    if ret_val isa JuliaInterpreter.BreakpointRef
        if ret_val.err === nothing
            JSONRPC.send(conn, stopped_notification_type, StoppedEventArguments("breakpoint", missing, 1, missing, missing, missing))
        else
            state.last_exception = ret_val.err
            error_msg = Base.invokelatest(string, ret_val.err)
            JSONRPC.send(conn, stopped_notification_type, StoppedEventArguments("exception", missing, 1, missing, error_msg, missing))
        end
    elseif ret_val isa Number
        JSONRPC.send(conn, stopped_notification_type, StoppedEventArguments("step", missing, 1, missing, missing, missing))
    elseif ret_val === nothing
        JSONRPC.send(conn, stopped_notification_type, StoppedEventArguments("step", missing, 1, missing, missing, missing))
    end
end