text:
module DaemonTask

using ...SimpleTasks.Types

export Info, run, make_key

"""
    DaemonTask.Result

This object contains the outcome of performing the task.
"""
type Result
    success::Bool
    outputs::Array{AbstractString, 1}
end

# Test to see if the execute function exists for this type
function can_execute(task_type::Type)
    if !(task_type <: DaemonTaskDetails)
        return false
    end

    prepare_methods = methods(prepare, Type[task_type, DatasourceService])
    if length(prepare_methods) == 0
        return false
    end
    execute_methods = methods(execute, Type[task_type, DatasourceService])
    if length(execute_methods) == 0
        return false
    end
    finalize_methods = methods(finalize, Type[task_type, DatasourceService])
    if length(finalize_methods) == 0
        return false
    end
    for execute_method in execute_methods
        sig_types = execute_method.sig.types
        if  length(sig_types) > 0 && sig_types[1] == task_type
            return true
        end
    end
    return false
end

"""
    make_key(task::DaemonTaskDetails, input::AbstractString)

Get the fully qualified input key which is the input appended to the
base_directory of the task
"""
function make_key(task::DaemonTaskDetails, input::AbstractString)
    return "$(task.basic_info.base_directory)/$input"
end

"""
    run(task::DaemonTaskDetails, datasource::DatasourceService)

Run the current task. 1. Prepare, 2. Execute 3. Finalize
"""
function run(task::DaemonTaskDetails, datasource::DatasourceService)
    prepare(task, datasource)
    result = execute(task, datasource)
    finalize(task, datasource, result)
    return result;
end

"""
    prepare(task::DaemonTaskDetails, datasource::DatasourceService)

prepare what is needed for the task
"""
function prepare(task::DaemonTaskDetails, datasource::DatasourceService)
    error("Prepare is unimplemented for this task $task")
end

"""
    execute(task::DaemonTaskDetails, datasource::DatasourceService)

Executes the given task. Must be overriden for new tasks
"""
function execute(task::DaemonTaskDetails, datasource::DatasourceService)
    error("Execute is unimplemented for this task $task")
end

"""
    finalize(daemon::DaemonService, task::DaemonTaskDetails,

After task has completed, perform this action
"""
function finalize(task::DaemonTaskDetails, datasource::DatasourceService)
    error("finalize is unimplemented for this task $task")
end

end # module DaemonTask

ext: jl

max_forks_repo_name: seung-lab/SimpleTasks


text:
using Random
using Distributions
using BenchmarkTools
using RandomDataODE

n_samples = 10^4;
Tmax = 6.;

# generate initial condition data
ϵ = 0.2;
xdist = Uniform(0.5-ϵ, 0.5 + ϵ);
ydist = Uniform(2.0-ϵ, 2.0 + ϵ);
Random.seed!(100)
u0data = [[rand(xdist), rand(ydist)] for _ in 1:n_samples];

f = sol->sol(Tmax)[1]; # observable function

@btime samples = run_obs_ensemble(u0data, Tmax, f);


text:
module AutoFlux

using ..NaiveGAflux

export fit
export ImageClassification, ImageClassifier

"""
    fit(x, y; cb)
    fit(data::Tuple; cb, mdir)

Return a population of models fitted to the given data.

The type of model will depend on the shape of `x`.

The following model types are currently supported
- 4D data -> [`ImageClassifier`](@ref)

Keyword `cb` can be used to supply a callback function which will be called each generation with the current population as input.

Keyword `mdir` is a directory which will be searched for serialized state from which the optimization will resume operations. Particularly useful in combination with `cb=persist`.
"""
function fit(x, y; cb=identity, mdir=missing)
    if ndims(x) == 4
        outsize = ndims(y) == 1 ? length(unique(y)) : size(y, 1)
        return fit(ImageClassifier(insize=size(x), outsize, mdir=defaultdir(mdir, "ImageClassifier")), x, y; cb=identity)
    end
    error("No model for $(ndims(x))D data")
end
fit((x,y)::Tuple; cb=identity, mdir=missing) = fit(x, y; cb=identity, mdir=mdir)

modeldir(::Missing, subdir) = defaultdir(subdir)
modeldir(d, subdir) = d
defaultdir(subdir, basedir = NaiveGAflux.modeldir) = joinpath(basedir, subdir)

include("imageclassification/ImageClassification.jl")
using .ImageClassification

end


text:
#Steps to read and modify the SKyTEM headerfile
#module to read the hdr file and its column number
using DataFrames
using CSV
using DelimitedFiles
using Statistics
using PyPlot

function read_survey_files(hdrfile::String,fname_dat::String,fname_specs_halt::String;
    #fname_dat= "",
    #fname_specs_halt="",
    frame_height = "",
    frame_dz = "",
    frame_dx = "",
    frame_dy = "",
    LM_Z = "",
    HM_Z = "",
    LM_Ã = "",
    HM_Ã = "",
    relerror = false,
    units=1e-12,
    figsize = (9,7),
    makesounding = false,
    dotillsounding = nothing,
    startfrom = 1,
    skipevery = 1,
    multnoise = 0.03,
    X = "",
    Y = "",
    Z = "",
    fid = "",
    linenum = "")
    df = CSV.read(hdrfile, DataFrame; header=false);
    ColNo,ColName = eachcol(df)
    # get column numbers for each required parameters
    X = parse(Int64,(df[:,1][findfirst(==(filter(in(Set(("Easting", "easting"))), ColName)[1]),ColName)]))
    Y = parse(Int64,(df[:,1][findfirst(==(filter(in(Set(("Northing", "northing"))), ColName)[1]),ColName)]))
    Z = parse(Int64,(df[:,1][findfirst(==(filter(in(Set(("ground_elevation_lidar_aem_merged","GPS_Alt","gps_height"))),ColName)[1]),ColName)]))
    fid = parse(Int64,(df[:,1][findfirst(==(filter(in(Set(("Fiducial","fiducial"))),ColName)[1]),ColName)]))
    linenum = parse(Int64,(df[:,1][findfirst(==(filter(in(Set(("line", "Line", "Line number","line number","Line_Number","Line number"))),ColName)[1]),ColName)]))
    frame_height = parse(Int64,(df[:,1][findfirst(==(filter(in(Set(("frame_height","Frame_height"))),ColName)[1]),ColName)]))
    frame_dx = parse(Int64,(df[:,1][findfirst(==(filter(in(Set(("frame_dx","Frame_dx"))),ColName)[1]),ColName)]))
    frame_dy = parse(Int64,(df[:,1][findfirst(==(filter(in(Set(("frame_dy","Frame_dy"))),ColName)[1]),ColName)]))
    frame_dz = parse(Int64,(df[:,1][findfirst(==(filter(in(Set(("frame_dz","Frame_dz"))),ColName)[1]),ColName)]))
    LM_Z = parse.(Int64,split(df[:,1][findfirst(==(filter(in(Set(("LM_Z","Z_LM"))),ColName)[1]),ColName)],"-"))
    HM_Z = parse.(Int64,split(df[:,1][findfirst(==(filter(in(Set(("HM_Z","Z_HM"))),ColName)[1]),ColName)],"-"))

    # println("X = $X", "\n","Y = $Y", "\n" ,"Z = $Z", "\n" ,
    # "fid = $fid","\n" ,"linenum = $linenum","\n",
    # "frame_height = $frame_height","\n",
    # "frame_dx = $frame_dx","\n",
    # "frame_dy = $frame_dy","\n",
    # "frame_dz = $frame_dz","\n",
    # "LM_Z = $LM_Z","\n",
    # "HM_Z = $HM_Z","\n",
    # "relerror = false","\n",
    # "units = 1e-12")

    #call read_survey_file
    @assert frame_height > 0
    @assert frame_dz > 0
    @assert frame_dx > 0
    @assert frame_dy > 0
    @assert all(LM_Z .> 0)
    @assert all(HM_Z .> 0)
    if relerror
        @assert all(LM_Ã .> 0)
        @assert all(HM_Ã .> 0)
    end
    @assert X > 0
    @assert Y > 0
    @assert Z > 0
    @assert linenum > 0
    @assert fid > 0
    @info "reading $fname_dat"
    if dotillsounding!= nothing
        soundings = readdlm(fname_dat)[startfrom:skipevery:dotillsounding,:]
    else
        soundings = readdlm(fname_dat)[startfrom:skipevery:end,:]
    end
    easting = soundings[:,X]
    northing = soundings[:,Y]
    topo = soundings[:,Z]
    fiducial = soundings[:,fid]
    whichline = soundings[:,linenum]
    d_LM = soundings[:,LM_Z[1]:LM_Z[2]]
    d_HM = soundings[:,HM_Z[1]:HM_Z[2]]
    if relerror
        Ã_LM = soundings[:,LM_Ã[1]:LM_Ã[2]]
        Ã_HM = soundings[:,HM_Ã[1]:HM_Ã[2]]
    end
    zTx = soundings[:,frame_height]
    zRx = -(zTx + soundings[:,frame_dz])
    zTx = -zTx
    rRx = sqrt.(soundings[:,frame_dx].^2 + soundings[:,frame_dy].^2)

    @info "reading $fname_specs_halt"
    include(fname_specs_halt)
    @assert size(d_LM, 2) == length(LM_times)
    @assert size(d_HM, 2) == length(HM_times)
    if !relerror
        @assert size(d_LM, 2) == length(LM_noise)
        @assert size(d_HM, 2) == length(HM_noise)
        LM_noise[:] .*= units
        HM_noise[:] .*= units
    else
        Ã_LM[:] .*= units
        Ã_HM[:] .*= units
    end
    d_LM[:]     .*= units
    d_HM[:]     .*= units
    f = figure(figsize=figsize)
    ax = Array{Any, 1}(undef, 4)
    ax[1] = subplot(2,2,1)
    nsoundings = size(soundings, 1)
    plot_dLM = permutedims(d_LM)
    plot_dLM[plot_dLM .<0] .= NaN
    pcolormesh(1:nsoundings, LM_times, log10.(plot_dLM), shading="nearest")
    xlabel("sounding #")
    cbLM = colorbar()
    cbLM.set_label("log d_LM")
    ylabel("LM time s")
    axLM = ax[1].twiny()
    axLM.semilogy(LM_noise, LM_times)
    axLM.set_xlabel("high alt noise")
    ax[2] = subplot(2,2,3,sharex=ax[1], sharey=ax[1])
    plot_dHM = permutedims(d_HM)
    plot_dHM[plot_dHM .<0] .= NaN
    pcolormesh(1:nsoundings, HM_times, log10.(plot_dHM), shading="nearest")
    xlabel("sounding #")
    cbHM = colorbar()
    cbHM.set_label("log d_HM")
    ylabel("HM time s")
    ax[2].invert_yaxis()
    axHM = ax[2].twiny()
    axHM.semilogy(HM_noise, HM_times)
    axHM.set_xlabel("high alt noise")
    ax[3] = subplot(2,2,2, sharex=ax[1])
    plot(1:nsoundings, zRx, label="Rx")
    plot(1:nsoundings, zTx, label="Tx")
    legend()
    xlabel("sounding #")
    ylabel("height m")
    ax[3].invert_yaxis()
    ax[4] = subplot(2,2,4, sharex=ax[1])
    ax[4].plot(1:nsoundings, rRx)
    xlabel("sounding #")
    ylabel("rRx m")
    plt.tight_layout()
    if makesounding
        s_array = Array{SkyTEMsoundingData, 1}(undef, nsoundings)
        for is in 1:nsoundings
            l, f = Int(whichline[is]), fiducial[is]
            @info "read $is out of $nsoundings"
            dlow, dhigh = vec(d_LM[is,:]), vec(d_HM[is,:])
            if !relerror
                ÃLM = sqrt.((multnoise*dlow).^2 + LM_noise.^2)
                ÃHM = sqrt.((multnoise*dhigh).^2 + HM_noise.^2)
            else
                ÃLM, ÃHM = vec(Ã_LM[is,:]), vec(Ã_HM[is,:])
            end
            s_array[is] = SkyTEMsoundingData(rRx=rRx[is], zRxLM=zRx[is], zTxLM=zTx[is],
                zRxHM=zRx[is], zTxHM=zTx[is], rTx=rTx, lowpassfcs=lowpassfcs,
                LM_times=LM_times, LM_ramp=LM_ramp,
                HM_times=HM_times, HM_ramp=HM_ramp,
                LM_noise=ÃLM, HM_noise=ÃHM, LM_data=dlow, HM_data=dhigh,
                sounding_string="sounding_$(l)_$f",
                X=easting[is], Y=northing[is], Z=topo[is], fid=f,
                linenum=l)
        end
        return s_array
    end
end

read_survey_files("em_with_noise_calibration_lines.hdr",
"/g/data/z67/nfm547/transD_GP/examples/SkyTEM1D/Menindee/coincidewithtempest.dat",
"/g/data/z67/nfm547/transD_GP/examples/SkyTEM1D/Menindee/electronics_halt.jl")


text:
abstract type RackettTranslationModel <: TranslationModel end

struct RackettTranslationParam <: EoSParam
    Vc::SingleParam{Float64}
end

@newmodelsimple RackettTranslation RackettTranslationModel RackettTranslationParam

export RackettTranslation
function RackettTranslation(components::Vector{String}; userlocations::Vector{String}=String[], verbose::Bool=false)
    params = getparams(components, ["properties/critical.csv"]; userlocations=userlocations, verbose=verbose)
    Vc = params["vc"]
    packagedparams = RackettTranslationParam(Vc)
    model = RackettTranslation(packagedparams, verbose=verbose)
    return model
end

function translation(model::CubicModel,V,T,z,translation_model::RackettTranslation)
    Tc = model.params.Tc.values
    Pc = model.params.Pc.values
    Vc = translation_model.params.Vc.values
    return @. 0.40768*R̄*Tc/Pc*(0.29441-Pc*Vc/(R̄*Tc))
end


text:
@testset "orc" begin

let ctx = Context()
    tm  = JITTargetMachine()
    orc = OrcJIT(tm)

    known_functions = Dict{String, OrcTargetAddress}()
    fnames = Dict{String, Int}()
    function lookup(name, ctx)
        name = unsafe_string(name)
        try
            if !haskey(fnames, name)
                fnames[name] = 0
            end
            fnames[name] += 1

            return get(known_functions, name, OrcTargetAddress(C_NULL)).ptr
        catch ex
            @error "Exception during lookup" exception=(ex, catch_backtrace())
            return UInt64(0)
        end
    end

    mod = LLVM.Module("jit", ctx)
    T_Int32 = LLVM.Int32Type(ctx)
    ft = LLVM.FunctionType(T_Int32, [T_Int32, T_Int32])
    fn = LLVM.Function(mod, "mysum", ft)
    linkage!(fn, LLVM.API.LLVMExternalLinkage)

    fname = mangle(orc, "wrapper")
    wrapper = LLVM.Function(mod, fname, ft)
    # generate IR
    Builder(ctx) do builder
        entry = BasicBlock(wrapper, "entry", ctx)
        position!(builder, entry)

        tmp = call!(builder, fn, [parameters(wrapper)...])
        ret!(builder, tmp)
    end

    triple!(mod, triple(tm))
    ModulePassManager() do pm
        add_library_info!(pm, triple(mod))
        add_transform_info!(pm, tm)
        run!(pm, mod)
    end
    verify(mod)

    mysum_name = mangle(orc, "mysum")
    known_functions[mysum_name] = OrcTargetAddress(@cfunction(+, Int32, (Int32, Int32)))

    f_lookup = @cfunction($lookup, UInt64, (Cstring, Ptr{Cvoid}))
    GC.@preserve f_lookup begin
        orc_mod = compile!(orc, mod, f_lookup, lazy=true) # will capture f_lookup

        addr = address(orc, fname)
        @test errormsg(orc) == ""
        addr2 = addressin(orc, orc_mod, fname)

        @test addr == addr2
        @test addr.ptr != 0
        @test !haskey(fnames, mysum_name)

        r = ccall(pointer(addr), Int32, (Int32, Int32), 1, 2) # uses f_lookup
        @test r == 3
    end

    @test haskey(fnames, mysum_name)
    @test fnames[mysum_name] == 1

    empty!(fnames)
    delete!(orc, orc_mod)
    dispose(orc)
end

let ctx = Context()
    tm = JITTargetMachine()
    orc = OrcJIT(tm)
    sym = mangle(orc, "SomeFunction")

    mod = LLVM.Module("jit", ctx)
    ft = LLVM.FunctionType(LLVM.VoidType(ctx))
    fn = LLVM.Function(mod, sym, ft)

    Builder(ctx) do builder
        entry = BasicBlock(fn, "entry")
        position!(builder, entry)
        ret!(builder)
    end
    verify(mod)

    obj = emit(tm, mod, LLVM.API.LLVMObjectFile)

    orc_m = add!(orc, MemoryBuffer(obj))
    addr = address(orc, sym)

    @test addr.ptr != 0
    delete!(orc, orc_m)

    toggle = Ref{Bool}(false)
    on()  = (toggle[] = true; nothing)
    off() = (toggle[] = false; nothing)

    # Note that `CFunction` objects can be GC'd (???)
    # and we capture them below.
    func_on = @cfunction($on, Cvoid, ())
    GC.@preserve func_on begin
        ptr = Base.unsafe_convert(Ptr{Cvoid}, func_on)

        create_stub!(orc, "mystub", OrcTargetAddress(ptr))
        addr = address(orc, "mystub")

        @test addr.ptr != 0
        @test toggle[] == false

        ccall(pointer(addr), Cvoid, ())
        @test toggle[] == true
    end

    func_off = @cfunction($off, Cvoid, ())
    GC.@preserve func_off begin
        ptr = Base.unsafe_convert(Ptr{Cvoid}, func_off)

        set_stub!(orc, "mystub", OrcTargetAddress(ptr))

        @test addr == address(orc, "mystub")
        @test toggle[] == true

        ccall(pointer(addr), Cvoid, ())
        @test toggle[] == false
    end

    dispose(orc)
end

# TODO:
# Test for `callback!`, currently unsure how to trigger that.

end


text:
  using PDBTools  #, ComplexMixtures
  function autocorr(trajectory::Trajectory)

    # number of atoms
    nprot = length(trajectory.x_solute)
    nsvt  = length(trajectory.x_solvent)
    nframes = trajectory.nframes

    # vector with time - ns
    delta = 0.01
    time = zeros(nframes)
    t1 = 0.
    for i in 1:nframes
      t1 = t1 + delta
      time[i] = t1
    end

    # matrix of distances, absorptions, number of events and sample space
    Dist      = Array{Float64}(undef,nframes,nsvt)
    domain    = zeros(Float64, nframes, nsvt)
    evals     = zeros(Float64, nframes, nsvt)
    sp        = zeros(Float64, nframes, nsvt)

    for iframe in 1:nframes

      # reading coordinates of next frame
      nextframe!(trajectory)

      # variables to compute the autocorrelation function
      solute    = trajectory.solute
      solvent   = trajectory.solvent
      x_solute  = trajectory.x_solute
      x_solvent = trajectory.x_solvent

      # box
      sides = getsides(trajectory,iframe)

      for i in 1:nsvt            # loop though solvent atoms
        # distance
        dist = +Inf
        dims = x_solvent[i]

        for p in 1:nprot

          x_this_solute = viewmol(i,x_solvent,solvent)
          pt = x_solute[p]
          D = distance(dims,pt,sides)                  # distance calculation

          # to calculate the protein atom that are at the closest distance to the sovlent
          if dist > D
            dist = D
          end

        end
        Dist[iframe,i] = dist
      end
    end

    closetraj(trajectory)
    return Dist, domain, evals, sp, nprot, nsvt, nframes, time

  end


text:
"""
    MPIData(data; kwargs...)

Wrapper used for signaling that this data is part of a distributed data structure and
communication should happen with MPI. `MPIData` can generally be used where an
[`AbstractDVec`](@ref) would be used otherwise. Unlike [`AbstractDVec`](@ref)s, `MPIData`
does not support indexing, or iteration over `keys`, `values`, and `pairs`.

Keyword arguments:
* `setup = mpi_point_to_point` - controls the communication stratgy
  * [`mpi_one_sided`](@ref) uses one-sided communication with remote memory access (RMA), sets [`MPIOneSided`](@ref) strategy.
  * [`mpi_point_to_point`](@ref) uses [`MPIPointTOPoint`](@ref) strategy.
  * [`mpi_all_to_all`](@ref) uses [`MPIAllToAll`](@ref) strategy.
  * [`mpi_no_exchange`](@ref) sets [`MPINoWalkerExchange`](@ref) strategy. Experimental. Use with caution!
* `comm = mpi_comm()`
* `root = mpi_root`
* The rest of the keyword arguments are passed to `setup`.
"""
struct MPIData{D,S}
    data::D # local data, e.g. a DVec
    comm::MPI.Comm
    root::Int32 # rank of root process
    isroot::Bool # true if running on root process
    s::S # type (struct) with further details needed for communication

    function MPIData(data::D, comm, root, s::S) where {D, S<:DistributeStrategy}
        return new{D,S}(data, comm, root, s.id == root, s)
    end
end
# convenient constructor with setup function
function MPIData(data; setup=mpi_point_to_point, comm=mpi_comm(), root=mpi_root, kwargs...)
    return setup(data, comm, root; kwargs...)
end

Base.eltype(md::MPIData) = eltype(md.data)
Base.valtype(md::MPIData) = valtype(md.data)
Base.keytype(md::MPIData) = keytype(md.data)
Rimu.localpart(md::MPIData) = md.data
Rimu.StochasticStyle(d::MPIData) = Rimu.StochasticStyle(d.data)

function Base.summary(io::IO, md::MPIData)
    data = nameof(typeof(md.data))
    len = length(md)
    style = StochasticStyle(md)
    strat = nameof(typeof(md.s))
    print(io, "MPIData($data) with $len entries, style = $style, strategy = $strat")
end
function Base.show(io::IO, md::MPIData)
    summary(io, md)
    limit, _ = displaysize()
    for (i, p) in enumerate(pairs(localpart(md)))
        if length(md) > i > limit - 4
            print(io, "\n  ⋮   => ⋮")
            break
        else
            print(io, "\n  ", p)
        end
    end
end

###
### Iterators
###
"""
    MPIDataIterator{I,M<:MPIData}

Iterator over `keys`, `values`, or `pairs` of a `dv::MPIData`. Unlike its name would
suggest, it does not actually support iteration. To perform computations with it, use
`mapreduce`, or its derivatives (`sum`, `prod`, `reduce`...), which will perform the
reduction accross MPI ranks.
"""
struct MPIDataIterator{I,M<:MPIData}
    iter::I
    data::M
end

function Base.iterate(it::MPIDataIterator, args...)
    error(
        "iterating over `::MPIData` is not supported. ",
        "Use `localpart` to iterate over the local part of the vector or `mapreduce` to ",
        "perform a reduction accross ranks",
    )
end

function Base.mapreduce(f, op, it::MPIDataIterator; kwargs...)
    res = mapreduce(f, op, it.iter; kwargs...)
    return MPI.Allreduce(res, op, it.data.comm)
end

Base.IteratorSize(::MPIDataIterator) = Base.SizeUnknown()
Base.pairs(data::MPIData) = MPIDataIterator(pairs(localpart(data)), data)
Base.keys(data::MPIData) = MPIDataIterator(keys(localpart(data)), data)
Base.values(data::MPIData) = MPIDataIterator(values(localpart(data)), data)

Rimu.localpart(it::MPIDataIterator) = it.iter

"""
    length(md::MPIData)

Compute the length of the distributed data on every MPI rank with
`MPI.Allreduce`. MPI syncronizing.
"""
Base.length(md::MPIData) = MPI.Allreduce(length(md.data), +, md.comm)

"""
    norm(md::MPIData, p=2)

Compute the norm of the distributed data on every MPI rank with `MPI.Allreduce`.
MPI syncronizing.
"""
function LinearAlgebra.norm(md::MPIData, p::Real=2)
    if p === 2
        return sqrt(sum(abs2, values(md)))
    elseif p === 1
        return float(sum(abs, values(md)))
    elseif p === Inf
        return float(mapreduce(abs, max, values(md); init=real(zero(valtype(md)))))
    else
        error("$p-norm of MPIData is not implemented.")
    end
end

"""
    walkernumber(md::MPIData)

Compute the walkernumber of the distributed data on every MPI rank with `MPI.Allreduce`.
MPI syncronizing.
"""
function Rimu.DictVectors.walkernumber(md::MPIData)
    return MPI.Allreduce(walkernumber(md.data), +, md.comm)
end

"""
    mpi_synchronize!(md::MPIData)

Synchronize `md`, ensuring its contents are distributed among ranks correctly.
"""
function mpi_synchronize!(md::MPIData)
    P = eltype(md)
    myrank = mpi_rank(md.comm)
    buffers = Vector{P}[P[] for _ in 1:mpi_size(md.comm)]
    for (add, val) in pairs(localpart(md))
        tr = targetrank(add, mpi_size(md.comm))
        if tr ≠ myrank
            push!(buffers[tr + 1], add => val)
            localpart(md)[add] = zero(valtype(md))
        end
    end
    mpi_communicate_buffers!(localpart(md), buffers, md.comm)
    return md
end

"""
    mpi_communicate_buffers!(target::AbstractDVec{K,V}, buffers::Vector{<:Vector{V}})

Use MPI to communicate the contents of `buffers` and sort them into `target`. The length
of `buffers` should be equal to [`mpi_size`](@ref).
"""
function mpi_communicate_buffers!(target, buffers, comm)
    myrank = mpi_rank(comm)
    recbuf = buffers[myrank + 1]
    datatype = MPI.Datatype(eltype(target))
    # Receive from lower ranks.
    for id in 0:(myrank - 1)
        resize!(recbuf, MPI.Get_count(MPI.Probe(id, 0, comm), datatype))
        MPI.Recv!(recbuf, id, 0, comm)
        for (add, value) in recbuf
            target[add] += value
        end
    end
    # Perform sends.
    for id in 0:(mpi_size(comm) - 1)
        id == myrank && continue
        MPI.Send(buffers[id + 1], id, 0, comm)
    end
    # Receive from higher ranks.
    for id in (myrank + 1):(mpi_size(comm) - 1)
        resize!(recbuf, MPI.Get_count(MPI.Probe(id, 0, comm), datatype))
        MPI.Recv!(recbuf, id, 0, comm)
        for (add, value) in recbuf
            target[add] += value
        end
    end
    return target
end

"""
    *(lop::AbstractHamiltonian, md::MPIData)
Allocating "Matrix"-"vector" multiplication with MPI-distributed "vector" `md`. The result is similar to
[`localpart(md)`](@ref) with all content having been communicated to the correct [`targetrank`](@ref).
MPI communicating.

See [`MPIData`](@ref).
"""
function Base.:*(lop, md::MPIData)
    T = promote_type(eltype(lop),valtype(md))
    P = Pair{keytype(md),T}
    buffers = Vector{P}[P[] for _ in 1:mpi_size(md.comm)]
    myrank = mpi_rank()

    result = similar(localpart(md), T)

    # Sort values into buffers and communicate.
    for (key, val) in pairs(localpart(md))
        result[key] += diagonal_element(lop, key)*val
        for (add, elem) in offdiagonals(lop, key)
            tr = targetrank(add, mpi_size(md.comm))
            if tr == myrank
                result[add] += elem * val
            else
                push!(buffers[tr + 1], add => elem * val)
            end
        end
    end
    mpi_communicate_buffers!(result, buffers, md.comm)

    return result
end

# Note: the following methods assume MPIDatas are distributed correctly.
function LinearAlgebra.dot(x, md::MPIData)
    return MPI.Allreduce(localpart(x)⋅localpart(md), +, md.comm)
end
function LinearAlgebra.dot(md::MPIData, x)
    return MPI.Allreduce(localpart(md)⋅localpart(x), +, md.comm)
end
function LinearAlgebra.dot(md_left::MPIData, md_right::MPIData)
    return MPI.Allreduce(localpart(md_left)⋅localpart(md_right), +, md_left.comm)
end

# Note: the following two methods work with x::DVec and assume `x` is the same on all ranks.
function LinearAlgebra.dot(x, lop, md::MPIData)
    return MPI.Allreduce(dot(x, lop, localpart(md)), +, md.comm)
end
function LinearAlgebra.dot(md::MPIData, lop, x)
    return MPI.Allreduce(dot(localpart(md), lop, x), +, md.comm)
end

"""
    copy_to_local(md::MPIData)

Collect all pairs in `md` from all ranks and store them in a local `AbstractDVec`.
"""
copy_to_local(md::MPIData) = copy_to_local!(similar(localpart(md)), md)

"""
    copy_to_local!(target, md::MPIData)

Collect all pairs in `md` from all ranks and store them in `target`.
In-place version of [`copy_to_local`](@ref).
"""
function copy_to_local!(target, md::MPIData)
    copy!(target, localpart(md))
    comm = md.comm
    myrank = mpi_rank(comm)
    datatype = MPI.Datatype(eltype(md))

    # Store all pairs to a buffer.
    sendbuf = collect(pairs(localpart(md)))
    recbuf = eltype(md)[]

    # Receive from lower ranks.
    for id in 0:(myrank - 1)
        resize!(recbuf, MPI.Get_count(MPI.Probe(id, 0, comm), datatype))
        MPI.Recv!(recbuf, id, 0, comm)
        for (add, value) in recbuf
            target[add] += value
        end
    end
    # Perform sends.
    for id in 0:(mpi_size(comm) - 1)
        id == myrank && continue
        MPI.Send(sendbuf, id, 0, comm)
    end
    # Receive from higher ranks.
    for id in (myrank + 1):(mpi_size(comm) - 1)
        resize!(recbuf, MPI.Get_count(MPI.Probe(id, 0, comm), datatype))
        MPI.Recv!(recbuf, id, 0, comm)
        for (add, value) in recbuf
            target[add] += value
        end
    end
    MPI.Barrier(comm)
    return target
end

function LinearAlgebra.dot(md_left::MPIData, lop, md_right::MPIData)
    # Idea: lop * md_right can be huge. It might be better to just collect the full left
    # vector and do the multiplication locally.
    left = copy_to_local(md_left)
    return dot(left, lop, md_right)
end

function Rimu.freeze(md::MPIData)
    mpi_synchronize!(md)
    return freeze(localpart(md))
end

function Rimu.all_overlaps(operators::Tuple, vecs::NTuple{N,MPIData}) where {N}
    local_vec_i = similar(localpart(vecs[1]))
    T = promote_type((valtype(v) for v in vecs)..., eltype.(operators)...)
    names = String[]
    values = T[]
    for i in 1:N, j in i+1:N
        push!(names, "c$(i)_dot_c$(j)")
        push!(values, dot(vecs[i], vecs[j]))
        copy_to_local!(local_vec_i, vecs[i])
        for (k, op) in enumerate(operators)
            push!(names, "c$(i)_Op$(k)_c$(j)")
            push!(values, dot(local_vec_i, op, vecs[j]))
        end
    end

    num_reports = (N * (N - 1) ÷ 2) * (length(operators) + 1)
    return SVector{num_reports,String}(names).data, SVector{num_reports,T}(values).data
end


text:
export TR1D_N, TR1D_s, TR1D_s2, TR1D_C2

""" `TR1D_N(ϕ, a, b, stp, α, β, η₁, η₂, red. aug, M)`
   - ϕ is a linemodel
   - kwargs from `bracket`
      - minimization in the interval [a,b]
      - solution to satisfy   α ≤ ϕ'(t\\*) ≤ β and ϕ(t\\*)<ϕ(a)
      - stp : stopping object; atol and rtol are bypassed and the conditions above apply

`TR1D_N` is a scalar trust region implementation using the Taylor quadratic model (Newton's method)
"""
function TR1D_N(ϕ; kwargs...)
    return TR1D(ϕ; M = Taylor2{Float64}(), kwargs...)
end

""" `TR1D_N(ϕ, a, b, stp, α, β, η₁, η₂, red. aug, M)`
   - ϕ is a linemodel
   - kwargs from `bracket`
      - minimization in the interval [a,b]
      - solution to satisfy   α ≤ ϕ'(t\\*) ≤ β and ϕ(t\\*)<ϕ(a)
      - stp : stopping object; atol and rtol are bypassed and the conditions above apply

`TR1D_s` is a scalar trust region implementation using the two-points secant  quadratic model
"""
function TR1D_s(ϕ; kwargs...)
    return TR1D(ϕ; M = secant{Float64}(), kwargs...)
end

""" `TR1D_C2(ϕ, a, b, stp, α, β, η₁, η₂, red. aug, M)`
   - ϕ is a linemodel
   - kwargs from `bracket`
      - minimization in the interval [a,b]
      - solution to satisfy   α ≤ ϕ'(t\\*) ≤ β and ϕ(t\\*)<ϕ(a)
      - stp : stopping object; atol and rtol are bypassed and the conditions above apply

`TR1D_C2` is a scalar trust region implementation using the two-points cubic model. Convergence not proved.
"""
function TR1D_C2(ϕ; kwargs...)
    return TR1D(ϕ; M = TwoPoints3{Float64}(), kwargs...)
end

""" `TR1D_C2(ϕ, a, b, stp, α, β, η₁, η₂, red. aug, M)`
   - ϕ is a linemodel
   - kwargs from `bracket`
      - minimization in the interval [a,b]
      - solution to satisfy   α ≤ ϕ'(t\\*) ≤ β and ϕ(t\\*)<ϕ(a)
      - stp : stopping object; atol and rtol are bypassed and the conditions above apply

`TR1D_s2` is a scalar trust region implementation using the two-points cubic model approximated by 2 secant sub-iterations. Convergence not proved.
"""
function TR1D_s2(ϕ; kwargs...)
    return TR1D(ϕ; M = secant3{Float64}(), kwargs...)
end


text:
#==
# Some basic numerical integration routines
> Code taken from: https://github.com/cmcbride/julia/blob/d54095465da8ad06d1d473f6709df0d599f02c7d/base/integrate.jl
Issues to try and fix:
- This assumes that eltype exists, which may not be true for iterables. It also assumes that zero exists, which is not true if e.g. eltype is an array.
- The right thing to do is to consider three cases separately:
    - zero or one elements: return sum(y)*(1.0*dx). The sum function will do the right thing,
      but you need the 1.0*dx to ensure that the type is the same as below. 2 elements: (y[1]+y[2])*(dx*0.5)
    - 3+ elements: iterate, but initialize the summand to y[1]*0.5 + y[2]
      However, you have to use start/next/done rather than length(y) or y[i] to get the length and elements, respectively.

# Future ideas
Port over as much as I can from scipy.integrate. It would be nice to have feature parity
with that package.
==#

# simple integration routines mostly for tabulated data

function trapz(y)
    trapz(y, one(eltype(y)))
end

# integration for uniform points (in x)
function trapz(y, dx::Number)
    r = zero(zero(eltype(y)) * zero(dx))
    r += (y[1] + y[end])/2
    r += sum(y[2:end-1])
    r * dx
end

# integration for non-uniform points (in x)
function trapz(y, x)
    n = length(y)
    if n != length(x)
        throw(ArgumentError("Input x,y must be of same length"))
    end
    r = zero(zero(eltype(x))*zero(eltype(y)))
    for i in 2:n
        r += (x[i] - x[i-1]) * (y[i] + y[i-1])
    end
    r/2.0 #TODO: this forces the output to be Float64
end

function simps(y)
    simps(y, one(eltype(y)))
end

function simps(y, dx::Number)
    n = length(y)
    if iseven(n)
        throw(ArgumentError("Simpson rule requires ODD length input (EVEN number of intervals)"))
    end
    r = zero(zero(eltype(y))*zero(dx))
    for i in 2:2:n-1
        r += y[i]
    end
    r *= 2.0 #TODO: this forces the output to be Float64
    for i in 3:2:n-2
        r += y[i]
    end
    r *= 2.0 #TODO: this forces the output to be Float64
    for i in (1,n)
        r += y[i]
    end
    r * dx / 3.0 #TODO: this forces the output to be Float64
end

function simps(y, x)
    n = length(y)
    if n != length(x)
        throw(ArgumentError("Input x,y must be of same length"))
    end
    if iseven(n)
        throw(ArgumentError("Simpson rule requires ODD length input (EVEN number of intervals)"))
    end

    # this is a quick generalization of the simpson's rule for
    # arbitrary separations in x
    r = zero(zero(eltype(x))*zero(eltype(y)))
    for i in 3:2:n
        d1 = x[i-1] - x[i-2]
        d2 = x[i] - x[i-1]
        h = x[i] - x[i-2]
        r += h * (
                y[i-2] * (2 - d2 / d1) +
                y[i-1] * h * h / (d1 * d2) +
                y[i] * (2 - d1 / d2)
            )
    end
    r / 6.0 #TODO: this forces the output to be Float64
end

max_stars_repo_name: gabrielgellner/Quadrature

max_issues_repo_name: gabrielgellner/Quadrature

max_forks_repo_name: gabrielgellner/Quadrature


text:
m = Model()

@variable(m, x ≥ 0)
@variable(m, y)

xu = UnitJuMP.UnitVariableRef(x, u"m/s")
yu = UnitJuMP.UnitVariableRef(y, u"km/hr")

@test _MA.@rewrite(xu) == UnitJuMP.UnitAffExpr(1x, u"m/s")
@test _MA.@rewrite(-xu) == UnitJuMP.UnitAffExpr(-x, u"m/s")
@test _MA.@rewrite(5xu) == UnitJuMP.UnitAffExpr(5x, u"m/s")
@test _MA.@rewrite(xu/5) == UnitJuMP.UnitAffExpr(0.2x, u"m/s")

@test _MA.@rewrite(xu + yu) == UnitJuMP.UnitAffExpr(x + y / 3.6, u"m/s")
@test _MA.@rewrite(xu - yu) == UnitJuMP.UnitAffExpr(x - y / 3.6, u"m/s")
@test _MA.@rewrite(-xu + yu) == UnitJuMP.UnitAffExpr(-x + y / 3.6, u"m/s")

speed = 10u"m/s"
@test _MA.@rewrite(speed*x) == UnitJuMP.UnitAffExpr(10x, u"m/s")
@test _MA.@rewrite(x*speed) == UnitJuMP.UnitAffExpr(10x, u"m/s")
@test _MA.@rewrite(speed*xu) == UnitJuMP.UnitAffExpr(10x, u"m^2/s^2")
@test _MA.@rewrite(xu*speed) == UnitJuMP.UnitAffExpr(10x, u"m^2/s^2")
@test _MA.@rewrite(x/speed) == UnitJuMP.UnitAffExpr(0.1x, u"s/m")


text:
using Abaco
using Test

interval=10
ts = nowts()

val1 = 100
val2 = 101
val3 = 102

expected_triggers = 3
actual_triggers = 0

#function onresult(ts, ne, name, value, inputs)
#    global actual_triggers
#    @debug "timestamp [$ts]: ne: [$ne] $name = $value"
#    if name === "offset_footprint"
#        @test value == x + val1 + val2 + val3
#    end
#    actual_triggers += 1
#end

function onresult(timestamp, ne, name, value::Abaco.PValue, inputs)
    global actual_triggers
    @debug "progressive timestamp [$ts]: ne: [$ne] $name = $value"
    if value.contribs == 1
        @test ne == "italy"
        @test timestamp == span(ts, interval)
        @test value.value == val1
    elseif value.contribs == 2
        @test value.value == val1 + val2
    elseif value.contribs == 3
        @test value.value == val1 + val2 + val3
    else
        # never here
        @test 1 == 0
    end
    actual_triggers += 1
end

abaco = abaco_init(onresult, interval=interval)

@debug "starting demo with time span: $interval, values timestamp: $ts"

setup_settings(abaco, "state", oncomplete=onresult)

formula(abaco, "total", "sum(region.city.footprint)", "state")

italy = node(abaco, "italy", "state")
veneto = node(abaco, italy, "veneto", "region")
trentino = node(abaco, italy, "trentino", "region")
belluno = node(abaco, veneto, "belluno", "city")
padova = node(abaco, veneto, "padova", "city")
rovereto = node(abaco, trentino, "rovereto", "city")
pergine = node(abaco, trentino, "pergine", "city")

#@info "TEST: veneto origins: $(abaco.origins[veneto.ne])"

#current_footprint = get_collected(abaco, italy.ne, "region.city.footprint")
#@debug "initial carbon footprints: $current_footprint"

ingest(abaco, ts, belluno.ne, Dict("footprint" => val1))

ingest(abaco, ts, padova.ne, Dict("footprint" => val2))
ingest(abaco, ts, rovereto.ne, Dict("footprint" => val3))

#current_footprint = get_collected(abaco, italy.ne, "region.city.footprint", ts)
#@debug "current carbon footprint: $current_footprint"

@test actual_triggers === expected_triggers


text:
# common types or functions

# the result type

immutable Result{T}
    U::Matrix{T}
    H::Matrix{T}
    niters::Integer
    converged::Bool

    function Result(U::Matrix{T}, H::Matrix{T},
                    niter::Integer,
                    converged::Bool)
        size(U, 2) == size(H, 1) ||
               throw(DimensionMismatch("Inner dimension of U and H mismatch."))
        new(U, H, niter, converged)
    end
end

immutable SVDResult{T}
    U::Matrix{T}
    H::Matrix{T}
    function SVDResult(U::Matrix{T}, H::Matrix{T})
        size(U, 2) == size(H, 1) ||
               throw(DimensionMismatch("Inner dimension of U and H mismatch."))
        new(U, H)
    end
end

# niters and converged are nothing for the SVD method

# the objective type

immutable Objective{T}
    absolute::T  # absolute error
    relative::T  # relative error

    function Objective(absolute::Real, relative::Real)
        absolute >= 0 || error("absolute must be non-negative.")
        relative >= 0 || error("relative must be non-negative.")

        new(absolute, relative)
    end
end

function evaluate_objv{T}(preU::Matrix{T}, U::Matrix{T})
    rel = norm(preU - U, Inf) / norm(preU, Inf)
    abs = norm(I - U'*U, Inf)
    return Objective{T}(rel, abs)
end

abstract PolarUpdater

abstract PolarAlg

# common algorithm skeleton for iterative updating methods

function common_iter!{T}(updater::PolarUpdater,
                         X::Matrix{T},
                         U::Matrix{T},
                         H::Matrix{T},
                         maxiter::Int,
                         verbose::Bool,
                         tol::T)

    preU = Array(T, size(X))
    copy!(U, X)
    converged = false
    t = 0
    if verbose
        @printf("%-5s    %-13s    %-13s\n", "Iter.", "Rel. err.",  "Obj.")
    end

    while !converged && t < maxiter
        t += 1
        copy!(preU, U)
        update_U!(updater, U)

        # determine convergence
        diff = vecnorm(preU - U)
        if diff < tol
            converged = true
        end

        # display infomation
        if verbose
            objv = evaluate_objv(preU, U)
            @printf("%5d    %13.6e    %13.6e\n",
                    t, objv.absolute, objv.relative)
        end
    end

    # compute H
    A_mul_B!(H, U', X)
    H = 0.5 * (H + H')
    return Result{T}(U, H, t, converged)
end

# Scaling iterative algorithm
function common_iter_scal!{T}(updater::PolarUpdater,
                              X::Matrix{T},
                              U::Matrix{T},
                              H::Matrix{T},
                              maxiter::Int,
                              verbose::Bool,
                              tol::T)

    preU = Array(T, size(X))
    copy!(U, X)
    converged = false
    t = 0
    if verbose
        @printf("%-5s    %-13s    %-13s\n", "Iter.", "Rel. err.",  "Obj.")
    end

    while !converged && t < maxiter
        t += 1
        copy!(preU, U)
        update_U!(updater, U)

        # determine convergence
        diff = vecnorm(preU - U)
        if diff < tol
            converged = true
        end

        # determine scaling
        reldiff = diff/vecnorm(U) # relative error
        if updater.scale && (reldiff < updater.scale_tol)
            updater.scale = false
        end

        # display infomation
        if verbose
            objv = evaluate_objv(preU, U)
            @printf("%5d    %13.6e    %13.6e\n",
                    t, objv.absolute, objv.relative)
        end
    end

    # compute H
    A_mul_B!(H, U', X)
    H = 0.5 * (H + H')
    return Result{T}(U, H, t, converged)
end

# Hybrid iteration algorithm
function common_iter_hybr!{T}(updater1::PolarUpdater,
                              updater2::PolarUpdater,
                              X::Matrix{T},
                              U::Matrix{T},
                              H::Matrix{T},
                              maxiter::Int,
                              verbose::Bool,
                              tol::T,
                              theta::T) # theta is the switch parameter

    preU = Array(T, size(X))
    copy!(U, X)
    converged = false
    switched = false
    t = 0
    if verbose
        @printf("%-5s    %-13s    %-13s\n", "Iter.", "Rel. err.",  "Obj.")
    end

    while !converged && t < maxiter
        t += 1
        copy!(preU, U)

        if switched
            update_U!(updater2, U)
        else
            obj = norm(I - U'*U, 1)
            if obj > theta # theta is the switch parameter
                update_U!(updater1, U)
            else
                switched = true
                update_U!(updater2, U)
            end
        end

        # determine convergence
        diff = vecnorm(preU - U)
        if diff < tol
            converged = true
        end

        # display infomation
        if verbose
            objv = evaluate_objv(preU, U)
            @printf("%5d    %13.6e    %13.6e\n",
                    t, objv.absolute, objv.relative)
        end
    end

    # compute H
    A_mul_B!(H, U', X)
    H = 0.5 * (H + H')
    return Result{T}(U, H, t, converged)
end


text:
function getMutationTypes(data::Array{ASCIIString,1})
#    sort(unique(vcat(map(x->split(x,';'), data))))
    data = sort(unique(mapreduce(x->split(x,';'),vcat, data)))
    data = map(x->convert(ASCIIString,x), data)
end

function getGenes(data::Array{ASCIIString,1})
    sort(unique(data))
end

function getSampleNames(data::Array{ASCIIString,1})
    sort(unique(data))
end

function outerProduct(vector1::Array{ASCIIString,1},vector2::Array{ASCIIString,1})
    num_vector1 = length(vector1)
    num_vector2 = length(vector2)

    data = Array{ASCIIString,1}(num_vector1*num_vector2)
    for i = 1:num_vector1
        for j = 1:num_vector2
            data[(i-1)*num_vector2+j] = string(vector1[i],":",vector2[j])
        end
    end

    sort(data)
end

function preprocess(input_file)
    data = readdlm(input_file, '\t', ASCIIString)

    mutation_types = getMutationTypes(data[:,3])
    num_mutation_types = length(mutation_types)
    genes = getGenes(data[:,2])
    num_genes = length(genes)

    gene_mutations = outerProduct(genes,mutation_types)
    num_gene_mutations = length(gene_mutations)
    row_names = gene_mutations
    feature_inds = Dict{ASCIIString,Int64}(zip(gene_mutations,1:num_gene_mutations))

    sample_names = getSampleNames(data[:,1])
    num_samples = length(sample_names)
    col_names = sample_names
    samplename_inds = Dict{ASCIIString,Int64}(zip(sample_names, 1:num_samples))

    #TODO here need to refactor
    samplename_feature_inds = Dict{Tuple{ASCIIString,ASCIIString},Int64}()
    for sample_name in sample_names
        for gene_mutation in gene_mutations
            samplename_feature_inds[(sample_name,gene_mutation)] =
                (feature_inds[gene_mutation]-1)*(num_samples) + samplename_inds[sample_name]
        end
    end

    info("There are $num_samples samples")
    info("There are $num_gene_mutations features")

    nrow,ncol = size(data)
    matrix = zeros(Float64, num_samples, num_gene_mutations)
    for i = 1:nrow
        sample_name, gene_name, mutation_type = data[i,:]
        key = (sample_name,string(gene_name,":",mutation_type))

        if in(key, keys(samplename_feature_inds))
            matrix[samplename_feature_inds[key]] = 1.0
        end
    end

    subfeat_inds = sum(matrix,1) .!= 0
    matrix = matrix[:,subfeat_inds]

    writecsv("../data/matrix.csv", matrix)
    writecsv("../data/sample_names.csv", sample_names)
    @show size(sample_names)
    @show size(subfeat_inds)
    @show size(gene_mutations)
    gene_mutations = gene_mutations[subfeat_inds[:]]

    gene_mutations = map(gene_mutations) do gemut
        len = length(gemut)
        if len <= 20
            return gemut
        else
            return gemut[len-20:end]
        end
    end

    writecsv("../data/gene_mutations.csv", gene_mutations)
end

preprocess("../data/guo_gene_heatmap.txt")


text:
using SeisXcorrelation
using Test

@testset "SeisXcorrelation.jl" begin
    # Write your own tests here.
    @test 1 == 1
end

max_stars_repo_name: kura-okubo/SeisXcorrelation

max_issues_repo_name: kura-okubo/SeisXcorrelation

max_forks_repo_name: kura-okubo/SeisXcorrelation


text:
using Test

using Catlab.Present
using Catlab.Theories: CatDescType, AttrDescType

# Schema
########

A, B = Ob(FreeSchema, :A), Ob(FreeSchema, :B)
f, g = Hom(:f, A, B), Hom(:g, B, A)

C = Data(FreeSchema.Data, :C)
x, y = Attr(:x, A, C), Attr(:y, B, C)

@test dom(x) == A
@test codom(x) == C
@test dom(compose(f,y)) == A
@test dom(compose(g,x)) == B

# CatDesc and AttrDesc
######################

@present TheoryDecGraph(FreeSchema) begin
  V::Ob
  E::Ob

  src::Hom(E,V)
  tgt::Hom(E,V)

  X::Data

  vdec::Attr(V,X)
  edec::Attr(E,X)
end

CD = CatDescType(TheoryDecGraph)
AD = AttrDescType(TheoryDecGraph)

@test dom(CD,:src) == :E
@test codom(CD,:tgt) == :V

@test dom(AD,:vdec) == :V
@test codom(AD,:edec) == :X


text:
"""
    transferoperator(triang::AbstractTriangulation;
                exact = false, parallel = true,
                n_pts = 200, sample_randomly = false)

Estimate the transfer operator from a triangulation.
"""
function transferoperator_triang(triang::AbstractTriangulation;
                            exact = false, parallel = true,
                            n_pts = 200, sample_randomly = false)
    if exact
        if parallel
            transferoperator_exact_p(triang)
        else
            transferoperator_exact(triang)
        end
    else
        transferoperator_approx(triang, n_pts = n_pts,
                                sample_randomly = sample_randomly)
    end
end


text:
module stats

export SimpleStats, set_state, update_state, get_stats, DeltaStats, test_state, edge_stats, directed_triplet_stats
using LinearAlgebra

mutable struct SimpleStats
    stats_function
    state

    function SimpleStats(stats_function)
        new(stats_function, undef)
    end
end

function set_state(stats :: SimpleStats, state)
    stats.state = copy(state)
end

function update_state(stats :: SimpleStats, update)
    i, x = update
    stats.state[i] = x
end

function test_state(stats :: SimpleStats, update)
    i, x = update
    old_x = stats.state[i]
    stats.state[i] = x
    test_stats = stats.stats_function(stats.state)
    stats.state[i] = old_x
    test_stats
end

function get_stats(stats :: SimpleStats)
    stats.stats_function(stats.state)
end

function get_stats(stats :: SimpleStats, state)
    stats.stats_function(state)
end

function Base.copy(stats :: SimpleStats)
    SimpleStats(stats.stats_function)
end

mutable struct DeltaStats
    stats_function
    delta_stats_function
    current_stats
    state

    function DeltaStats(stats_function, delta_stats_function)
        new(stats_function, delta_stats_function, undef, undef)
    end
end

function set_state(stats :: DeltaStats, state)
    stats.state = copy(state)
    stats.current_stats = stats.stats_function(state)
end

function update_state(stats :: DeltaStats, update)
    i, x = update
    stats.current_stats = stats.delta_stats_function(
        stats.state,
        stats.current_stats,
        update
    )
    stats.state[i] = x
end

function test_state(stats :: DeltaStats, update)
    i, x = update
    stats.delta_stats_function(
        stats.state,
        stats.current_stats,
        update
    )
end

function get_stats(stats :: DeltaStats)
    stats.current_stats
end

function get_stats(stats :: DeltaStats, state)
    stats.stats_function(state)
end

function Base.copy(stats :: DeltaStats)
    DeltaStats(stats.stats_function, stats.delta_stats_function)
end

s_string = [
"122223333444456",  # 0
"010000011111223",  # 1
"001001311232246",  # 2
"000101010211123",  # 3
"000011001112123",  # 4
"000001000212036",  # 5
"000000100010012",  # 6
"000000010210236",  # 7
"000000001012236",  # 8
"000000000100013",  # 9
"000000000010026",  # 10
"000000000001013",  # 11
"000000000000113",  # 12
"000000000000016",  # 13
"000000000000001"]
teto = vcat([[parse(Integer, c) for c ∈ r]' for r ∈ s_string]...)
tote = round.(Integer, teto^-1)

directed_triplet_stats = DeltaStats(
    function(G)
        s = zeros(Integer, 15)
        A = G.adjacency
        n, _ = size(A)
        m = sum(A)
        g_sq = A^2
        g_cu = g_sq * A
        tr_g_sq = tr(g_sq)
        s[1] = (n - 2) * m
        s[2] = (n - 2) * tr_g_sq ÷ 2
        s[3] = sum(g_sq) - tr_g_sq
        div = A' * A
        conv = A * A'
        s[4] = (sum(div) - m) ÷ 2
        s[5] = (sum(conv) - m) ÷ 2
        s[6] = sum(A .* g_sq)
        s[7] = tr(g_cu) ÷ 3
        g_sym = A .* A'
        bidegi = sum(g_sym, dims=2)
        s[8] = sum(A .* bidegi) - sum(g_sym)
        s[9] = sum(A .* bidegi') - sum(g_sym)
        s[10] = sum(g_sym .* conv) ÷ 2
        s[11] = sum(g_sym .* g_sq)
        s[12] = sum(g_sym .* div) ÷ 2
        g_sym_sq = g_sym^2
        s[13] = (sum(g_sym_sq) - tr(g_sym_sq)) ÷ 2
        s[14] = sum(g_sym_sq .* A)
        s[15] = tr(g_sym_sq * g_sym) ÷ 6
        c = tote * s
        c ./ (n * (n - 1) * (n - 2))
    end,
    function(G, s, up)
        A = G.adjacency
        n, _ = size(A)
        i, x = up
        u, v = i

        if A[u, v] == x
            return s
        end

        u, v = i
        Duv = 2 * A[u, v] - 1
        δ = zeros(Integer, 15)
        δ[1] = n - 2
        δ[2] = (n - 2) * A[v, u]
        outu = sum(A[u, :])
        inu = sum(A[:, u])
        outv = sum(A[v, :])
        inv = sum(A[:, v])
        δ[3] = inu + outv - 2 * A[v, u]
        δ[4] = outu - A[u, v]
        δ[5] = inv - A[u, v]

        common_post = A[u, :] .* A[v, :]
        common_pre = A[:, u] .* A[:, v]
        δ[6] = sum(common_post) + sum(common_pre) + dot(A[u, :], A[:, v])
        δ[7] = dot(A[v, :], A[:, u])

        bidegu = dot(A[u, :], A[:, u])
        bidegv = dot(A[v, :], A[:, v])
        δ[8] = bidegu + A[v, u] * (outu + outv - 2 * A[u, v] - 2 * A[v, u] + 1)
        δ[9] = bidegv + A[v, u] * (inu + inv - 2 * A[u, v] - 2 * A[v, u] + 1)

        δ[10] = dot(A[u, :], common_pre) + A[v, u] * sum(common_post)
        δ[11] = dot(common_post, A[:, u]) + dot(A[v, :], common_pre)
        δ[12] = dot(common_post, A[:, v]) + A[v, u] * sum(common_pre)

        δ[13] = A[v, u] * (bidegu + bidegv - 2 * A[u, v])
        δ[14] = dot(common_post, common_pre)

        if A[v, u] != 0
            δ[11] += dot(A[u, :], A[:, v]) + dot(A[v, :], A[:, u])
            δ[14] += dot(common_post, A[:, v]) + dot(common_post, A[:, u]) +
                     dot(A[u, :], common_pre) + dot(A[v, :], common_pre)
            δ[15] = dot(common_post, common_pre)
        end

        δ = tote * δ
        s - Duv * δ ./ (n * (n - 1) * (n - 2))
    end
)

end


text:
import Base: -, <, copysign, flipsign, convert

struct Double{T<:IEEEFloat} <: Number
    hi::T
    lo::T
end
Double(x::T) where {T<:IEEEFloat} = Double(x, zero(T))

(::Type{T})(x::Double{T}) where {T<:IEEEFloat} = x.hi + x.lo

@inline trunclo(x::Float64) = reinterpret(Float64, reinterpret(UInt64, x) & 0xffff_ffff_f800_0000) # clear lower 27 bits (leave upper 26 bits)
@inline trunclo(x::Float32) = reinterpret(Float32, reinterpret(UInt32, x) & 0xffff_f000) # clear lowest 12 bits (leave upper 12 bits)

@inline function splitprec(x::IEEEFloat)
    hx = trunclo(x)
    hx, x - hx
end

@inline function dnormalize(x::Double{T}) where {T}
    r = x.hi + x.lo
    Double(r, (x.hi - r) + x.lo)
end

@inline flipsign(x::Double{T}, y::T) where {T<:IEEEFloat} = Double(flipsign(x.hi, y), flipsign(x.lo, y))

@inline scale(x::Double{T}, s::T) where {T<:IEEEFloat} = Double(s * x.hi, s * x.lo)

@inline (-)(x::Double{T}) where {T<:IEEEFloat} = Double(-x.hi, -x.lo)

@inline function (<)(x::Double{T}, y::Double{T}) where {T<:IEEEFloat}
    x.hi < y.hi
end

@inline function (<)(x::Double{T}, y::Number) where {T<:IEEEFloat}
    x.hi < y
end

@inline function (<)(x::Number, y::Double{T}) where {T<:IEEEFloat}
    x < y.hi
end

# quick-two-sum x+y
@inline function dadd(x::T, y::T) where {T<:IEEEFloat} #WARNING |x| >= |y|
    s = x + y
    Double(s, (x - s) + y)
end

@inline function dadd(x::T, y::Double{T}) where {T<:IEEEFloat} #WARNING |x| >= |y|
    s = x + y.hi
    Double(s, (x - s) + y.hi + y.lo)
end

@inline function dadd(x::Double{T}, y::T) where {T<:IEEEFloat} #WARNING |x| >= |y|
    s = x.hi + y
    Double(s, (x.hi - s) + y + x.lo)
end

@inline function dadd(x::Double{T}, y::Double{T}) where {T<:IEEEFloat} #WARNING |x| >= |y|
    s = x.hi + y.hi
    Double(s, (x.hi - s) + y.hi + y.lo + x.lo)
end

@inline function dsub(x::Double{T}, y::Double{T}) where {T<:IEEEFloat} #WARNING |x| >= |y|
    s = x.hi - y.hi
    Double(s, (x.hi - s) - y.hi - y.lo + x.lo)
end

@inline function dsub(x::Double{T}, y::T) where {T<:IEEEFloat} #WARNING |x| >= |y|
    s = x.hi - y
    Double(s, (x.hi - s) - y + x.lo)
end

@inline function dsub(x::T, y::Double{T}) where {T<:IEEEFloat} #WARNING |x| >= |y|
    s = x - y.hi
    Double(s, (x - s) - y.hi - y.lo)
end

@inline function dsub(x::T, y::T) where {T<:IEEEFloat} #WARNING |x| >= |y|
    s = x - y
    Double(s, (x - s) - y)
end

# two-sum x+y  NO BRANCH
@inline function dadd2(x::T, y::T) where {T<:IEEEFloat}
    s = x + y
    v = s - x
    Double(s, (x - (s - v)) + (y - v))
end

@inline function dadd2(x::T, y::Double{T}) where {T<:IEEEFloat}
    s = x + y.hi
    v = s - x
    Double(s, (x - (s - v)) + (y.hi - v) + y.lo)
end

@inline dadd2(x::Double{T}, y::T) where {T<:IEEEFloat} = dadd2(y, x)

@inline function dadd2(x::Double{T}, y::Double{T}) where {T<:IEEEFloat}
    s = x.hi + y.hi
    v = s - x.hi
    Double(s, (x.hi - (s - v)) + (y.hi - v) + x.lo + y.lo)
end

@inline function dsub2(x::T, y::T) where {T<:IEEEFloat}
    s = x - y
    v = s - x
    Double(s, (x - (s - v)) + (-y - v))
end

@inline function dsub2(x::T, y::Double{T}) where {T<:IEEEFloat}
    s = x - y.hi
    v = s - x
    Double(s, (x - (s - v)) + (-y.hi - v) - y.lo)
end

@inline function dsub2(x::Double{T}, y::T) where {T<:IEEEFloat}
    s = x.hi - y
    v = s - x.hi
    Double(s, (x.hi - (s - v)) + (-y - v) + x.lo)
end

@inline function dsub2(x::Double{T}, y::Double{T}) where {T<:IEEEFloat}
    s = x.hi - y.hi
    v = s - x.hi
    Double(s, (x.hi - (s - v)) + (-y.hi - v) + x.lo - y.lo)
end

if FMA_FAST

    # two-prod-fma
    @inline function dmul(x::T, y::T) where {T<:IEEEFloat}
        z = x * y
        Double(z, fma(x, y, -z))
    end

    @inline function dmul(x::Double{T}, y::T) where {T<:IEEEFloat}
        z = x.hi * y
        Double(z, fma(x.hi, y, -z) + x.lo * y)
    end

    @inline dmul(x::T, y::Double{T}) where {T<:IEEEFloat} = dmul(y, x)

    @inline function dmul(x::Double{T}, y::Double{T}) where {T<:IEEEFloat}
        z = x.hi * y.hi
        Double(z, fma(x.hi, y.hi, -z) + x.hi * y.lo + x.lo * y.hi)
    end

    # x^2
    @inline function dsqu(x::T) where {T<:IEEEFloat}
        z = x * x
        Double(z, fma(x, x, -z))
    end

    @inline function dsqu(x::Double{T}) where {T<:IEEEFloat}
        z = x.hi * x.hi
        Double(z, fma(x.hi, x.hi, -z) + x.hi * (x.lo + x.lo))
    end

    # sqrt(x)
    @inline function dsqrt(x::Double{T}) where {T<:IEEEFloat}
        zhi = _sqrt(x.hi)
        Double(zhi, (x.lo + fma(-zhi, zhi, x.hi)) / (zhi + zhi))
    end

    # x/y
    @inline function ddiv(x::Double{T}, y::Double{T}) where {T<:IEEEFloat}
        invy = 1 / y.hi
        zhi = x.hi * invy
        Double(zhi, (fma(-zhi, y.hi, x.hi) + fma(-zhi, y.lo, x.lo)) * invy)
    end

    @inline function ddiv(x::T, y::T) where {T<:IEEEFloat}
        ry = 1 / y
        r = x * ry
        Double(r, fma(-r, y, x) * ry)
    end

    # 1/x
    @inline function drec(x::T) where {T<:IEEEFloat}
        zhi = 1 / x
        Double(zhi, fma(-zhi, x, one(T)) * zhi)
    end

    @inline function drec(x::Double{T}) where {T<:IEEEFloat}
        zhi = 1 / x.hi
        Double(zhi, (fma(-zhi, x.hi, one(T)) + -zhi * x.lo) * zhi)
    end

else

    #two-prod x*y
    @inline function dmul(x::T, y::T) where {T<:IEEEFloat}
        hx, lx = splitprec(x)
        hy, ly = splitprec(y)
        z = x * y
        Double(z, ((hx * hy - z) + lx * hy + hx * ly) + lx * ly)
    end

    @inline function dmul(x::Double{T}, y::T) where {T<:IEEEFloat}
        hx, lx = splitprec(x.hi)
        hy, ly = splitprec(y)
        z = x.hi * y
        Double(z, (hx * hy - z) + lx * hy + hx * ly + lx * ly + x.lo * y)
    end

    @inline dmul(x::T, y::Double{T}) where {T<:IEEEFloat} = dmul(y, x)

    @inline function dmul(x::Double{T}, y::Double{T}) where {T<:IEEEFloat}
        hx, lx = splitprec(x.hi)
        hy, ly = splitprec(y.hi)
        z = x.hi * y.hi
        Double(z, (((hx * hy - z) + lx * hy + hx * ly) + lx * ly) + x.hi * y.lo + x.lo * y.hi)
    end

    # x^2
    @inline function dsqu(x::T) where {T<:IEEEFloat}
        hx, lx = splitprec(x)
        z = x * x
        Double(z, (hx * hx - z) + lx * (hx + hx) + lx * lx)
    end

    @inline function dsqu(x::Double{T}) where {T<:IEEEFloat}
        hx, lx = splitprec(x.hi)
        z = x.hi * x.hi
        Double(z, (hx * hx - z) + lx * (hx + hx) + lx * lx + x.hi * (x.lo + x.lo))
    end

    # sqrt(x)
    @inline function dsqrt(x::Double{T}) where {T<:IEEEFloat}
        c = _sqrt(x.hi)
        u = dsqu(c)
        Double(c, (x.hi - u.hi - u.lo + x.lo) / (c + c))
    end

    # x/y
    @inline function ddiv(x::Double{T}, y::Double{T}) where {T<:IEEEFloat}
        invy = 1 / y.hi
        c = x.hi * invy
        u = dmul(c, y.hi)
        Double(c, ((((x.hi - u.hi) - u.lo) + x.lo) - c * y.lo) * invy)
    end

    @inline function ddiv(x::T, y::T) where {T<:IEEEFloat}
        ry = 1 / y
        r = x * ry
        hx, lx = splitprec(r)
        hy, ly = splitprec(y)
        Double(r, (((-hx * hy + r * y) - lx * hy - hx * ly) - lx * ly) * ry)
    end

    # 1/x
    @inline function drec(x::T) where {T<:IEEEFloat}
        c = 1 / x
        u = dmul(c, x)
        Double(c, (one(T) - u.hi - u.lo) * c)
    end

    @inline function drec(x::Double{T}) where {T<:IEEEFloat}
        c = 1 / x.hi
        u = dmul(c, x.hi)
        Double(c, (one(T) - u.hi - u.lo - c * x.lo) * c)
    end

end


text:
module CelesteModSampleEntity

using ..Ahorn, Maple

@mapdef Entity "CelesteMod/SampleEntity" SampleEntity(
    x::Integer, y::Integer,
    sampleProperty::Integer=0
)

const placements = Ahorn.PlacementDict(
    "Sample Entity (CelesteMod)" => Ahorn.EntityPlacement(
        SampleEntity,
    )
)

sprite = "objects/CelesteMod/sampleEntity/idle00"

function Ahorn.selection(entity::SampleEntity)
    x, y = Ahorn.position(entity)
    return Ahorn.getSpriteRectangle(sprite, x, y)
end

Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::SampleEntity, room::Maple.Room) = Ahorn.drawSprite(ct, sprite, 0, 0)

end

max_stars_repo_name: swoolcock/CelesteModTemplate

max_issues_repo_name: swoolcock/CelesteModTemplate

max_forks_repo_name: swoolcock/CelesteModTemplate


text:
"""
Given the inverse matrix iM = [A B; C D] of a matrix M = [U V; W X]
(in block form), this returns inv(X). This is efficient if dim(A) ≪ dim(X),
because it requires only the inversion of A. The indices `i`  (an iterator
of integers) determine the sub-matrix X = M[i,i]. See:
Bernstein, Dennis S. 2009. Matrix Mathematics: Theory, Facts,
and Formulas. 2nd ed. Princeton Univ. Press. Eqs. 2.8.16-17.

Arguments:
    iM  matrix inverse of M
    i   index vector of the block X = M[i,i] that will be inverted
    iA  (optional) inverse of A if you have it.
"""
function blockinv end

function blockinv(iM::AbstractMatrix, i)
    o = setdiff(1 : size(iM, 1), i)
    blockinv(iM, i, inv(iM[o,o]))
end

function blockinv(iM::Hermitian, i)::Hermitian
    o = setdiff(1 : size(iM, 1), i)
    blockinv(iM, i, inv(Hermitian(iM[o,o])))
end

function blockinv(iM::Hermitian, i, iA::Hermitian)::Hermitian
    @assert allunique(i) && all(1 .≤ i .≤ size(iM, 1))
    o = setdiff(1 : size(iM, 1), i)
    @assert size(iA) == (length(o), length(o))

    B = iM[o,i]
    D = Hermitian(iM[i,i])

    Hermitian(D - Hermitian(B' * iA * B))
end

function blockinv(iM::AbstractMatrix, i, iA::AbstractMatrix)
    @assert allunique(i) && all(1 .≤ i .≤ size(iM, 1))
    o = setdiff(1 : size(iM, 1), i)
    @assert size(iA) == (length(o), length(o))

    B = iM[o,i]
    C = iM[i,o]
    D = iM[i,i]

    D - C * iA * B
end


text:
using EngEconomics

# Given
A = 10000
gPrime = 0.01
g = sqrt(1 + gPrime) - 1
n = 2
i = .18
m = 12

ic = i / m
nc = 2 * m

# Find: Amount of profit after 2 years
P = A * presentWorthConversionFactorGeometric(g, ic, nc)


text:
#################
### Key class ###
#################

immutable Key
  root    :: ModInt{12}
  ismajor :: Bool
end

const key_dict = TwoWayDict(
    Key(0,  true) => "C",
    Key(1,  true) => "Db",
    Key(2,  true) => "D",
    Key(3,  true) => "Eb",
    Key(4,  true) => "E",
    Key(5,  true) => "F",
    Key(6,  true) => "Gb",
    Key(7,  true) => "G",
    Key(8,  true) => "Ab",
    Key(9,  true) => "A",
    Key(10, true) => "Bb",
    Key(11, true) => "B",
    Key(0, false) => "c",
    Key(1, false) => "c#",
    Key(2, false) => "d",
    Key(3, false) => "eb",
    Key(4, false) => "e",
    Key(5, false) => "f",
    Key(6, false) => "f#",
    Key(7, false) => "g",
    Key(8, false) => "g#",
    Key(9, false) => "a",
    Key(10,false) => "bb",
    Key(11,false) => "b")

Key(root::ModInt{12}, mode::AbstractString) =
    Key(root, mode == "major" ? true : false)
Key(root::Int, mode::AbstractString) =
    Key(ModInt{12}(root), mode == "major" ? true : false)
Key(s::AbstractString) = key_dict[s]

show(io::IO, k::Key) = print(io, key_dict[k])
root(key::Key) = key.root
ismajor(key::Key) = key.ismajor
isminor(k::Key) = !ismajor(k)
mode(k::Key)    = ismajor(k) ? "major" : "minor"
change_mode(k::Key) = Key(root(k), !ismajor(k))

function uses_flats(k::Key)
    if ismajor(k)
        root(k) * 7 > ModInt{12}(5)
    else
        (root(k) + 3) * 7 > ModInt{12}(5)
    end
end

+(k::Key, i) = Key(root(k) + i, ismajor(k))
-(k::Key, i) = Key(root(k) - i, ismajor(k))

###########################
### Scalae degree class ###
###########################

immutable ScaleDeg
    deg :: ModInt{7}
end

const roman_dict = TwoWayDict(
    "I"  => ModInt{7}(0),
    "II" => ModInt{7}(1),
    "III"=> ModInt{7}(2),
    "IV" => ModInt{7}(3),
    "V"  => ModInt{7}(4),
    "VI" => ModInt{7}(5),
    "VII"=> ModInt{7}(6))

ScaleDeg(str::String) = ScaleDeg(roman_dict[str])

numeral(s::ScaleDeg) = roman_dict[s.deg]
show(io::IO, s::ScaleDeg) = print(io, numeral(s))

+(s::ScaleDeg, i::Int) = ScaleDeg(s.deg + i)
-(s::ScaleDeg, i::Int) = ScaleDeg(s.deg - i)

const major_pattern = (0,2,4,5,7,9,11)
const minor_pattern = (0,2,3,5,7,8,10)

function root(s::ScaleDeg, k::Key) :: ModInt{12}
    pattern = ismajor(k) ? major_pattern : minor_pattern
    pattern[s.deg] + root(k)
end

const natural_tone_dict = TwoWayDict(
    "C" => ModInt{12}(0),
    "D" => ModInt{12}(2),
    "E" => ModInt{12}(4),
    "F" => ModInt{12}(5),
    "G" => ModInt{12}(7),
    "A" => ModInt{12}(9),
    "B" => ModInt{12}(11))

function spelled_root(s::ScaleDeg, k::Key) :: String
    pc = root(s, k) # pitch class of the root
    if pc in major_pattern # then it is the pitch class of a natural tone
        natural_tone_dict[pc]
    elseif uses_flats(k)
        natural_tone_dict[pc+1] * "b"
    else # key uses sharps
        natural_tone_dict[pc-1] * "#"
    end
end

const major_modes = ("major","minor","minor","major","major","minor","minor")
const minor_modes = ("minor","minor","major","minor","minor","major","major")
mode(s::ScaleDeg, k::Key) = (ismajor(k) ? major_modes : minor_modes)[s.deg]

modulate(s::ScaleDeg, k::Key) = Key(root(s, k), mode(s, k))

const major_triad_forms = ("maj","min","min","maj","maj","min","dim")
const minor_triad_forms = ("min","hdim","maj","min","min","maj","maj")
chord_form(s::ScaleDeg, k::Key) = (ismajor(k) ? major_triad_forms : minor_triad_forms)[s.deg]

terminal(s::ScaleDeg, k::Key, form::String=chord_form(s)) = string(Int(root(s, k)), form)

######################
### Category class ###
######################

abstract type AbstractCategory end

const chord_form_classification = Dict(
    "13b9" => "maj",
    "7sus4" => "maj",
    "o7" => "dim",
    "7b9" => "maj",
    "sus4" => "sus4",
    "m11" => "min",
    "m" => "min",
    "7" => "maj",
    "9#11" => "maj",
    "%7" => "hdim",
    "^7#11" => "maj",
    "13" => "maj",
    "+7" => "aug",
    "^6" => "maj",
    "+^7" => "aug",
    "9" => "maj",
    "^9#11" => "maj",
    "+" => "aug",
    "" => "maj",
    "m9" => "min",
    "m7" => "min",
    "^13" => "maj",
    "mb6" => "min",
    "11" => "maj",
    "13#9" => "maj",
    "o" => "dim",
    "^7" => "maj",
    "13" => "maj",
    "m6" => "min",
    "m6" => "min",
    "7b13" => "maj",
    "+9" => "maj",
    "7alt" => "alt",
    "7#11" => "maj",
    "9" => "maj",
    "13#11" => "maj",
    "7b5" => "maj",
    "7b9" => "maj",
    "m^7" => "min",
    "^9" => "maj"
)

immutable Terminal <: AbstractCategory
    symbol :: String
    root   :: ModInt{12}
    form   :: String
end
==(t1::Terminal, t2::Terminal) = t1.root == t2.root && t1.form == t2.form
hash(t::Terminal) = hash("Terminal", hash(t.root, hash(t.form)))

Terminal(s::ScaleDeg, k::Key, form::String=chord_form(s, k)) =
    Terminal(spelled_root(s, k) * form, root(s, k), form)
function Terminal(str::AbstractString)
    if ismatch(r"^(\d+)(.*)", str)
        m = match(r"^(\d*)(.*)", str)
        Terminal(str, parse(m[1]), m[2])
    elseif ismatch(r"^([A-G][b#]*)(.*)", str)
        m = match(r"^([A-G][b#]*)(.*)", str)
        pc = pitch_class(m[1])
        form = chord_form_classification[m[2]]
        Terminal(str, pc, form)
    end
end
show(io::IO, t::Terminal) = print(io, t.symbol)

function make_terminals(scaledegs, key)
   lst = split(scaledegs)
   [Terminal(ScaleDeg(String(s)),Key(key)) for s in lst]
end

immutable Category <: AbstractCategory
  scaledeg   :: ScaleDeg
  key        :: Key
  isdomsubst :: Bool
end
Category(str::String, k::Key, b::Bool) = Category(ScaleDeg(str), k, b)
Category(deg::ScaleDeg, str::String, b::Bool) = Category(deg, Key(str), b)
Category(str1::String, str2::String, b::Bool) = Category(ScaleDeg(str1), Key(str2), b)

show(io::IO, cat::Category) =
    print(io, scaledeg(cat), "_{", key(cat), "}", isdomsubst(cat) ? "T" : "F")

scaledeg(cat::Category) = cat.scaledeg
key(cat::Category) = cat.key
isdomsubst(cat::Category) = cat.isdomsubst

keyroot(cat::Category) = root(key(cat))
keymode(cat::Category) = mode(key(cat))
root(cat::Category) = root(scaledeg(cat), key(cat))
modulate(cat::Category) = modulate(scaledeg(cat), key(cat))
change_mode(cat::Category) = Category(scaledeg(cat), change_mode(key(cat)), isdomsubst(cat))

Terminal(cat::Category) = Terminal(scaledeg(cat), key(cat))
Terminal(cat::Category, form::AbstractString) = Terminal(scaledeg(cat), key(cat), form)

notdomsubst(cat::Category) = Category(scaledeg(cat), key(cat), false)

# diatonic transposition
transpose(cat::Category, i) = Category(scaledeg(cat) + i, key(cat), isdomsubst(cat))


text:
using Test

using ProximalOperators

using ProximalOperators:
    ArrayOrTuple,
    is_prox_accurate,
    is_separable,
    is_convex,
    is_singleton,
    is_cone,
    is_affine,
    is_set,
    is_smooth,
    is_quadratic,
    is_generalized_quadratic,
    is_strongly_convex,
    is_positively_homogeneous,
    is_support

function call_test(f, x::ArrayOrTuple{R}) where R <: Real
    try
        fx = f(x)
        @test typeof(fx) == R
        return fx
    catch e
        if !isa(e, MethodError)
            return nothing
        end
    end
end

# tests equality of the results of prox, prox! and prox_naive
function prox_test(f, x::ArrayOrTuple{R}, gamma=R(1)) where R <: Real
    y, fy = prox(f, x, gamma)

    @test typeof(fy) == R

    y_prealloc = zero(x)
    fy_prealloc = prox!(y_prealloc, f, x, gamma)

    @test typeof(fy_prealloc) == R

    y_naive, fy_naive = ProximalOperators.prox_naive(f, x, gamma)

    @test typeof(fy_naive) == R

    rtol = if ProximalOperators.is_prox_accurate(f) sqrt(eps(R)) else 1e-4 end

    if ProximalOperators.is_convex(f)
        @test all(isapprox.(y_prealloc, y, rtol=rtol, atol=100*eps(R)))
        @test all(isapprox.(y_naive, y, rtol=rtol, atol=100*eps(R)))
        if ProximalOperators.is_set(f)
            @test fy_prealloc == 0
        end
        @test isapprox(fy_prealloc, fy, rtol=rtol, atol=100*eps(R))
        @test isapprox(fy_naive, fy, rtol=rtol, atol=100*eps(R))
    end

    if !ProximalOperators.is_set(f) || ProximalOperators.is_prox_accurate(f)
        f_at_y = call_test(f, y)
        if f_at_y !== nothing
            @test isapprox(f_at_y, fy, rtol=rtol, atol=100*eps(R))
        end
    end

    return y, fy
end

# tests equality of the results of prox, prox! and prox_naive
function gradient_test(f, x::ArrayOrTuple{R}, gamma=R(1)) where R <: Real
    grad_fx, fx = gradient(f, x)
    @test typeof(fx) == R
    return grad_fx, fx
end

# test predicates consistency
# i.e., that more specific properties imply less specific ones
# e.g., the indicator of a subspace is the indicator of a set in particular
function predicates_test(f)
    # quadratic => generalized_quadratic && smooth
    @test !is_quadratic(f) || (is_generalized_quadratic(f) && is_smooth(f))
    # (singleton || cone || affine) => set
    @test !(is_singleton(f) || is_cone(f) || is_affine(f)) || is_set(f)
    # cone => positively homogeneous
    @test !is_cone(f) || is_positively_homogeneous(f)
    # (convex && positively homogeneous) <=> (convex && support)
    @test (is_convex(f) && is_positively_homogeneous(f)) == (is_convex(f) && is_support(f))
    # strongly_convex => convex
    @test !is_strongly_convex(f) || is_convex(f)
end

@testset "Sanity checks" begin
    @test isempty(detect_unbound_args(ProximalOperators))
end

@testset "Utilities" begin
    include("test_symmetricpacked.jl")
end

@testset "Functions" begin
    include("test_cubeNormL2.jl")
    include("test_huberLoss.jl")
    include("test_indAffine.jl")
    include("test_leastSquares.jl")
    include("test_logisticLoss.jl")
    include("test_quadratic.jl")
    include("test_linear.jl")
    include("test_indHyperslab.jl")
    include("test_graph.jl")
    include("test_normL1plusL2.jl")
end

include("test_calls.jl")
include("test_indPolyhedral.jl")

@testset "Gradients" begin
    include("test_gradients.jl")
end

@testset "Calculus rules" begin
    include("test_calculus.jl")
    include("test_epicompose.jl")
    include("test_moreauEnvelope.jl")
    include("test_precompose.jl")
    include("test_pointwiseMinimum.jl")
    include("test_postcompose.jl")
    include("test_regularize.jl")
    include("test_separableSum.jl")
    include("test_slicedSeparableSum.jl")
    include("test_sum.jl")
end

@testset "Equivalences" begin
    include("test_equivalences.jl")
end

include("test_optimality_conditions.jl")

@testset "Hardcoded" begin
    include("test_results.jl")
end

@testset "Demos" begin
    include("../demos/lasso.jl")
    include("../demos/rpca.jl")
end


text:
function which_room(env::FourRooms, state)
    frp = FourRoomsContParams
    room = -1
    if state[1] < 6
        # LEFT
        if state[2] < 6
            # TOP
            room = frp.ROOM_TOP_LEFT
        else
            # Bottom
            room = frp.ROOM_BOTTOM_LEFT
        end
    else
        # RIGHT
        if state[2] < 7
            # TOP
            room = frp.ROOM_TOP_RIGHT
        else
            # Bottom
            room = frp.ROOM_BOTTOM_RIGHT
        end
    end
    return room
end


text:
# Place here configuration options that will be set for all environments
if isdefined(:config)
end

max_stars_repo_name: essenciary/genie-demo-hello-world

max_issues_repo_name: essenciary/genie-demo-hello-world

max_forks_repo_name: essenciary/genie-demo-hello-world


text:
@testset "Test 2d vem poisson" begin

    using jFEMTools
    import Tensors

    mesh = unitSquareMesh(RectangleCell, (1,1))

    # forcing function
    rhs(x::Tensors.Vec{2}) = 2*π^2*sin(π*x[1])*sin(π*x[2])

    # Test for order = 1
    degree = 1;
    dim = 2;
    element = PoissonVirtualElement(dim,degree);
    Vh = VEMFunctionSpace(mesh,element)
    u = TrialFunction(Vh)
    dofs = DofHandler(mesh, [u]);
    operators = VEMOperators(dofs, u;load = rhs);
    # Test
    d = sqrt(2)
    Bref = 1/4*[1 1 1 1;
                -d d -d d;
                -d -d d d]
    @test Bref ≈ operators.elements[1].B

    Dref = 1/4*[4 -d -d;
                4 d -d;
                4 d  d;
                4 -d d]

    @test Dref ≈ operators.elements[1].D

    Gref = 1/2*[2 0 0;
                0 1 0;
                0 0 1]

    @test Gref ≈ operators.elements[1].G

    # Test order 2 implementation
    degree = 2
    dim = 2
    element = PoissonVirtualElement(dim,degree)
    u = TrialFunction(VEMFunctionSpace(mesh,element))
    dofs = DofHandler(mesh, [u])
    operators = VEMOperators(dofs, u;load = rhs)
    # Test
    d = sqrt(2)
    Bref = 1/12*[0 0 0 0 0 0 0  0 12;
                 -d d d -d 0 4*d 0 -4*d 0;
                 -d -d d d -4*d 0 4*d 0 0;
                 1 1 1 1 0 4 0 4 -12;
                 1 -1 1 -1 0 0 0 0 0;
                 1 1 1 1 4 0 4 0 -12]
    @test Bref ≈ operators.elements[1].B

    Dref = 1/24*[24 -6*d -6*d 3 3 3;
                24 6*d -6*d 3 -3 3;
                24 6*d 6*d 3 3 3;
                24 -6*d 6*d 3 -3 3;
                24 0 -6*d 0 0 3;
                24 6*d 0 3 0 0;
                24 0 6*d 0 0 3;
                24 -6*d 0 3 0 0;
                24 0 0 1 0 1]

    @test Dref ≈ operators.elements[1].D

    Gref = 1/24*[24 0 0 1 0 1;
                0 12 0 0 0 0;
                0 0 12 0 0 0;
                0 0 0 2 0 0;
                0 0 0 0 1 0;
                0 0 0 0 0 2]

    @test Gref ≈ operators.elements[1].G
end

max_stars_repo_name: Paulms/jFEMToools

max_issues_repo_name: Paulms/jFEMToools

max_forks_repo_name: Paulms/jFEMToools


text:
# Disable updating registry on `add` (still runs on `up`)
using Pkg: Pkg
Pkg.UPDATED_REGISTRY_THIS_SESSION[] = true

# use more colors for displaying Pkg conflicts
if isdefined(Pkg.Resolve, :CONFLICT_COLORS)
    append!(Pkg.Resolve.CONFLICT_COLORS,  [21:51; 55:119; 124:142; 160:184; 196:220])
end

# Setup Revise.jl
atreplinit() do repl
    try
        @eval using Revise
        @async Revise.wait_steal_repl_backend()
    catch err
        println("Error starting Revise $err")
    end
end

max_stars_repo_name: oxinabox/dotfiles

max_issues_repo_name: oxinabox/dotfiles

max_forks_repo_name: oxinabox/dotfiles


text:
maximum(abs.(u[4,1,1,:].- mean(u[4,1,1,:],1))) ./ abs.(mean(u[4,1,1,:],1)
mean(u[4,1,1,:],1))
plot(abs(u[4,1,1,:]))#first harmonic

u = getMeasurements(f, false, frames=measFGFrameIdx(f),
                   fourierTransform=true, bgCorrection=true,
                    tfCorrection=false)


text:
"""
Return an array containing a rank ordering of the processable attributes in the data file.
For each attribute, its index, name, type, rank value, and optimized slices (for continuous attributes)
are included. The array can be printed with printattributeranking().

function generateattributeranking(dataFile)
	table, codes, classValues = generatetraintesttable(dataFile)
	classes = unique(classValues)
	@show table codes classValues
	# loop thru the attributes
	for i in eachindex(codes)
		d = dropnull(table[i])
		uniques=unique(d)
		@show i uniques length(uniques)
		# loop thru unique values
		for j in eachindex(uniques)
			# make a list of all values for each class
			# note the use of get() to enable the == operator; the default "" will always fail
			@show uniques[j]
			# loop thru each class
			for k in eachindex(classes)
				@show classes[k]
				# loop thru each case
				for m in eachindex(classValues)
					# test whether
					if get(classes[k],"") == get(classValues[m],"")
						@show [get(v,"") == uniques[j] for v in table[i]]
					end
				end
			end
		end
	end
end
"""
function generateattributeranking(dataFile)
	table, codes, classValues = generatetraintesttable(dataFile)
	classes = unique(classValues)
	labels = names(table)
	# @show labels table
	A = falses(length(codes),length(classValues),length(classes),length(classValues))
	B = zeros(Int32,length(codes),length(classValues),length(classes))
	for i in indices(A,1) # attributes
		uniques = unique(dropnull(table[i]))
		values = table[i]
		# @show i labels[i] values uniques
		for j in eachindex(uniques) # unique values for each attribute
			# @show j uniques[j]
			for k in indices(A,3) # classes
				# @show k classes[k]
				for m in indices(A,4) # cases
					# @show m classValues[m] values[m]
					# set an element for a particular class k and unique value j of an attribute i
					# to true if the class for its case m is the same as class k
					equalsvaluesuniques(x::T,y::T) where {T<:Any} = x == y
					equalsvaluesuniques(x::T,y::T) where {T<:Integer} = get(x,0) == get(y,0)
					equalsvaluesuniques(x::T,y::T) where {T<:Rational} = get(x,0.) == get(y,0.)
					equalsvaluesuniques(x::Nullable{String},y::String) = get(x,"") == y
					equalsvaluesuniques(x::Nullable{Int64},y::Int64) = get(x,0) == y
					equalsvaluesuniques(x::Nullable{Float64},y::Float64) = get(x,0.) == y
					# @show equalsvaluesuniques(values[m], uniques[j])
					A[i,j,k,m] = get(classes[k],"") == get(classValues[m],"") && equalsvaluesuniques(values[m], uniques[j])
				end
				# @show A[i,j,k,:] sum(A[i,j,k,:])
				# sum along the m dimension; put the result in B
				B[i,j,k] = sum(A[i,j,k,:])
			end
		end
	end
	# @show size(A) A
	# @show size(B) B
	# @show size(B,1)
	# for each attribute, get the absolute values of the differences between each class
	C = Array{Int64,2}(length(codes), 2)
	for i in indices(B,1) # attributes
		rank = 0
		for j in indices(B,2) # cases
			for k in indices(B,3) # classes
#rank1 = [abs(x - y) for (x,y) in zip(B)]
				for m in k:size(B,3)
					# @show i j k m
					rank += abs(B[i,j,k] - B[i,j,m])
				end
			end
		end
		# @show i labels[i] rank
		C[i,1] = rank
		C[i,2] = i
	end
	D = sortrows(C,rev=true)
	# @show labels table D
	return labels, codes, table, D
end


text:
using OrdinaryDiffEq
using Trixi

"""
    initial_condition_jeans_instability(x, t,
                                        equations::Union{CompressibleEulerEquations2D,
                                                         HyperbolicDiffusionEquations2D})

The classical Jeans instability taken from
- Michael Schlottke-Lakemper, Andrew R. Winters, Hendrik Ranocha, Gregor J. Gassner (2020)
  A purely hyperbolic discontinuous Galerkin approach for self-gravitating gas dynamics
  [arXiv: 2008.10593](https://arxiv.org/abs/2008.10593)
- Dominik Derigs, Andrew R. Winters, Gregor J. Gassner, Stefanie Walch (2016)
  A Novel High-Order, Entropy Stable, 3D AMR MHD Solver with Guaranteed Positive Pressure
  [arXiv: 1605.03572](https://arxiv.org/abs/1605.03572)
- Flash manual https://flash.uchicago.edu/site/flashcode/user_support/flash_ug_devel.pdf
in CGS (centimeter, gram, second) units.
"""
function initial_condition_jeans_instability(x, t,
                                            equations::CompressibleEulerEquations2D)
  # Jeans gravitational instability test case
  # see Derigs et al. https://arxiv.org/abs/1605.03572; Sec. 4.6
  # OBS! this uses cgs (centimeter, gram, second) units
  # periodic boundaries
  # domain size [0,L]^2 depends on the wave number chosen for the perturbation
  # OBS! Be very careful here L must be chosen such that problem is periodic
  # typical final time is T = 5
  # gamma = 5/3
  dens0  = 1.5e7 # g/cm^3
  pres0  = 1.5e7 # dyn/cm^2
  delta0 = 1e-3
  # set wave vector values for pertubation (units 1/cm)
  # see FLASH manual: https://flash.uchicago.edu/site/flashcode/user_support/flash_ug_devel.pdf
  kx = 2.0*pi/0.5 # 2π/λ_x, λ_x = 0.5
  ky = 0.0   # 2π/λ_y, λ_y = 1e10
  k_dot_x = kx*x[1] + ky*x[2]
  # perturb density and pressure away from reference states ρ_0 and p_0
  dens = dens0*(1.0 + delta0*cos(k_dot_x))                 # g/cm^3
  pres = pres0*(1.0 + equations.gamma*delta0*cos(k_dot_x)) # dyn/cm^2
  # flow starts as stationary
  velx = 0.0 # cm/s
  vely = 0.0 # cm/s
  return prim2cons((dens, velx, vely, pres), equations)
end

function initial_condition_jeans_instability(x, t,
                                             equations::HyperbolicDiffusionEquations2D)
  # gravity equation: -Δϕ = -4πGρ
  # Constants taken from the FLASH manual
  # https://flash.uchicago.edu/site/flashcode/user_support/flash_ug_devel.pdf
  rho0   = 1.5e7
  delta0 = 1e-3

  phi = rho0*delta0 # constant background pertubation magnitude
  q1  = 0.0
  q2  = 0.0
  return (phi, q1, q2)
end

initial_condition = initial_condition_jeans_instability

###############################################################################
# semidiscretization of the compressible Euler equations
gamma = 5/3
equations_euler = CompressibleEulerEquations2D(gamma)

polydeg = 3
solver_euler = DGSEM(polydeg, flux_hll)

coordinates_min = (0, 0)
coordinates_max = (1, 1)
mesh = TreeMesh(coordinates_min, coordinates_max,
                initial_refinement_level=4,
                n_cells_max=10_000)

semi_euler = SemidiscretizationHyperbolic(mesh, equations_euler, initial_condition, solver_euler)

###############################################################################
# semidiscretization of the hyperbolic diffusion equations
equations_gravity = HyperbolicDiffusionEquations2D()

solver_gravity = DGSEM(polydeg, flux_lax_friedrichs)

semi_gravity = SemidiscretizationHyperbolic(mesh, equations_gravity, initial_condition, solver_gravity,
                                            source_terms=source_terms_harmonic)

###############################################################################
# combining both semidiscretizations for Euler + self-gravity
parameters = ParametersEulerGravity(background_density=1.5e7, # aka rho0
                                    gravitational_constant=6.674e-8, # aka G
                                    cfl=1.6,
                                    resid_tol=1.0e-4,
                                    n_iterations_max=1000,
                                    timestep_gravity=timestep_gravity_carpenter_kennedy_erk54_2N!)

semi = SemidiscretizationEulerGravity(semi_euler, semi_gravity, parameters)

###############################################################################
# ODE solvers, callbacks etc.
tspan = (0.0, 5.0)
ode = semidiscretize(semi, tspan);

summary_callback = SummaryCallback()

stepsize_callback = StepsizeCallback(cfl=1.0)

save_solution = SaveSolutionCallback(interval=10,
                                     save_initial_solution=true,
                                     save_final_solution=true,
                                     solution_variables=cons2prim)

save_restart = SaveRestartCallback(interval=100,
                                   save_final_restart=true)

analysis_interval = 100
alive_callback = AliveCallback(analysis_interval=analysis_interval)

Trixi.pretty_form_utf(::Val{:energy_potential}) = "∑e_potential"
Trixi.pretty_form_ascii(::Val{:energy_potential}) = "e_potential"

function Trixi.analyze(::Val{:energy_potential}, du, u_euler, t, semi::SemidiscretizationEulerGravity)

  u_gravity = Trixi.wrap_array(semi.cache.u_ode, semi.semi_gravity)

  mesh, equations_euler, dg, cache = Trixi.mesh_equations_solver_cache(semi.semi_euler)
  _, equations_gravity, _, _ = Trixi.mesh_equations_solver_cache(semi.semi_gravity)

  e_potential = Trixi.integrate_via_indices(u_euler, mesh, equations_euler, dg, cache, equations_gravity, u_gravity) do u, i, j, element, equations_euler, dg, equations_gravity, u_gravity
    u_euler_local   = Trixi.get_node_vars(u_euler,   equations_euler,   dg, i, j, element)
    u_gravity_local = Trixi.get_node_vars(u_gravity, equations_gravity, dg, i, j, element)
    # OBS! subtraction is specific to Jeans instability test where rho0 = 1.5e7
    return (u_euler_local[1] - 1.5e7) * u_gravity_local[1]
  end
  return e_potential
end

analysis_callback = AnalysisCallback(semi_euler, interval=analysis_interval,
                                     save_analysis=true,
                                     extra_analysis_integrals=(energy_total, energy_kinetic, energy_internal, Val(:energy_potential)))

callbacks = CallbackSet(summary_callback, stepsize_callback,
                        save_restart, save_solution,
                        analysis_callback, alive_callback)

###############################################################################
# run the simulation
sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
            dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
            save_everystep=false, callback=callbacks);
summary_callback() # print the timer summary
println("Number of gravity subcycles: ", semi.gravity_counter.ncalls_since_readout)


text:
using CpelTdm, Test

@testset "Bioinformatics" begin
    # Aligning strand
    @test CpelTdm.get_align_strand(true,UInt16(99),UInt16(147)) == "OT"
    @test CpelTdm.get_align_strand(true,UInt16(83),UInt16(163)) == "OB"
    @test CpelTdm.get_align_strand(true,UInt16(147),UInt16(99)) == "CTOT"
    @test CpelTdm.get_align_strand(true,UInt16(163),UInt16(83)) == "CTOB"
    # Mean coverage
    xobs=[[1,-1] for i=1:10]; append!(xobs,[[1,0] for i=1:10]);
    @test CpelTdm.mean_cov(xobs)==15.0
    @test CpelTdm.get_obs_per_cpg(xobs)==[20,10]
end

@testset "Inference" begin
    # Energy function
    Ux_fun = CpelTdm.create_Ux([2,2],[1.0,1.0],1.0);
    @test Ux_fun([1,1,1,1]) == -7.0
    # Matrices
    @test isapprox(CpelTdm.get_W(4,0.0,0.0),[4.0 4.0; 4.0 4.0];atol=2)
    @test isapprox(CpelTdm.get_V([0.0,0.0],0.0),[1.0 1.0; 1.0 1.0];atol=2)
    @test CpelTdm.get_u(0.0) == [1.0,1.0]
    @test CpelTdm.comp_Z([4],[0.0],0.0) ≈ 16
    @test CpelTdm.comp_g([1],[0.0],0.0,0.0,0.0) == 2.0
    @test isapprox(CpelTdm.comp_lkhd([1,1,0,1,1],[5],[1.0],1.0),0.9536;atol=2)
    # MML
    n=[1,1,1]; θ=[1.0,-1.0,1.0,0.0]; ∇logZ = CpelTdm.get_∇logZ(n,θ);
    @test isapprox(CpelTdm.comp_mml(n,∇logZ),0.62693236)
    # NME
    n=[5]; α=[0.0]; β=0.0; θ=vcat([α,β]...); ∇logZ=CpelTdm.get_∇logZ(n,θ);
    @test isapprox(CpelTdm.comp_nme(n,θ,∇logZ),1.0)
    # CMD
    n=[20]; θ1=[-5.0,5.0]; θ2=[5.0,5.0];
    ∇logZ1s = CpelTdm.get_∇logZ(n,θ1); ∇logZ2s = CpelTdm.get_∇logZ(n,θ2);
    @test isapprox(CpelTdm.comp_cmd(n,θ1,θ2,∇logZ1s,∇logZ2s),1.0;atol=2)
end

@testset "Hypothesis Testing" begin
    # Unmatched test
    n=[10]; θ1s=fill([0.0,0.0],5); θ2s=fill([0.0,2.5],5);
    tmml_test,tnme_test,tpdm_test = CpelTdm.unmat_tests(n,θ1s,θ2s)
    @test tmml_test == (0.0, 1.0)
    @test tnme_test == (0.84782978, 0.007936507936507936)
    @test tpdm_test == (0.43786215630140846, 0.007936507936507936)
    # Matched test
    n=[10]; θ1s=fill([0.0,0.0],5); θ2s=fill([0.0,2.5],5);
    tmml_test,tnme_test,tpdm_test = CpelTdm.mat_tests(n,θ1s,θ2s)
    @test tmml_test == (0.0, 1.0)
    @test tnme_test == (0.84782978, 0.0625)
    @test tpdm_test[1] == 0.4378621563014087
end


text:
using Base: Integer

export RobotBasedProtocol, positions, postMoveWaitTime, numCooldowns, robotVelocity, switchBrakes, preMoveAction, postMoveAction

abstract type RobotBasedProtocol <: Protocol end
abstract type RobotBasedProtocolParams <: ProtocolParams end

positions(protocol::RobotBasedProtocol)::Union{Positions, Missing} = protocol.params.positions
postMoveWaitTime(protocol::RobotBasedProtocol)::typeof(1.0u"s") = protocol.params.postMoveWaitTime
numCooldowns(protocol::RobotBasedProtocol)::Integer = protocol.params.numCooldowns
robotVelocity(protocol::RobotBasedProtocol)::typeof(1.0u"m/s") = protocol.params.robotVelocity
switchBrakes(protocol::RobotBasedProtocol)::Bool = protocol.params.switchBrakes

function execute(protocol::RobotBasedProtocol)
  scanner_ = scanner(protocol)
  robot = getRobot(scanner_)

  positions_ = positions(protocol)
  vel = robotVelocity(protocol)
  switchBrakes_ = switchBrakes(protocol)

  enable(robot)
  if !isReferenced(robot)
    doReferenceDrive(robot)
  end
  movePark(robot)

  for (index, pos) in enumerate(positions_)
    for coord in pos
      @debug coord
    end
    # Cooldown pause
    numCooldowns_ = numCooldowns(protocol)
    if numCooldowns_ > 0 && index == round(Int, length(positions)/numCooldowns_)
      println("Cooled down? Enter \"yes\"")
      while readline() != "yes"
        println("Cooled down? Enter \"yes\"")
      end
    end

    @debug pos

    preMoveAction(protocol, upreferred.(pos))
    moveAbs(robot, upreferred.(pos), vel)
    sleep(ustrip(u"s", postMoveWaitTime(protocol)))

    # if hasBrake(robot) && switchBrakes_
    #   setBrake(robot, false)
    # end

    postMoveAction(protocol, pos)

    # if hasBrake(robot) && switchBrakes_
    #   setBrake(robot, true)
    # end
  end

  movePark(robot)
end

@mustimplement preMoveAction(protocol::RobotBasedProtocol)
@mustimplement postMoveAction(protocol::RobotBasedProtocol)

"Create the params struct from a dict. Typically called during scanner instantiation."
function createRobotBasedProtocolParams(ProtocolType::DataType, dict::Dict{String, Any})
  @debug "" ProtocolType dict
  for (key, value) in dict
    println("$key => $value")
  end
  @assert ProtocolType <: RobotBasedProtocolParams "The supplied type `$type` cannot be used for creating robot based protocol params, since it does not inherit from `ProtocolType`."

  # Split between main section fields and channels, which are dictionaries
  positionsDict = Dict{String, Any}()
  for (key, value) in dict["positions"]
    key = String([i == 1 ? uppercase(c) : c for (i, c) in enumerate(key)]) # No indexing to prevent possible errors with unicode
    positionsDict["positions"*key] = tryuparse.(value)
  end

  # Remove key in order to process the rest with the standard function
  delete!(dict, "positions")

  splattingDict = dict_to_splatting(dict)
  splattingDict[:positions] = Positions(positionsDict)

  try
    return ProtocolType(;splattingDict...)
  catch e
    if e isa UndefKeywordError
      throw(ScannerConfigurationError("The required field `$(e.var)` is missing in your configuration "*
                                      "for a device with the params type `$ProtocolType`."))
    elseif e isa MethodError
      @warn e.args e.world e.f
      throw(ScannerConfigurationError("A required field is missing in your configuration for a device "*
                                      "with the params type `$ProtocolType`. Please check "*
                                      "the causing stacktrace."))
    else
      rethrow()
    end
  end
end

include("RobotBasedMagneticFieldStaticProtocol.jl")
#include("RobotBasedMagneticFieldSweepProtocol.jl")
#include("RobotBasedSystemMatrixProtocol.jl")


text:
function dataexperiment(filling)
    ## internal tank ("sloshing related" parameters):
    # [pslosh pbeam prigid] = dataexperiment(filling)
	pslosh = Dict()
    pslosh["g"] = 9.8; # m/s^2
    pslosh["rho"] = 999; # water density: kg/m^3      #    +- 1
    L = 0.470;    # tank internal length (m)       #   +- 0.001
    R = 0.105/2; # tank internal radius (m)        #   +- 0.0002
    e = filling;
    #"equivalent" rectangular tank (same surface area and volume):
    theta = 2*acos(1-2*e);
    ls = 2*R*sqrt(1-(2*e-1)^2);
    Vcyl = R^2/2*(theta-sin(theta))*L;

    pslosh["a"] = L;           # rectangular tank length
    pslosh["b"] = ls;          # rectangular tank width
    pslosh["h"] = Vcyl/ls/L;   # rectangular tank liquid height

    pslosh["R"] = R;
    pslosh["H"] = e*2*R;
    ## structural parameters
    #plate parameters   - data in Ref. Robu2012,  'Simultaneous Hinf Vibration
    # Control of Fluid/Plate via Reduced-Order Controller')
    # --- IEEE Transactions on Control Sys Tech:
    E = 75e9; #Young modulus
    mu = 0.33;  #Poisson ratio
    rho = 2970; #density
    L = 1.36; #length
    h = 0.005;  #thickness       +- 0.00005
    l = 0.16;   #width           +- 0.001

    mrb = 0; # rigid body (tank) mass (here it is zero since we are adding
    # all rigid bodies when coupling with sloshing ("D" matrix))
    Irb = 0; # rigid body (tank) inertia (along axis x and y)
    xrb = 1.28; # position of rigid body (along axis x)        +- 0.001

	pbeam = Dict()
    pbeam["E"] = E; pbeam["mu"] = mu; pbeam["L"] = L; pbeam["rho"] = rho; pbeam["h"] = h;
    pbeam["l"] = l; pbeam["mrb"] = mrb; pbeam["Irb"] = Irb; pbeam["xrb"] = xrb;

    # piezoelectric data (data in ref. Robu 2012)
    rhop = 7800; # density
    Ep = 67e9; # young modulus
    tp = 5e-4; # thickness
    d31 = -2.1e-10; # piezoelectric constant

    pbeam["Ep"] = Ep; pbeam["tp"] = tp; pbeam["d31"] = d31; pbeam["rhop"] = rhop;

    ## rigid mass and inertias
    #(tank, fixation rings and fluid inertia in "y" direction)
    #tank:

    # tank data - data in Robu2012/IEEE
    rhoplast = 1180; # plastic density (kg/m^3)
    Der = 0.11;      # tank external diameter (m)
    Dir = 0.105;     # tank internal diameter (m)
    Lr = 0.50;       # tank length (m)              +- 0.001

    # fixation rings data - data measured/approximated
    # obs.: external diameter and thickness doesn't represent the exact
    # values. actually, the fixation rings have several holes/screws.
    # Def, Di and hf are approximate values. They were adjusted so that
    # the TOTAL fixation ring masses are the measured ones
    #   (588.6g)
    rhofix = 2970;       # fixation rings density (= plate,kg/m^3)
    Def = 0.1445 + 0.15*0; # fixation rings external diameter (m)
    Dif = 0.11;            # fixation rings internal diameter (m)
    hf = 0.02876+0.0325*0; # fixation rings thickness (m)
    distf = (h+hf)/2;      # fixation rings position (relative to plate mean axis) (m)

    Dh = Der;        # PLATE hole diameter (m)
    massrlat = rhoplast* pi * (Der^2-Dir^2)/4 * Lr; # tank lateral mass
    massrtip = 0.313/2;     # each tip mass (measured)
                            # OBS.: the tip mass was obtained by measuring
                            # the tank mass and comparing with the predicted
                            # lateral tank mass. In order to confirm if the
                            # lateral tank mass was correct, I also
                            # compared the results with a 30cm tank (in
                            # addition to the 50cm one).
                            # The tip consists of: plastic disc + glue +
                            # plug: not so simple geometry, but the effects
                            # on the system mass/inertia can be
                            # approximated as a point mass.
    massfix = rhofix*pi*(Def^2-Dif^2)/4 * hf;         # each fixation rings total
    masshole = rho*pi*Dh^2/4 * h;                     # plate hole mass - it will reduce the rigid unit mass!
    mass = (massrlat + massrtip*2+2*massfix-masshole);# total rigid added mass (kg)

    # calculate rigid body inertia (x and y direction, since tank and
    # fixation rings are symmetric):
    Inertialat = (massrlat/12*(3*(Der^2+Dir^2)/4+Lr^2)); # tank lateral inertia
    Inertiatip = massrtip*(Lr/2-0.01)^2; # tank tip inertia
    Inertiafix = massfix*(3*(Dif^2+Def^2)/4+hf^2)/12+massfix*distf^2; # fixation rings inertia
    Inertiahole = masshole*(Dh^2)/16;   # inertia of hole in the plate
    Inertia = (Inertialat + 2*Inertiatip + 2*Inertiafix - Inertiahole); #(kg.m^2)

    # "frozen" fluid inertia (y direction)
    massfluid = pslosh["a"]*pslosh["b"]*pslosh["h"]*pslosh["rho"];
    IFLUID = massfluid*pslosh["a"]^2/12;    #(kg.m^2)

	prigid = Dict()
    prigid["xrb"] = xrb; #position of rigid units (m)
    prigid["mass"] = mass; #tank + fixation rings mass (kg)
    prigid["Inertia"] = Inertia; #tank + fixation rings inertia ("y" and "z" directions) (kg.m^2)
    prigid["Ifluid"] = IFLUID; # "frozen" fluid inertia in "y" direction (kg.m^2)
    prigid["massfluid"] = massfluid; # "frozen" fluid mass

    ## auxiliary equations and parameters: most equations are implemented on
    # other codes
    # they are repeated here only for verification purposes

    # beam structural parameters
    J = 0.33*l*h^3;              # torsion constant
    IY = rho/12*(l^2 + h^2)*l*h; # section moment of inertia (X direction)  # VERIFICAR!!!
    Iz = l*h^3/12;               # section area inertia (X direction)
    Ix = l^3*h/12;               # section area inertia (X direction)
    G = E/2/(1+mu);              # shear modulus

   return (pslosh, pbeam, prigid)
end

## measures:
# large tank: 50 cm  (47 cm liquid)
# tank filled up to 34.3 cm (+-0.2):
#    mass = 3820.3 g +- 0.1 g
# tank filled 100 #:
#    mass = 4938.0 g +- 0.4 g
# tank filled 0#:
#    mass = 821.8 g +- 0.1 g
#
# small tank: 30 cm (27 cm liquid)
#  empty tank:
#    606.7 g +- 0.1 g
#
# small tank: 30 cm (27 cm liquid): broken one
#    578.1 g +- 0.1 g
#

## new measures:
# TANK WEIGHT:
#   * empty: 818.2 gramas (erro do aparelho: 0.1, mas gotas extras
#   podem adicionar uns 0.2 g)
#   * 0 +-0.1cm:      839.8  +-0.1g
#   * 20.2 +- 0.2cm: 2551.2  +-0.2g
#   * 35.0 +- 0.1cm: 3879.2  +-0.2g
#   * 47.0 +- 0.1cm: 4933.4  +-0.2g
#   * full: 4939.9 +- 0.1 g (mas nao estava 100# cheio, impossivel
#   fazer isso)
#
# TANK medidas dos bouchons individuais: 15.9, 15.8, 16.7, 16.2
# (gramas)
#
# TANK dimensoes:
# L = 30.0 +-0.1 cm
# Li = 27.0 +- 0.1cm
# espessura do disco interno de plastico: 0.70 +-0.05 cm
# distancia entre a parte interna do disco e o fim do bouchon: 4.5cm
# distancia entre a parte interna do disco e a extremidade: 1.5cm
#
#
# ANEIS de fixacao:
#  1) lado direito (visto de frente)
#       561.4 g (sem parafusos)
#       568.2 g (parafuso central)
#       580.6 g (parafuso central+interm)
#       596.2 g (parafuso central+interm+sup)
#   2) lado esquerdo
#       565.3 g (sem parafusos)
#       580.9 g (parafuso superior)
#
#
# peso dos parafusos do anel:
#   central: 3.3g
#   interm: 6.2g
#   sup: 15.7g
#
# dimensoes do anel:
#  diametro exterior 1: 16.0cm, 16.0cm, 16.0cm +- 0.1
#  diametro exterior 2 (anel 1): 145.46mm, 145.52mm,145.45mm,145.40mm
#  diametro exterior 2 (anel 2): 144.62mm, 144.56mm,144.51mm
#  diametro interior (anel 1): 110.53mm, 110.8mm,110.81mm,110.72mm
#  diametro interior (anel 2): 110.73mm, 110.90mm,110.62mm
#  altura placa fina: 5.16 5.22 5.02 4.97 4.48mm
#  altura total anel: 30.15 30.01 30.13 29.79mm
#
#
# DIMENSOES DA PLACA
#    comprimento: 136.3 +-0.2 cm
#    largura: 16.0 +- 0.1 cm
#    espessura: 4.95mm, 4.98mm, 4.92mm, 4.92mm, 4.84mm, 5.02mm, 4.88mm,
#    4.86mm
#
#
# POSICAO DO FURO NA PLACA:
#  relativo a ponta da placa:
#    entre 2.4cm e 13.65cm
#  diametro interno do furo:
#          112.05 +-0.3mm

max_stars_repo_name: flavioluiz/phjulia

max_issues_repo_name: flavioluiz/phjulia

max_forks_repo_name: flavioluiz/phjulia


text:
const Vec{Dim,T} = SVector{Dim,T}
const Vec3 = Vec{3,Float}
const Vec2 = Vec{2,Float}
const Vec1 =  Vec{1,Float}