text:
using DataFrames, Gadfly

function plot_traj_with_steps(traj,grid_path;Δt=0.01,pad=0.01)
    t_range = get_start_time(traj):Δt:get_end_time(traj)
    cw = grid_path.cellwidth/2 - pad
    xpts = [map(w->w.pt.x, grid_path.waypoints)..., grid_path.start_pt.x]
    ypts = [map(w->w.pt.y, grid_path.waypoints)..., grid_path.start_pt.y]
    plot(
        layer( x=[get_position(traj,t).x for t in t_range], y=[get_position(traj,t).y for t in t_range],
            Geom.path, Theme(default_color="black") ),
        layer( x=[get_start_pt(seg).x for seg in traj.segments], y=[get_start_pt(seg).y for seg in traj.segments],
            size=[8pt], Geom.point, Theme(default_color="red") ),
        layer( x=[get_end_pt(seg).x for seg in traj.segments], y=[get_end_pt(seg).y for seg in traj.segments],
            size=[8pt],Geom.point, Theme(default_color="yellow") ),
        layer( xmin=xpts .- cw, xmax=xpts .+ cw, ymin=ypts .- cw, ymax=ypts .+ cw, Geom.rect ),
        Coord.cartesian(fixed=true)
    )
end
function compare_velocity_profiles(traj1,traj2;Δt=0.05)
    t0 = max(get_start_time(traj1),get_start_time(traj2))
    tf = min(get_end_time(traj1),get_end_time(traj2))
    @assert t0 < tf "trajectories do not overlap"
    t_vec = collect(get_start_time(base_traj):Δt:get_end_time(base_traj))
    df = DataFrame(
        t=t_vec,
        v1=[norm(get_vel(traj1,t)) for t in t_vec],
        v2=[norm(get_vel(traj2,t)) for t in t_vec],
    )
    plot(
        layer(df,x="t",y="v2",Geom.path,Theme(default_color="red")),
        layer(df,x="t",y="v1",Geom.path,Theme(default_color="black")),
        Guide.manual_color_key("",["v1","v2"],["black","red"])
    )
end
function plot_vel_and_yaw_rate(traj;Δt=0.1)
    t_vec = collect(get_start_time(traj):Δt:get_end_time(traj))
    df = DataFrame(
        t=t_vec,
        theta=[atan(get_heading(traj,t)) for t in t_vec],
        s=[get_dist(traj,t) for t in t_vec],
        v=[norm(get_vel(traj,t)) for t in t_vec] )
    vstack( plot(df,x="t", y="theta",Geom.path),
        plot(df,x="t",y="v",Geom.path,Theme(default_color="red")) )
end
function plot_accel_vel_pos(base_traj,traj,t_vec,accel,vel,pos)
    df = DataFrame(
        t = t_vec,
        a = [accel..., 0.0],
        v = vel,
        v1=map(t->norm(get_vel(base_traj,t)),t_vec),
        s = pos
    )
    df2 = DataFrame(
        t = [[get_start_time(seg) for seg in traj.segments]..., get_end_time(traj)],
        s = [0, cumsum([get_length(seg) for seg in traj.segments])...]
    )
#     vstack(
        plot(
            layer(df,x="t",y="a",Geom.step,Theme(default_color="green")),
            layer(df,x="t",y="v",Geom.path,Theme(default_color="red")),
            layer(df,x="t",y="v1",Geom.path,Theme(default_color="black")),
            Guide.manual_color_key("",["accel","vel_opt","vel"],["green","red","black"])
        )
#         plot(
#             layer(df,x="t",y="s",Geom.path,Theme(default_color="blue")),
#             layer(df2,x="t",y="s",Geom.point,Theme(default_color="orange")),
#             Guide.manual_color_key("",["distance","constraints"],["blue","orange"])
#         )
#     )
end
function summarize_simulation(traj,states,cmds,time_vec)
    if length(time_vec) < length(states)
        time_vec = [time_vec..., 2*time_vec[end]-time_vec[end-1]]
    end
    errors = Vector{Vector{Float64}}()
    for (t,state) in zip(time_vec,states)
        target_pt = get_trajectory_point_by_time(traj, t)
        target = [target_pt.pos.x, target_pt.pos.y, atan(target_pt.heading)]
        push!(errors, [target[1]-state[1],target[2]-state[2],get_angular_offset(state[3],target[3])])
    end
    df = DataFrame(
        t=time_vec,
        w=[map(u->u[1],cmds)...,0.0],
        v=[map(u->u[2],cmds)...,0.0],
        wr=map(t->get_yaw_rate(traj,t),time_vec),
        vr=map(t->norm(get_vel(traj,t)),time_vec),
        x=map(s->s[1],states[1:length(time_vec)]),
        y=map(s->s[2],states[1:length(time_vec)]),
        θ=map(s->wrap_to_pi(s[3]),states[1:length(time_vec)]),
        xr=map(t->get_position(traj,t).x,time_vec),
        yr=map(t->get_position(traj,t).y,time_vec),
        θr=map(t->atan(get_heading(traj,t)),time_vec),
        ex=map(e->e[1],errors),
        ey=map(e->e[2],errors),
        eθ=map(e->e[3],errors)
    )
    hstack(
            vstack(
                plot(
                    layer(df,x="x",y="y",Geom.path,Theme(default_color="blue")),
                    layer(x=[states[1][1]],y=[states[1][2]],Geom.point,Theme(default_color="blue")),
                    layer(df,x="xr",y="yr",Geom.path,Theme(default_color="red")),
                    layer(x=[get_position(traj,time_vec[1]).x],
                        y=[get_position(traj,time_vec[1]).y],Geom.point,Theme(default_color="red")),
                    Coord.Cartesian(fixed=true),
                    Guide.manual_color_key("",["true","ref"],["blue","red"])
                )
            ),
            vstack(
                plot(
                    layer(df,x="t",y="w",Geom.path,Theme(default_color="red")),
                    layer(df,x="t",y="wr",Geom.path,Theme(default_color="black")),
                    Guide.manual_color_key("",["w","wref"],["red","black"]) ),
                plot(
                    layer(df,x="t",y="v",Geom.path,Theme(default_color="blue")),
                    layer(df,x="t",y="vr",Geom.path,Theme(default_color="black")),
                    Guide.manual_color_key("",["v","vref"],["blue","black"]) ),
                plot(
                    layer(df,x="t",y="ex",Geom.path,Theme(default_color="blue")),
                    layer(df,x="t",y="ey",Geom.path,Theme(default_color="red")),
                    layer(df,x="t",y="eθ",Geom.path,Theme(default_color="green")),
                    Guide.manual_color_key("",["ex","ey","eθ"],["blue","red","green"]) )
            )
        )
end


text:
module Ex5_17
using FaultDetectionTools
using DescriptorSystems
using LinearAlgebra
using Polynomials
using Test

# Example 5.17 - Solution of an H2 AMMP
println("Example 5.17")

# define system with control, noise and actuator fault inputs
A = [-.8 0 0; 0 -.5 .6; 0 -.6 -.5];
Bu = [1 1;1 0;0 1]; Bw = 0.25*[0 0;0 1;1 0]; Bf = Bu;
C = [0 1 1; 1 1 0]; Du = zeros(2,2); Dw = zeros(2,2);
p, mu = size(Du); mf = mu; mw = size(Dw,2);

# setup synthesis model with additive actuator faults with Bf = Bu, Df = Du
sysf = fdimodset(dss(A,[Bu Bw],C,[Du Dw]), c = 1:mu, n = mu .+ (1:mw),f = 1:mu);

# define Mr(s) = 10/(s+10)
s = rtf("s")
Mr = FDFilterIF(dss([10/(s+10) 0; 0 10/(s+10)]),0,0,mf)

Q, R, info = ammsyn(sysf, Mr; H2syn = true, atol = 1.e-7, nullspace = false, reltol = 5.e-4, sdeg = -10, normalize = "dcgain");

println("gamma_opt0 = $(info.gammaopt0)")  # optimal performance for the initial problem
println("gamma_opt  = $(info.gammaopt)")   # optimal performance
println("gamma_sub  = $(info.gammasub)")   # suboptimal performance

# check suboptimal solution
Ge = [sysf.sys[:,[sysf.faults;sysf.noise]]; zeros(mu,mf+mw)]; Me = [info.M*Mr.sys zeros(mf,mw)];
@test iszero(R.sys-Q.sys*Ge, atol = 1.e-7) &&
      fdimmperf(R, info.M*Mr,2; atolinf = 1.e-7) ≈ info.gammasub && gl2norm(Me-Q.sys*Ge)[1] ≈ info.gammasub

end # module

ext: jl

max_stars_repo_name: andreasvarga/FaultDetectionTools

max_issues_repo_name: andreasvarga/FaultDetectionTools

max_forks_repo_name: andreasvarga/FaultDetectionTools


text:
module SeaPearl

using Random
using ReinforcementLearning
const RL = ReinforcementLearning

include("abstract_types.jl")

include("trailer.jl")
include("CP/CP.jl")
#include("MOI_wrapper/MOI_wrapper.jl")
include("datagen/datagen.jl")
include("experiment/experiment.jl")

greet() = print("Hello World!")

end # module


text:
using Test
using Permutations

@testset "Permutation" begin
    p = Permutation(7:-1:1)
    i = Permutation(7)
    @test p*p == i
    @test p == inv(p)
    a =[2,3,1,6,7,8,5,4,9]
    p = Permutation(a)
    @test order(p) == 6
    @test p^6 == Permutation(9)
    @test p^-5 == p
    @test p.data == a

    q = Permutation([1,5,3,9,4,8,6,7,2])
    @test p*q != q*p
    @test q[2] == 5
    @test q(2) == 5

    P = Matrix(p)
    Q = Matrix(q)
    @test P*Q == Matrix(p*q)
    @test P' == Matrix(inv(p))

    row = [1,3,5,2,4,6]
    p = Permutation(row)
    @test length(p) == 6
    @test p[1] == 1
    @test length(longest_increasing(p)) == 4
    @test length(longest_decreasing(p)) == 2
    @test length(fixed_points(p)) == 2
    @test sign(p) == -1
    @test hash(p) == hash(p.data)

    M = two_row(p)
    @test M[2,:] == row

    @test Permutation(6,1) == Permutation(6)

    p = RandomPermutation(10)
    @test p*inv(p) == Permutation(10)
    @test p'*p == Permutation(10)
    @test reverse(reverse(p)) == p

    a,b,c,d = Permutation([2,1,3,4]),Permutation([1,2,4,3]),
                Permutation([1,3,2,4]),Permutation([2,1,3,4])

    @test Matrix(a*b*c*d) == Matrix(a)*Matrix(b)*Matrix(c)*Matrix(d)

    p = RandomPermutation(12)
    d = dict(p)
    for k=1:12
        @test p(k)==d[k]
    end

    p = RandomPermutation(12)
    pp = p*p
    q = sqrt(pp)
    @test q*q == pp

end

@testset "CoxeterDecomposition" begin
    n = 10
    p = RandomPermutation(n)
    c = CoxeterDecomposition(p)
    @test p == Permutation(c)

    #s₂ = CoxeterGenerator(n, 2)
    #@test Permutation(s₂) == Permutation([1; 3; 2; 4:n])
    #@test s₂*c isa CoxeterDecomposition
    #@test c*s₂ isa CoxeterDecomposition
    #@test s₂*s₂ isa CoxeterDecomposition
    #@test s₂*c == CoxeterDecomposition(s₂)*c
    #@test Permutation(s₂*c) == Permutation(CoxeterDecomposition(s₂)*c) == Permutation(s₂)*p

    @test inv(c)*c == CoxeterDecomposition(n, Int[])
    @test inv(c) == c'

    @test CoxeterDecomposition(5, [3,4,1]) == CoxeterDecomposition(5, [1,3,4])
    @test CoxeterDecomposition(5, [2,1,3,4,1]) == CoxeterDecomposition(5, [2,3,4])
    @test CoxeterDecomposition(5, [1,3,4,3,4,3,4,1]) == CoxeterDecomposition(5, Int[])

    @test Permutation(CoxeterDecomposition(5, Int[])).data ≈ 1:5
    @test Permutation(CoxeterDecomposition(5, [1])).data == [2; 1; 3:5]
    @test inv(CoxeterDecomposition(6, [1,3,5])) == CoxeterDecomposition(6, [1,3,5])
end


text:
set_theme!(
    rowgap = 5,
    colgap = 5,
    fontsize = 16,
    backgroundcolor = :black, # "#212946",
    Axis = (
        xgridcolor = :white,
        ygridcolor = :white,
        backgroundcolor = :black, # "#212946",
        xgridstyle=:dash,
        ygridstyle=:dash,
        xgridwidth=0.35,
        ygridwidth=0.35,
        xtickalign=0,
        ytickalign=0,
        xtickcolor = :white,
        ytickcolor = :white,
        bottomspinecolor = :white,
        topspinecolor =:white,
        leftspinecolor =:white,
        rightspinecolor =:white,
        xlabelcolor=:white,
        ylabelcolor = :white,
        xticklabelcolor = :white,
        yticklabelcolor = :white,
        xgridvisible = true,
        ygridvisible = true),
    Colorbar = (
        topspinevisible = false,
        rightspinevisible = false,
        bottomspinevisible = false,
        leftspinevisible = false,
        width = 12,
        height = Relative(4/4),
        labelcolor = :white,
        ticklabelcolor = :white,
        tickalign = 1,
        tickcolor = :white,
        labelpadding = -5,
        ticklabelsize = 16),
    Legend = (#titlecolor = :white # not working or missing
        tellheight = false,
        tellwidth = false,
        halign = :right,
        valign = :top,
        labelsize = 16,
        linewidth = 2,
        markersize = 7,
        margin = (10, 10, 10, 10),
        framevisible = true,
        framecolor = (:white, 0.25),
        bgcolor = (:grey, 0.25),
        labelcolor = :white,
        titlecolor = :white)
)

max_stars_repo_name: aivodji/MakieNotebooks

max_issues_repo_name: aivodji/MakieNotebooks

max_forks_repo_name: aivodji/MakieNotebooks


text:
using CSV
using MultiChannelFactorAnalysis
using DataFrames
using Statistics
using Distributions
using LinearAlgebra
#using Arpack
using ArnoldiMethod

function fitlargeFA(data, Ls, p, pjs)
    X, Xobs, Xchan = MultiChannelFactorAnalysis.stack_and_view(data)
    k = p+sum(pjs)
    K = maximum(size(X))
    decomp, history = partialschur(transpose(X) * X; nev=k, which=LM())
    ls, evecs = partialeigen(decomp)
    evecs = reverse(evecs, dims=2)
    Fhat = sqrt(size(X)[2]) * evecs
    Chat = X * Fhat / size(X)[2]
    H = Chat[:, 1:p]
    G = Chat[:, (p+1):end]
    #Gs = MultiChannelFactorAnalysis.extract_blocks(Chat[:, (p+1):end] , Ls, pjs;diag=true)
    #H = MultiChannelFactorAnalysis.flip_signs(H * MultiChannelFactorAnalysis.LTOrthog(H, p))
    #Gs = [MultiChannelFactorAnalysis.flip_signs(G * MultiChannelFactorAnalysis.LTOrthog(G, size(G)[2])) for G in Gs]
    #G = MultiChannelFactorAnalysis.blockdiag(Gs)

    return H, G
end

function fitIGFM(data, Ls, p, pjs, delta)
    H, G = fitlargeFA(data, Ls, p , pjs)
    N = sum(Ls)
    T = length(data)
    m = length(Ls)
    #delta * log(1/mest * 1/log(p)* time )*(1/sqrt(time) * sqrt(log(p)) + 1/sqrt(p))
    τ = delta * log(1/m *1/log(N)  * T ) * (1/sqrt(T) * sqrt(log(N)) + 1/sqrt(N))
    Gthresh = G
    Gthresh[abs.(Gthresh) .< τ] .= 0
    rbounds = vcat([0], cumsum(Ls))
    cbounds = vcat([0], cumsum(pjs))

    #block_sums = [[ sum(abs.(Gthresh[rstart:rstop, cstart:cstop])) for  (cstart, cstop) in zip(cbounds .+ 1, cbounds[2:end])]
    #        for (rstart, rstop) in zip(rbounds .+ 1, rbounds[2:end])]

    #permutes = [argmax(block_sums[i, :][1]) for i in 1:size(block_sums)[1]]
    #perm_Gthresh = blockdiag()
    return H, Gthresh
end
Ls = [120, 120, 120]
p = 1
pjs = [1, 1, 1]
N=450

function loading_trial(N, Ls, p, pjs; factor_func=MultiChannelFactorAnalysis.simFactor)
    Hs, Gs, Sigs = MultiChannelFactorAnalysis.powerStructLoadings(Ls, p, pjs)
    H = vcat(Hs...)
    G = MultiChannelFactorAnalysis.blockdiag(Gs)
    edists = [ MvNormal(Sig) for (L, Sig) in zip(Ls, Sigs)]
    data = MultiChannelFactorAnalysis.simData(N, Hs, Gs, edists)
    H1, G1, Sig1 = MultiChannelFactorAnalysis.fitMFA(data, Ls, p, pjs; maxiter=20)
    #G1s = MultiChannelFactorAnalysis.extract_blocks(G1, Ls, pjs; diag=true)
    #H1 = MultiChannelFactorAnalysis.flip_signs(H1)
    #G1s = [MultiChannelFactorAnalysis.flip_signs(G * MultiChannelFactorAnalysis.LTOrthog(G, size(G)[2])) for G in G1s]
    #G1 = MultiChannelFactorAnalysis.blockdiag(G1s)
    H2, G2 = fitlargeFA(data, Ls, p, pjs)
    R = H * transpose(H) + G * transpose(G)
    R1 = H1 * transpose(H1) + G1 * transpose(G1)
    R2 = H2 * transpose(H2) + G2 * transpose(G2)
    return norm(R - R1)^2/norm(R)^2, norm(R - R2)^2/norm(R)^2
end

function get_trial_results(times, d)
    Ls = [d, floor(Int32,d/2), floor(Int32, d/4)]
    p = 2
    pjs = repeat([1], 3)
    log1s = []
    log2s = []
    for N in 50:50:450
        print(N)
        r1s = []
        r2s = []
        for i in 1:times
            r1, r2 = loading_trial(N, Ls, p, pjs)
            push!(r1s, log(r1))
            push!(r2s, log(r2))
        end
        push!(log1s, mean(r1s))
        push!(log2s, mean(r2s))
    end
    return (log1s, log2s)
end

ds = [10, 20, 40, 60, 100]
fixed_p_loss3 = []
fixed_p_loss4 = []
for d in ds
    print(d)
    v1, v2 = get_trial_results(times, d)
    push!(fixed_p_loss3, v1)
    push!(fixed_p_loss4, v2)
end

df1 = DataFrame(MCFMlosses = vcat(fixed_p_loss3...),
                IGFMlosses = vcat(fixed_p_loss4...),
                d    =  repeat(ds, inner=[9]),
                N          = repeat(collect(50:50:450), outer=[5])
)

CSV.write("/home/gray/research/ngsim/FactorComparison_pfixed.csv", df1)
3
pnratios = [1/15, 2/15, 3/15, 4/15, 5/15, 10/15, 15/15]
times = 10
fixed_p_loss1 = []
fixed_p_loss2 = []
for rat in pnratios
    print(rat)
    s1, s2 = get_trial_results2(times, rat)
    push!(fixed_p_loss1, s1)
    push!(fixed_p_loss1, s2)
end

df1 = DataFrame(MCFMlosses = vcat(fixed_p_loss1[1], fixed_p_loss1[3], fixed_p_loss1[5], fixed_p_loss1[7], fixed_p_loss1[9]),
                IGFMlosses = vcat(fixed_p_loss1[2], fixed_p_loss1[4], fixed_p_loss1[6],fixed_p_loss1[8], fixed_p_loss1[10]),
                pnratio    = repeat([1/15, 2/15, 3/15, 4/15, 5/15], inner=[9]),
                N          = repeat(collect(50:50:450), outer=[5])
)

CSV.write("/home/gray/research/ngsim/FactorComparison_pn_grow.csv", df1)

function get_trial_results2(times, pnratio)
    p = 1
    pjs = repeat([1], 3)
    log1s = []
    log2s = []
    for N in 50:50:450
        print(N)
        d = floor(Int32, N * pnratio)
        Ls = repeat([d], 3)
        r1s = []
        r2s = []
        for i in 1:times
            r1, r2 = loading_trial2(N, Ls, p, pjs)
            push!(r1s, log(r1))
            push!(r2s, log(r2))
        end
        push!(log1s, mean(r1s))
        push!(log2s, mean(r2s))
    end
    return (log1s, log2s)
end

function get_trial_results3(times, d)
    p = 1
    pjs = repeat([1], 3)
    log1s = []
    log2s = []
    for N in 50:50:450
        print(N)
        Ls = repeat([d], 3)
        r1s = []
        r2s = []
        for i in 1:times
            r1, r2 = loading_trial2(N, Ls, p, pjs)
            push!(r1s, log(r1))
            push!(r2s, log(r2))
        end
        push!(log1s, mean(r1s))
        push!(log2s, mean(r2s))
    end
    return (log1s, log2s)
end

function loading_trial2(N, Ls, p, pjs; factor_func=MultiChannelFactorAnalysis.simFactor)
    Hs, Gs, Sigs = MultiChannelFactorAnalysis.LyuouLoadings1(Ls, p, pjs)
    H = vcat(Hs...)
    G = MultiChannelFactorAnalysis.blockdiag(Gs)
    edists = [ MvNormal(Sig) for (L, Sig) in zip(Ls, Sigs)]
    data = MultiChannelFactorAnalysis.simData(N, Hs, Gs, edists)
    H1, G1, Sig1 = MultiChannelFactorAnalysis.fitMFA(data, Ls, p, pjs; maxiter=15)
    #G1s = MultiChannelFactorAnalysis.extract_blocks(G1, Ls, pjs; diag=true)
    #H1 = MultiChannelFactorAnalysis.flip_signs(H1)
    #G1s = [MultiChannelFactorAnalysis.flip_signs(G * MultiChannelFactorAnalysis.LTOrthog(G, size(G)[2])) for G in G1s]
    #G1 = MultiChannelFactorAnalysis.blockdiag(G1s)
    H2, G2 = fitIGFM(data, Ls, p, pjs, 1)
    R = H * transpose(H) + G * transpose(G)
    R1 = H1 * transpose(H1) + G1 * transpose(G1)
    R2 = H2 * transpose(H2) + G2 * transpose(G2)
    return norm(R - R1)^2/norm(R)^2, norm(R - R2)^2/norm(R)^2
end


text:
# types.jl - common types and aliases used throught the code

# name of operation in a :call node - either symbol or Module.symbol
const OpName = Union{Symbol, Expr}

# Wrapper around Expr adding `ex.head` to type parameters thus adding convenient
# type dispatching
mutable struct ExH{H}
    head::Symbol
    args::Vector
end

Base.show(io::IO, ex::ExH) = print(io, "ExH{:$(ex.head)}($(Expr(ex)))")

Base.convert(::Type{ExH}, ex::Expr) = ExH{ex.head}(ex.head, ex.args)
ExH(ex::Expr) = ExH{ex.head}(ex.head, ex.args)
Expr(ex::ExH) = Expr(ex.head, ex.args...)


text:
struct SolverState{TF<:AbstractFloat}
    iter::Int
    f::TF
    rss::TF
    xnorm::TF
    idx::Int
    step::TF
end

struct SolverTrace{TF<:AbstractFloat}
    states::Vector{SolverState{TF}}
end

@fieldequal SolverTrace

SolverTrace(TF::Type) = SolverTrace(Vector{SolverState{TF}}(undef, 0))

function show(io::IO, st::SolverState, show_thread::Bool=false)
    if show_thread
        @printf io " thread %2d:  iter=%-8d  f=%-8.3g  rss=%-8.3g  xnorm=%-8.3g  idx=%-8d  step=%-8.3g\n" Threads.threadid() st.iter st.f st.rss st.xnorm st.idx st.step
    else
        @printf io " iter=%-8d  f=%-8.3g  rss=%-8.3g  xnorm=%-8.3g  idx=%-8d  step=%-8.3g\n" st.iter st.f st.rss st.xnorm st.idx st.step
    end
end

push!(tr::SolverTrace, st::SolverState) = push!(tr.states, st)

getindex(tr::SolverTrace, i::Integer) = getindex(tr.states, i)

function show(io::IO, tr::SolverTrace)
    for state in tr.states
        print(io, state)
    end
end

function update!(tr::SolverTrace, iter, f, rss, xnorm, idx, step,
        store_trace::Bool, show_trace::Bool, show_thread::Bool)
    st = SolverState(iter, f, rss, xnorm, idx, step)
    store_trace && push!(tr, st)
    show_trace && show(stdout, st, show_thread)
    return tr
end

update!(::Nothing, args...) = nothing

struct SolverResult{TF<:AbstractFloat}
    x0::Vector{TF}
    sol::Vector{TF}
    fit::Vector{TF}
    iter::Int
    f::TF
    rss::TF
    λ::TF
    xnorm::TF
    dfnorm::TF
    ftol::TF
    f_converged::Bool
    dxnorm::TF
    xtol::TF
    x_converged::Bool
    trace::Union{SolverTrace{TF},Nothing}
end

@fieldequal SolverResult

converged(r::SolverResult) = r.f_converged || r.x_converged

function show(io::IO, r::SolverResult)
    print(io, typeof(r).name.name, "(converged=", converged(r), ", iter=$(r.iter)")
    print(IOContext(io, :compact=>true), ", rss=$(r.rss), λ=$(r.λ), xnorm=$(r.xnorm))")
end

function show(io::IO, ::MIME"text/plain", r::SolverResult)
    print(io, "convexfit ", converged(r) ? "converged " : "did not converge ")
    print(io, "after ", r.iter, " iteration")
    r.iter > 1 && print(io, "s")
    @printf io ":\n |f(x) - f(x')| = %.1e < %.1e: %s\n" r.dfnorm r.ftol r.f_converged
    @printf io " |x - x'| = %.1e < %.1e: %s\n" r.dxnorm r.xtol r.x_converged
    @printf io " f=%-8.3g  rss=%-8.3g  λ=%-8.3g  xnorm=%-8.3g" r.f r.rss r.λ r.xnorm
end


text:
function apply_boundary_conditions_on_φ_axis!(  rbpot::AbstractArray{T, 4}, rbi::Int,
                                                ax::DiscreteAxis{T, :periodic, :periodic}, int::Interval{:closed, :open, T} )::Nothing where {T}
    @inbounds rbpot[:,   1, :, rbi] .= view(rbpot, :, size(rbpot, 2) - 1, :, rbi) # cycling boundary
    @inbounds rbpot[:, end, :, rbi] .= view(rbpot, :,                  2, :, rbi) # cycling boundary
    nothing
end

function apply_boundary_conditions_on_φ_axis!(  rbpot::AbstractArray{T, 4}, rbi::Int,
                                                ax::DiscreteAxis{T, :reflecting, :reflecting}, int::Interval{:closed, :closed, T} )::Nothing where {T}
    @inbounds rbpot[:,   1, :, rbi] .= view(rbpot, :,                  3, :, rbi) # cycling boundary
    @inbounds rbpot[:, end, :, rbi] .= view(rbpot, :, size(rbpot, 2) - 2, :, rbi) # cycling boundary
    nothing
end

function apply_boundary_conditions_on_r_axis!(  rbpot::AbstractArray{T, 4}, rbi::Int, ax::DiscreteAxis{T, :r0, :infinite}, int::Interval{:closed, :closed, T},
                                                grid_boundary_factors::NTuple{2, T})::Nothing where {T}
    @inbounds rbpot[:, :, end, rbi] .= grid_boundary_factors[2] .* view(rbpot, :, :, size(rbpot, 3) - 2, rbi) # infinity boundaries in r
    nothing
end
function apply_boundary_conditions_on_r_axis!(  rbpot::AbstractArray{T, 4}, rbi::Int, ax::DiscreteAxis{T, :r0, :reflecting}, int::Interval{:closed, :closed, T},
                                                grid_boundary_factors::NTuple{2, T})::Nothing where {T}
    @inbounds rbpot[:, :, end, rbi] .= view(rbpot, :, :, size(rbpot, 3) - 2, rbi) # infinity boundaries in r
    nothing
end
function apply_boundary_conditions_on_r_axis!(  rbpot::AbstractArray{T, 4}, rbi::Int, ax::DiscreteAxis{T, :r0, :fixed}, int::Interval{:closed, :closed, T},
                                                grid_boundary_factors::NTuple{2, T})::Nothing where {T}
    nothing
end

function apply_boundary_conditions!(pssrb::PotentialSimulationSetupRB{T, Cylindrical}, update_even_points::Val{even_points}, only2d::Val{only_2d}) where {T, even_points, only_2d}
    rbi::Int = even_points ? rb_even::Int : rb_odd::Int
    # Radial Axis
    apply_boundary_conditions_on_r_axis!( pssrb.potential, rbi, pssrb.grid.axes[1], pssrb.grid.axes[1].interval, pssrb.grid_boundary_factors[1])
    # # Azimutal-Axis
    if !only_2d
        apply_boundary_conditions_on_φ_axis!( pssrb.potential, rbi, pssrb.grid.axes[2], pssrb.grid.axes[2].interval)
    end
    # Cylindrical Z-Axis -> same as Cartesian X-Axis
    apply_boundary_conditions_on_x_axis!( pssrb.potential, rbi, pssrb.grid.axes[3], pssrb.grid.axes[3].interval, pssrb.grid_boundary_factors[3])
    nothing
end


text:
using LinearAlgebra

function diff_svd(A,dA) #for a single parameter
    A_svd = svd(A)
    U = A_svd.U
    V = A_svd.V
    S = A_svd.S

    dS = U'*dA*V
    s = size(dS,1)

    Up = triu(dS,1) #extract the upper triangular part of the dS without the diagonal
    Lo = tril(dS,-1) #extract the lower triangular part of the dS matrix without diagonal

    L2_bar = zeros(s,s)

    for i = 2:s,j=1:(i-1)
        L2_bar[i,j] = (Up[j,i]*S[j] + Lo[i,j]*S[i])/(S[i]^2- S[j]^2)
    end

    dV = V*(L2_bar' - L2_bar)
    return diag(dS),dV
end

function jacob_svd(A,dA)
    p = size(dA,2)
    if isempty(size(A))
        n = 1
        m = 1
    else
        n,m = size(A)
    end

    dV = zeros(n*m,p)
    dS = zeros(max(n,m),p)

    for k in 1:p
        dA_unvec = dA[:,k]
        dA_unvec = reshape(dA_unvec,n,m)
        dS_temp,dV_temp = diff_svd(A,dA_unvec)
        dS[:,k] = dS_temp
        dV[:,k] = vec(dV_temp)
    end

    return dS,dV

end

mutable struct svd_kalman_diff
    # ALL V MATRICES FROM THE SVD ARE STORE TRANSPOSED I.E. Vt
    A # state transition
    G #Observable transition

    #shock of observable variance
    S_Q
    V_Q

    #variance of the shock of the state
    S_R
    V_R

    state_mean

    #variance of the state

    S_var
    V_var

    dA
    dG

    dS_Q
    dV_Q

    dS_R
    dV_R

    dS_var
    dV_var

    diff_state

    D_re #this will always be the square root of what we need
    e_bar

end

function svd_kalman_diff(A,G,Q,R,dA,dG,dR,dQ)

    m,n = size(G)
    st_mean = zeros(n)
    st_var = zeros(n,n)

    dS_Q, dV_Q = jacob_svd(Q,dQ)
    dS_R, dV_R = jacob_svd(R,dR)

    diff_state = zeros(n^2,size(dA,2))

    svd_Q = svd(Q)
    svd_R = svd(R)

    S_R = diagm(sqrt.(svd_R.S))
    S_Q = diagm(sqrt.(svd_Q.S)[1:rank(Q)])

    V_R = svd_R.Vt
    V_Q = svd_Q.Vt

    dS_R = inv(kron(I(m),S_R) + kron(S_R,I(m)))*dS_R
    dS_Q = inv(kron(I(rank(Q)),S_Q) + kron(S_Q,I(rank(Q))))*dS_Q[1:rank(Q),:]

    D_re = zeros(m,m)
    e_bar = zeros(m)

    return svd_kalman_diff(A,G,S_Q,V_Q,S_R,V_R,st_mean,st_var,st_var,dA,dG,dS_Q,dV_Q,dS_R,dV_R,dS_Q,dV_Q,D_re,e_r)
end

function start!(svd_kalman_diff::svd_kalman_diff,x0,v0,dv0)
    m,n = size(svd_kalman_diff.G)

    svd_kalman_diff.state_mean = x0

    svd_v0 = svd(v0)
    S_var = diagm(sqrt.(svd_v0.S))
    svd_kalman_diff.S_var = S_var
    svd_kalman_diff.V_var = svd_v0.Vt

    svd_kalman_diff.diff_state = zeros(size(x0,1))

    dS_var,dV_var = jacob_svd(v0,dv0)

    dS_var = inv(kron(I(n),S_var) + kron(S_var,I(n)))*dS_var

    svd_kalman_diff.dS_var = dS_var
    svd_kalman_diff.dV_var = dV_var

    return nothing
end

function forecast!(svd_kalman_diff::svd_kalman_diff)
    m,n = size(svd_kalman_diff.G)
    Q = svd_kalman_diff.S_Q*svd_kalman_diff.V_Q
    A = svd_kalman_diff.A
    old_vari = svd_kalman_diff.S_var*svd_kalman_diff.V_var
    old_state = svd_kalman_diff.state_mean

    dSvar = svd_kalman_diff.diff_Svar
    dQvar = svd_kalman_diff.diff_Vvar

    new_state = A*old_state

    pre_var = [old_vari*A';Q]

    svd_var = svd(pre_var)

    svd_kalman_diff.state_mean = new_state
    svd_kalman_diff.S_var = diagm(svd_var.S)
    svd_kalman_diff.V_var = svd_var.Vt

    b1 = kron(A*svd_kalman_diff.V_var,I(n))*svd_kalman_diff.dS_var + kron(A,svd_kalman_diff.S_var)*svd_kalman_diff.dV_var + kron(I(n),old_vari)*commutation_matrix(n,n)*svd_kalman_diff.dA
    b2 = kron(svd_kalman_diff.V_Q,I(n))*svd_kalman_diff.dS_Q + kron(I,svd_kalman_diff.S_Q)*svd_kalman_diff.dV_Q
    dVar = [b1;b2]

    svd_kalman_diff.Svar,svd_kalman_diff.Vvar = jacob_svd(diagm(svd_var.S)*svd_var.Vt,dVar)

    svd_kalman_diff.diff_state = kron(new_state',I(n))*svd_kalman_diff.dA + A*svd_kalman_diff.diff_state

    return nothing
end

function update_state!(svd_kalman_diff::svd_kalman_diff,y)
    R = svd_kalman_diff.S_R*svd_kalman_diff.V_R
    old_var = svd_kalman_diff.S_var*svd_kalman_diff.V_var
    G = svd_kalman_diff.G
    old_state = svd_kalman_diff.state_mean
    d_var = kron(svd_kalman_diff.V_var',I(n))*svd_kalman_diff.dS_var + kron(I(n),svd_kalman_diff.S_var)*svd_kalman_diff.dV_var
    m,n = size(G)

    re_mat = [R;old_var*G]
    svd_re = svd(re_mat)

    K_bar = (old_var'*old_var)*G'*svd_re.Vt
    er = *(y - G*old_state)
    e_bar = svd_re.Vt*er

    svd_kalman_diff.e_bar = e_bar
    svd_kalman_diff.D_re = diagm(svd_re.S)

    new_state = old_state + K_bar*diagm(svd_re.S.^(-2))*e_bar

    K = K_bar*diagm(svd_re.S.^(-2))*svd_re.Vt

    pre_var = [old_var*(I - K*G)';R*K']

    svd_var = svd(pre_var)

    svd_kalman_diff.S_var = diagm(svd_var.S)
    svd_kalman_diff.V_var = svd_var.Vt

    b1 = kron(svd_kalman_diff.V_R,I(n))*svd_kalman_diff.dS_R + kron(I(n),svd_kalman_diff.S_R)*commutation_matrix(n,n)*svd_kalman_diff.dV_R

    b2 = kron(G*svd_kalman_diff.V_var,I(n))*svd_kalman_diff.dS_var + kron(G,svd_kalman_diff.S_var)*commutation_matrix(n,n)*svd_kalman_diff.dV_var + kron(I,svd_kalman_diff.S*svd_kalman_diff.V)*commutation_matrix(m,n)*svd_kalman_diff.dG

    diff_re = [b1;b2]
    diff_svd_re = jacob_svd(re_mat,diff_re)

    diff_k_bar = kron(svd_re.Vt*G,I(n))*d_var + kron(svd_re.Vt,old_var)*svd_kalman_diff.dG + kron(I(n),old_var*G')*diff_svd_re.dV

    diff_k = kron(svd_re.V*diagm(svd_re.S .^(-2)),I(n))*diff_k_bar - kron(svd_re.V,K_bar)*(kron(diagm(svd_re.S .^(-2)),diagm(svd_re.S)) + kron(diagm(svd_re.S),diagm(svd_re.S .^(-2))))*diff_svd_re.dS + kron(I(n),K_bar*diagm(svd_re.S .^(-2)))*diff_svd_re.dV

    e_diff = -(kron(old_state',I(m))*svd_kalman_diff.dG + G*svd_kalman_diff.diff_state)
    e_bar_diff = kron(er',I(m))*diff_svd_re.dV + svd_re.Vt*e_diff

    svd_kalman_diff.diff_state = svd_kalman_diff.diff_state + kron(e_bar'*diagm(svd_re.S .^(-2)),I(n))*diff_k_bar - kron(e_bar',K_bar)*(kron(diagm(svd_re.S .^(-2)),diagm(svd_re.S)) + kron(diagm(svd_re.S),diagm(svd_re.S .^(-2))))*diff_svd_re.dS + K_bar*diagm(svd_re.S .^(-2))*e_bar_diff

    b0 = I(n) - K*G
    b1 = kron(b0*svd_kalman_diff.V_var,I(n))*svd_kalman_diff.dS_var + kron(b0,svd_kalman_svd.S_var)*svd_kalman_diff.dV_var - kron(K,svd_kalman_diff.S_var*svd_kalman_diff.V_var)*commutation_matrix(m,n)*svd_kalman_diff.dG - kron(I(n),svd_kalman_diff.S_var*svd_kalman_diff.V_var*G')*commutation_matrix(n,m)*dK
    b2 = kron(K*svd_kalman_diff.V_R',I(m))*svd_kalman_diff.dS_R + kron(K,svd_kalman_diff.S_R)*svd_kalman_res.dQ_R + kron(I(n),svd_kalman_diff.D_R*svd_kalman_diff.Q_R)*commutation_matrix(n,m)*dK

    diff_var = [b1;b2]

    var_pre = [old_var*b0';svd_kalman_diff.S_R*svd_kalman_diff.V_R*K']

    svd_var = svd(var_pre)

    svd_kalman_diff.S_var = diagm(svd_var.S)
    svd_kalman_diff.V_var = svd_var.Vt

    svd_kalman_diff.dS_var,svd_kalman_diff.dV_var = svd_jacob(var_pre,diff_var)

    return nothing

end

function update_and_forecast!(svd_kalman_diff::svd_kalman_diff,y)
    update_state!(svd_kalman_diff,y)
    forecast!(svd_kalman_diff)
    return nothing
end

max_stars_repo_name: danmrc/DSGE-HMC

max_issues_repo_name: danmrc/DSGE-HMC

max_forks_repo_name: danmrc/DSGE-HMC


text:
#=
using SimpleValueGraphs.AbstractTuples

import SimpleValueGraphs: default_eltype

function testset_toplogical_equivalent(g::SimpleGraph, gv::ValGraph)
    @testset "Topological equivalent" begin
        @test all(lr -> lr[1] == lr[2], zip(g.fadjlist, gv.fadjlist))
    end
end

function testset_toplogical_equivalent(g::SimpleDiGraph, gv::ValOutDiGraph)
    @testset "Topological equivalent" begin
        @test all(lr -> lr[1] == lr[2], zip(g.fadjlist, gv.fadjlist))
    end
end

function testset_toplogical_equivalent(g::SimpleDiGraph, gv::ValDiGraph)
    @testset "Topological equivalent" begin
        @test all(lr -> lr[1] == lr[2], zip(g.fadjlist, gv.fadjlist))
        @test all(lr -> lr[1] == lr[2], zip(g.badjlist, gv.badjlist))
    end
end

=#

@testset "constructor type stability" begin

    @inferred ValGraph(1)
    @inferred ValDiGraph(1)
    @inferred ValOutDiGraph(1)

    @inferred ValGraph{Int8}(1)
    @inferred ValDiGraph{UInt8}(1)
    @inferred ValOutDiGraph{Int16}(1)

    @inferred ValGraph{Int8, Tuple{Int64, String}, Tuple{}}(1, vertexval_init=undef)
    @inferred ValDiGraph{UInt8, Tuple{}, @NamedTuple{a::Union{Missing, Int}, b::Vector}}(1)
    @inferred ValOutDiGraph{Int16, @NamedTuple{a::Int, b::Int}, Tuple{String}}(1, vertexval_init=undef)
end

#  ------------------------------------------------------
#  Constructors from other value graphs
#  ------------------------------------------------------

@testset "ValGraph{$V_OUT}(::ValGraph{$V_IN}" for
            V_IN ∈ (Int64, Int32, Int8, UInt8), V_OUT ∈ (Int64, Int32, Int8, UInt8)

    g1 = ValGraph{V_IN}(5;
        vertexval_types=(Int8,),
        vertexval_init=v -> (v,),
        edgeval_types=(a=Int8, b=Int16),
    )
    add_edge!(g1, 1, 2, (a=Int8(1), b=Int16(2)))
    add_edge!(g1, 2, 3, (a=Int8(2), b=Int16(3)))
    add_edge!(g1, 5, 5, (a=Int8(51), b=Int16(52)))

    g2 = ValGraph{V_OUT}(g1)
    @test g2 isa ValGraph{V_OUT, vertexvals_type(g1), edgevals_type(g1), graphvals_type(g1)}
    @test ne(g2) == ne(g1)
    @test g2.fadjlist == g1.fadjlist
    @test g2.vertexvals == g1.vertexvals
    @test g2.edgevals == g1.edgevals
    @test g2.graphvals == g1.graphvals
end

@testset "ValGraph(::ValGraph{$V_IN}" for V_IN ∈ (Int64, Int32, Int8, UInt8)

    g1 = ValGraph{V_IN}(5;
        vertexval_types=(Int8,),
        vertexval_init=v -> (v,),
        edgeval_types=(Int8, Int16),
    )
    add_edge!(g1, 1, 2, (Int8(1), Int16(2)))
    add_edge!(g1, 2, 3, (Int8(2), Int16(3)))
    add_edge!(g1, 5, 5, (Int8(51), Int16(52)))

    g2 = ValGraph(g1)
    @test g2 isa ValGraph{eltype(g1), vertexvals_type(g1), edgevals_type(g1), graphvals_type(g1)}
    @test ne(g2) == ne(g1)
    @test g2.fadjlist == g1.fadjlist
    @test g2.vertexvals == g1.vertexvals
    @test g2.edgevals == g1.edgevals
    @test g2.graphvals == g1.graphvals
end

@testset "ValDiGraph{$V_OUT}(::ValGraph{$V_IN}" for
            V_IN ∈ (Int64, Int32, Int8, UInt8), V_OUT ∈ (Int64, Int32, Int8, UInt8)

    g1 = ValGraph{V_IN}(5;
        vertexval_types=(x=Int8,),
        vertexval_init=v -> (x=Int8(v),),
        edgeval_types=(a=Int8, b=Int16),
    )
    add_edge!(g1, 1, 2, (a=Int8(1), b=Int16(2)))
    add_edge!(g1, 2, 3, (a=Int8(2), b=Int16(3)))
    add_edge!(g1, 5, 5, (a=Int8(51), b=Int16(52)))

    g2 = ValDiGraph{V_OUT}(g1)
    @test g2 isa ValDiGraph{V_OUT, vertexvals_type(g1), edgevals_type(g1), graphvals_type(g1)}
    @test ne(g2) == 5
    @test g2.fadjlist == g2.badjlist == [[2], [1, 3], [2], Int[], [5]]
    @test g2.vertexvals == g1.vertexvals
    @test g2.edgevals == g2.redgevals == g1.edgevals
    @test g2.graphvals == g1.graphvals
end

@testset "ValDiGraph(::ValGraph{$V_IN}" for V_IN ∈ (Int64, Int32, Int8, UInt8)

    g1 = ValGraph{V_IN}(5;
        vertexval_types=(x=Int8,),
        vertexval_init=v -> (x=Int8(v),),
        edgeval_types=(a=Int8, b=Int16),
    )
    add_edge!(g1, 1, 2, (a=Int8(1), b=Int16(2)))
    add_edge!(g1, 2, 3, (a=Int8(2), b=Int16(3)))
    add_edge!(g1, 5, 5, (a=Int8(51), b=Int16(52)))

    g2 = ValDiGraph(g1)
    @test g2 isa ValDiGraph{eltype(g1), vertexvals_type(g1), edgevals_type(g1), graphvals_type(g1)}
    @test ne(g2) == 5
    @test g2.fadjlist == g2.badjlist == [[2], [1, 3], [2], Int[], [5]]
    @test g2.vertexvals == g1.vertexvals
    @test g2.edgevals == g2.redgevals == g1.edgevals
    @test g2.graphvals == g1.graphvals
end

@testset "ValOutDiGraph{$V_OUT}(::ValGraph{$V_IN}" for
            V_IN ∈ (Int64, Int32, Int8, UInt8), V_OUT ∈ (Int64, Int32, Int8, UInt8)

    g1 = ValGraph{V_IN}(5;
        vertexval_types=(x=Int8,),
        vertexval_init=v -> (x=Int8(v),),
        edgeval_types=(a=Int8, b=Int16),
    )
    add_edge!(g1, 1, 2, (a=Int8(1), b=Int16(2)))
    add_edge!(g1, 2, 3, (a=Int8(2), b=Int16(3)))
    add_edge!(g1, 5, 5, (a=Int8(51), b=Int16(52)))

    g2 = ValOutDiGraph{V_OUT}(g1)
    @test g2 isa ValOutDiGraph{V_OUT, vertexvals_type(g1), edgevals_type(g1), graphvals_type(g1)}
    @test ne(g2) == 5
    @test g2.fadjlist == [[2], [1, 3], [2], Int[], [5]]
    @test g2.vertexvals == g1.vertexvals
    @test g2.edgevals == g1.edgevals
    @test g2.graphvals == g1.graphvals
end

@testset "ValOutDiGraph(::ValGraph{$V_IN}" for V_IN ∈ (Int64, Int32, Int8, UInt8)

    g1 = ValGraph{V_IN}(5;
        vertexval_types=(x=Int8,),
        vertexval_init=v -> (x=Int8(v),),
        edgeval_types=(a=Int8, b=Int16),
    )
    add_edge!(g1, 1, 2, (a=Int8(1), b=Int16(2)))
    add_edge!(g1, 2, 3, (a=Int8(2), b=Int16(3)))
    add_edge!(g1, 5, 5, (a=Int8(51), b=Int16(52)))

    g2 = ValOutDiGraph(g1)
    @test g2 isa ValOutDiGraph{eltype(g1), vertexvals_type(g1), edgevals_type(g1), graphvals_type(g1)}
    @test ne(g2) == 5
    @test g2.fadjlist == [[2], [1, 3], [2], Int[], [5]]
    @test g2.vertexvals == g1.vertexvals
    @test g2.edgevals == g1.edgevals
    @test g2.graphvals == g1.graphvals
end

@testset "ValDiGraph{$V_OUT}(::ValDiGraph{$V_IN}" for
            V_IN ∈ (Int64, Int32, Int8, UInt8), V_OUT ∈ (Int64, Int32, Int8, UInt8)

    g1 = ValDiGraph{V_IN}(5;
        vertexval_types=(x=Int, y=Int),
        vertexval_init=v -> (x=Int(v), y=Int(10 * v)),
        edgeval_types=(a=Int8, b=Int16),
    )

    add_edge!(g1, 1, 2, (a=Int8(1), b=Int16(2)))
    add_edge!(g1, 2, 3, (a=Int8(2), b=Int16(3)))
    add_edge!(g1, 3, 2, (a=Int8(3), b=Int16(2)))
    add_edge!(g1, 5, 5, (a=Int8(51), b=Int16(52)))

    g2 = ValDiGraph{V_OUT}(g1)
    @test g2 isa ValDiGraph{V_OUT, vertexvals_type(g1), edgevals_type(g1), graphvals_type(g1)}
    @test ne(g2) == ne(g1)
    @test g2.fadjlist == g1.fadjlist
    @test g2.badjlist == g1.badjlist
    @test g2.vertexvals == g1.vertexvals
    @test g2.edgevals == g1.edgevals
    @test g2.redgevals == g1.redgevals
    @test g2.graphvals == g1.graphvals
end

@testset "ValDiGraph(::ValDiGraph{$V_IN}" for V_IN ∈ (Int64, Int32, Int8, UInt8)

    g1 = ValDiGraph{V_IN}(5;
        vertexval_types=(x=Int, y=Int),
        vertexval_init=v -> (x=Int(v), y=Int(10 * v)),
        edgeval_types=(a=Int8, b=Int16),
    )

    add_edge!(g1, 1, 2, (a=Int8(1), b=Int16(2)))
    add_edge!(g1, 2, 3, (a=Int8(2), b=Int16(3)))
    add_edge!(g1, 3, 2, (a=Int8(3), b=Int16(2)))
    add_edge!(g1, 5, 5, (a=Int8(51), b=Int16(52)))

    g2 = ValDiGraph(g1)
    @test g2 isa ValDiGraph{eltype(g1), vertexvals_type(g1), edgevals_type(g1), graphvals_type(g1)}
    @test ne(g2) == ne(g1)
    @test g2.fadjlist == g1.fadjlist
    @test g2.badjlist == g1.badjlist
    @test g2.vertexvals == g1.vertexvals
    @test g2.edgevals == g1.edgevals
    @test g2.redgevals == g1.redgevals
    @test g2.graphvals == g1.graphvals
end

@testset "ValOutDiGraph{$V_OUT}(::ValOutDiGraph{$V_IN}" for
            V_IN ∈ (Int64, Int32, Int8, UInt8), V_OUT ∈ (Int64, Int32, Int8, UInt8)

    g1 = ValOutDiGraph{V_IN}(5;
        vertexval_types=(x=Int, y=Int),
        vertexval_init=v -> (x=Int(v), y=Int(10 * v)),
        edgeval_types=(a=Int8, b=Int16),
    )

    add_edge!(g1, 1, 2, (a=Int8(1), b=Int16(2)))
    add_edge!(g1, 2, 3, (a=Int8(2), b=Int16(3)))
    add_edge!(g1, 3, 2, (a=Int8(3), b=Int16(2)))
    add_edge!(g1, 5, 5, (a=Int8(51), b=Int16(52)))

    g2 = ValOutDiGraph{V_OUT}(g1)
    @test g2 isa ValOutDiGraph{V_OUT, vertexvals_type(g1), edgevals_type(g1), graphvals_type(g1)}
    @test ne(g2) == ne(g1)
    @test g2.fadjlist == g1.fadjlist
    @test g2.vertexvals == g1.vertexvals
    @test g2.edgevals == g1.edgevals
    @test g2.graphvals == g1.graphvals
end

@testset "ValOutDiGraph(::ValOutDiGraph{$V_IN}" for V_IN ∈ (Int64, Int32, Int8, UInt8)

    g1 = ValDiGraph{V_IN}(5;
        vertexval_types=(x=Int, y=Int),
        vertexval_init=v -> (x=Int(v), y=Int(10 * v)),
        edgeval_types=(a=Int8, b=Int16),
    )

    add_edge!(g1, 1, 2, (a=Int8(1), b=Int16(2)))
    add_edge!(g1, 2, 3, (a=Int8(2), b=Int16(3)))
    add_edge!(g1, 3, 2, (a=Int8(3), b=Int16(2)))
    add_edge!(g1, 5, 5, (a=Int8(51), b=Int16(52)))

    g2 = ValOutDiGraph(g1)
    @test g2 isa ValOutDiGraph{eltype(g1), vertexvals_type(g1), edgevals_type(g1), graphvals_type(g1)}
    @test ne(g2) == ne(g1)
    @test g2.fadjlist == g1.fadjlist
    @test g2.vertexvals == g1.vertexvals
    @test g2.edgevals == g1.edgevals
    @test g2.graphvals == g1.graphvals
end

@testset "ValOutDiGraph{$V_OUT}(::ValDiGraph{$V_IN}" for
            V_IN ∈ (Int64, Int32, Int8, UInt8), V_OUT ∈ (Int64, Int32, Int8, UInt8)

    g1 = ValOutDiGraph{V_IN}(5;
        vertexval_types=(x=Int, y=Int),
        vertexval_init=v -> (x=Int(v), y=Int(10 * v)),
        edgeval_types=(a=Int8, b=Int16),
    )

    add_edge!(g1, 1, 2, (a=Int8(1), b=Int16(2)))
    add_edge!(g1, 2, 3, (a=Int8(2), b=Int16(3)))
    add_edge!(g1, 3, 2, (a=Int8(3), b=Int16(2)))
    add_edge!(g1, 5, 5, (a=Int8(51), b=Int16(52)))

    g2 = ValOutDiGraph{V_OUT}(g1)
    @test g2 isa ValOutDiGraph{V_OUT, vertexvals_type(g1), edgevals_type(g1), graphvals_type(g1)}
    @test ne(g2) == ne(g1)
    @test g2.fadjlist == g1.fadjlist
    @test g2.vertexvals == g1.vertexvals
    @test g2.edgevals == g1.edgevals
    @test g2.graphvals == g1.graphvals
end

@testset "ValOutDiGraph(::ValDiGraph{$V_IN}" for V_IN ∈ (Int64, Int32, Int8, UInt8)

    g1 = ValDiGraph{V_IN}(5;
        vertexval_types=(x=Int, y=Int),
        vertexval_init=v -> (x=Int(v), y=Int(10 * v)),
        edgeval_types=(a=Int8, b=Int16),
    )

    add_edge!(g1, 1, 2, (a=Int8(1), b=Int16(2)))
    add_edge!(g1, 2, 3, (a=Int8(2), b=Int16(3)))
    add_edge!(g1, 3, 2, (a=Int8(3), b=Int16(2)))
    add_edge!(g1, 5, 5, (a=Int8(51), b=Int16(52)))

    g2 = ValOutDiGraph(g1)
    @test g2 isa ValOutDiGraph{eltype(g1), vertexvals_type(g1), edgevals_type(g1), graphvals_type(g1)}
    @test ne(g2) == ne(g1)
    @test g2.fadjlist == g1.fadjlist
    @test g2.vertexvals == g1.vertexvals
    @test g2.edgevals == g1.edgevals
    @test g2.graphvals == g1.graphvals
end


text:
function Dimer_Eq_of_M(du, u, p, t)
    du[1] = (-4) * (2 + u[3]^2 + (-2) * 2^(1 / 2) * u[4] + u[4]^2 + 2 * u[3] * ((-1) * 2^(1 / 2) + u[4])   + (u[1] + (-1) * u[2])^2)^(-2) * (2 + u[3]^2 + 2 * 2^(1 / 2) * u[4] + u[4]^2 + 2 * u[3] * (   2^(1 / 2) + u[4]) + (u[1] + (-1) * u[2])^2)^(-2) * (u[4]^2 + u[2]^2) * (2 + u[3]^2 + 2 *   u[3] * (2^(1 / 2) + (-1) * u[4]) + (-2) * 2^(1 / 2) * u[4] + u[4]^2 + (u[1] + u[2])^2)^(-2)     *(2 + u[3]^2 + 2 * 2^(1 / 2) * u[4] + u[4]^2 + (-2) * u[3] * (2^(1 / 2) + u[4]) + (u[1] + u[2])^2)^(-2) * ((u[3]^2 + u[1]^2) * (2 + (-2) * 2^(1 / 2) * (u[3] + u[4]) + (u[3] + u[4])^2 +   (u[1] + (-1) * u[2])^2)^(-1) * (2 + 2 * 2^(1 / 2) * (u[3] + u[4]) + (u[3] + u[4])^2 + (u[1] + (  -1) * u[2])^2)^(-1) * (u[4]^2 + u[2]^2) * (2 + 2 * 2^(1 / 2) * (u[3] + (-1) * u[4]) + (   u[3] + (-1) * u[4])^2 + (u[1] + u[2])^2)^(-1) * (2 + (u[3] + (-1) * u[4])^2 + 2 * 2^(1 / 2)   * ((-1) * u[3] + u[4]) + (u[1] + u[2])^2)^(-1))^(-1 / 2) * (3 * u[3]^9 + 24 * u[3]^6 *   u[4] * u[1] * u[2] + 4 * u[3]^7 * ((-4) + (-2) * u[4]^2 + 3 * u[1]^2 + 2 * u[2]^2) + 8 *   u[3]^4 * u[4] * u[1] * u[2] * (4 + (-2) * u[4]^2 + 7 * u[1]^2 + 2 * u[2]^2) + 4 * u[3]^3 *   u[1]^2 * (12 + 3 * u[4]^4 + 3 * u[1]^4 + (-4) * u[2]^2 + 3 * u[2]^4 + 2 * u[1]^2 * ((-2) +   u[2]^2) + (-2) * u[4]^2 * ((-2) + u[1]^2 + u[2]^2)) + 8 * u[4] * u[1]^3 * u[2] * (u[4]^4 + (   u[1]^2 + (-1) * u[2]^2)^2 + 2 * u[4]^2 * (2 + u[1]^2 + u[2]^2) + (-4) * (3 + u[1]^2 +   u[2]^2)) + 2 * u[3]^5 * (12 + 3 * u[4]^4 + 9 * u[1]^4 + (-4) * u[2]^2 + 3 * u[2]^4 + 8 *   u[1]^2 * ((-2) + u[2]^2) + (-2) * u[4]^2 * ((-2) + 4 * u[1]^2 + u[2]^2)) + (-8) *   u[3]^2 * u[4] * u[1] * u[2] * (u[4]^4 + (-5) * u[1]^4 + u[2]^4 + (-8) * u[1]^2 * (2 + u[2]^2)   + (-4) * (3 + u[2]^2) + 2 * u[4]^2 * (2 + 4 * u[1]^2 + u[2]^2)) + u[3] * ((-1) * u[4]^8 +  3 * u[1]^8 + 8 * u[1]^2 * ((-2) + u[2]^2) * (2 + u[2]^2)^2 + (-1) * (2 + u[2]^2)^4 + (  -4) * u[4]^6 * ((-2) + 2 * u[1]^2 + u[2]^2) + (-2) * u[1]^4 * (20 + 36 * u[2]^2 + 5 *   u[2]^4) + (-2) * u[4]^4 * (12 + 5 * u[1]^4 + (-4) * u[2]^2 + 3 * u[2]^4 + 4 * u[1]^2 * ((  -2) + u[2]^2)) + 4 * u[4]^2 * ((-1) * ((-2) + u[2]^2) * (2 + u[2]^2)^2 + u[1]^4 * (  18 + 23 * u[2]^2) + 2 * u[1]^2 * (4 + 20 * u[2]^2 + u[2]^4))));
    du[2] = 4 * (u[3]^2 + u[1]^2) * (2 + u[3]^2 + (-2) * 2^(1 / 2) * u[4] + u[4]^2 + 2 * u[3] * ((-1) *   2^(1 / 2) + u[4]) + (u[1] + (-1) * u[2])^2)^(-2) * (2 + u[3]^2 + 2 * 2^(1 / 2) * u[4] +   u[4]^2 + 2 * u[3] * (2^(1 / 2) + u[4]) + (u[1] + (-1) * u[2])^2)^(-2) * (2 + u[3]^2 + 2 *   u[3] * (2^(1 / 2) + (-1) * u[4]) + (-2) * 2^(1 / 2) * u[4] + u[4]^2 + (u[1] + u[2])^2)^(-2)     *(2 + u[3]^2 + 2 * 2^(1 / 2) * u[4] + u[4]^2 + (-2) * u[3] * (2^(1 / 2) + u[4]) + (u[1] + u[2])^2)^(-2) * ((u[3]^2 + u[1]^2) * (2 + (-2) * 2^(1 / 2) * (u[3] + u[4]) + (u[3] + u[4])^2 +   (u[1] + (-1) * u[2])^2)^(-1) * (2 + 2 * 2^(1 / 2) * (u[3] + u[4]) + (u[3] + u[4])^2 + (u[1] + (   -1) * u[2])^2)^(-1) * (u[4]^2 + u[2]^2) * (2 + 2 * 2^(1 / 2) * (u[3] + (-1) * u[4]) + (   u[3] + (-1) * u[4])^2 + (u[1] + u[2])^2)^(-1) * (2 + (u[3] + (-1) * u[4])^2 + 2 * 2^(1 / 2)   * ((-1) * u[3] + u[4]) + (u[1] + u[2])^2)^(-1))^(-1 / 2) * ((-3) * u[4]^9 + (-24) *   u[3] * u[4]^6 * u[1] * u[2] + 8 * u[3] * u[4]^4 * u[1] * u[2] * ((-4) + 2 * u[3]^2 + (-2) *   u[1]^2 + (-7) * u[2]^2) + 4 * u[4]^7 * (4 + 2 * u[3]^2 + (-2) * u[1]^2 + (-3) * u[2]^2) +   u[4] * ((u[3]^4 + 2 * u[3]^2 * ((-2) + u[1]^2) + (2 + u[1]^2)^2)^2 + 8 * (u[3]^6 +   u[3]^4 * ((-2) + u[1]^2) + (-1) * ((-2) + u[1]^2) * (2 + u[1]^2)^2 + (-1) * u[3]^2 *   (4 + 20 * u[1]^2 + u[1]^4)) * u[2]^2 + 2 * (20 + 5 * u[3]^4 + 36 * u[1]^2 + 5 * u[1]^4 + (   -2) * u[3]^2 * (18 + 23 * u[1]^2)) * u[2]^4 + (-3) * u[2]^8) + (-4) * u[4]^3 *   u[2]^2 * (12 + 3 * u[3]^4 + (-4) * u[1]^2 + 3 * u[1]^4 + 2 * ((-2) + u[1]^2) * u[2]^2 +   3 * u[2]^4 + (-2) * u[3]^2 * ((-2) + u[1]^2 + u[2]^2)) + (-8) * u[3] * u[1] * u[2]^3 * (   u[3]^4 + (u[1]^2 + (-1) * u[2]^2)^2 + 2 * u[3]^2 * (2 + u[1]^2 + u[2]^2) + (-4) * (3 +   u[1]^2 + u[2]^2)) + (-2) * u[4]^5 * (12 + 3 * u[3]^4 + (-4) * u[1]^2 + 3 * u[1]^4 + 8 * ((   -2) + u[1]^2) * u[2]^2 + 9 * u[2]^4 + (-2) * u[3]^2 * ((-2) + u[1]^2 + 4 * u[2]^2)) + 8 *   u[3] * u[4]^2 * u[1] * u[2] * ((-12) + u[3]^4 + (-4) * u[1]^2 + u[1]^4 + (-8) * (2 + u[1]^2)   * u[2]^2 + (-5) * u[2]^4 + 2 * u[3]^2 * (2 + u[1]^2 + 4 * u[2]^2)));
    du[4] = 4 * (2 + u[3]^2 + (-2) * 2^(1 / 2) * u[4] + u[4]^2 + 2 * u[3] * ((-1) * 2^(1 / 2) + u[4]) + ( u[1] + (-1) * u[2])^2)^(-2) * (2 + u[3]^2 + 2 * 2^(1 / 2) * u[4] + u[4]^2 + 2 * u[3] * ( 2^(1 / 2) + u[4]) + (u[1] + (-1) * u[2])^2)^(-2) * (u[4]^2 + u[2]^2) * (2 + u[3]^2 + 2 * u[3] * (2^(1 / 2) + (-1) * u[4]) + (-2) * 2^(1 / 2) * u[4] + u[4]^2 + (u[1] + u[2])^2)^(-2)     *(2 + u[3]^2 + 2 * 2^(1 / 2) * u[4] + u[4]^2 + (-2) * u[3] * (2^(1 / 2) + u[4]) + (u[1] + u[2])^2)^(-2) * ((u[3]^2 + u[1]^2) * (2 + (-2) * 2^(1 / 2) * (u[3] + u[4]) + (u[3] + u[4])^2 + (u[1] + (-1) * u[2])^2)^(-1) * (2 + 2 * 2^(1 / 2) * (u[3] + u[4]) + (u[3] + u[4])^2 + (u[1] + ( -1) * u[2])^2)^(-1) * (u[4]^2 + u[2]^2) * (2 + 2 * 2^(1 / 2) * (u[3] + (-1) * u[4]) + ( u[3] + (-1) * u[4])^2 + (u[1] + u[2])^2)^(-1) * (2 + (u[3] + (-1) * u[4])^2 + 2 * 2^(1 / 2) * ((-1) * u[3] + u[4]) + (u[1] + u[2])^2)^(-1))^(-1 / 2) * (3 * u[1]^9 + 24 * u[3] * u[4] *u[1]^6 * u[2] + 4 * u[1]^7 * (4 + 3 * u[3]^2 + 2 * u[4]^2 + (-2) * u[2]^2) + 8 * u[3] * u[4] * u[1]^4 * u[2] * ((-4) + 7 * u[3]^2 + 2 * u[4]^2 + (-2) * u[2]^2) + 4 * u[3]^2 * u[1]^3 * ( 12 + 3 * u[3]^4 + 4 * u[4]^2 + 3 * u[4]^4 + (-2) * (2 + u[4]^2) * u[2]^2 + 3 * u[2]^4 + 2 * u[3]^2 * (2 + u[4]^2 + (-1) * u[2]^2)) + 2 * u[1]^5 * (12 + 9 * u[3]^4 + 4 * u[4]^2 + 3 * u[4]^4 + (-2) * (2 + u[4]^2) * u[2]^2 + 3 * u[2]^4 + 8 * u[3]^2 * (2 + u[4]^2 + (-1) * u[2]^2)) + 8 * u[3] * u[4] * u[1]^2 * u[2] * (5 * u[3]^4 + 8 * u[3]^2 * ((-2) + u[4]^2 + (-1) * u[2]^2) + 4 * (3 + (-1) * u[4]^2 + u[2]^2) + (-1) * (u[4]^2 + u[2]^2)^2) + 8 * u[3]^3 * u[4] * u[2] * (u[3]^4 + 4 * ((-3) + u[4]^2 + (-1) * u[2]^2) + 2 * u[3]^2 * (2 + ( -1) * u[4]^2 + u[2]^2) + (u[4]^2 + u[2]^2)^2) + u[1] * (3 * u[3]^8 + (-1) * (u[4]^4 + 2 * u[4]^2 * ((-2) + u[2]^2) + (2 + u[2]^2)^2)^2 + (-2) * u[3]^4 * (20 + 5 * u[4]^4 + 36 * u[2]^2 + 5 * u[2]^4 + (-2) * u[4]^2 * (18 + 23 * u[2]^2)) + 8 * u[3]^2 * (u[4]^6 + u[4]^4 * ((-2) + u[2]^2) + (-1) * ((-2) + u[2]^2) * (2 + u[2]^2)^2 + (-1) * u[4]^2 * (4 + 20 * u[2]^2 + u[2]^4))));
    du[3] = (-4) * (u[3]^2 + u[1]^2) * (2 + u[3]^2 + (-2) * 2^(1 / 2) * u[4] + u[4]^2 + 2 * u[3] * (( -1) * 2^(1 / 2) + u[4]) + (u[1] + (-1) * u[2])^2)^(-2) * (2 + u[3]^2 + 2 * 2^(1 / 2) * u[4] + u[4]^2 + 2 * u[3] * (2^(1 / 2) + u[4]) + (u[1] + (-1) * u[2])^2)^(-2) * ((-8) * u[3] * u[4]^3 * u[1] * (u[3]^4 + 4 * ((-3) + u[4]^2 + (-1) * u[1]^2) + 2 * u[3]^2 * (2 + (-1) * u[4]^2 + u[1]^2) + (u[4]^2 + u[1]^2)^2) + (u[3]^8 + (-3) * u[4]^8 + 8 * u[4]^2 * ((-2) + u[1]^2) * (2 + u[1]^2)^2 + (2 + u[1]^2)^4 + 4 * u[3]^6 * ((-2) + (-2) * u[4]^2 + u[1]^2) + 2 * u[3]^4 * (12 + 8 * u[4]^2 + 5 * u[4]^4 + (-4) * (1 + u[4]^2) * u[1]^2 + 3 * u[1]^4) + 2 * u[4]^4 * (20 + 36 * u[1]^2 + 5 * u[1]^4) + 4 * u[3]^2 * (((-2) + u[1]^2) * (2 + u[1]^2)^2 + (-1) * u[4]^4 * (18 + 23 * u[1]^2) + 2 * u[4]^2 * (4 + 20 * u[1]^2 + u[1]^4))) * u[2] + 8 * u[3] * u[4] * u[1] * (u[3]^4 + (-5) * u[4]^4 + u[1]^4 + 8 * u[4]^2 * (2 + u[1]^2) + (-4) * (3 + u[1]^2) + 2 * u[3]^2 * (2 + (-4) * u[4]^2 + u[1]^2)) * u[2]^2 + ( -4) * u[4]^2 * (12 + 3 * u[3]^4 + 4 * u[4]^2 + 3 * u[4]^4 + (-2) * (2 + u[4]^2) * u[1]^2 + 3 * u[1]^4 + 2 * u[3]^2 * (2 + u[4]^2 + (-1) * u[1]^2)) * u[2]^3 + (-8) * u[3] * u[4] * u[1] * (2 * u[3]^2 + 7 * u[4]^2 + (-2) * (2 + u[1]^2)) * u[2]^4 + (-2) * (12 + 3 * u[3]^4 + 9 * u[4]^4 + (-4) * u[1]^2 + 3 * u[1]^4 + u[3]^2 * (4 + 8 * u[4]^2 + (-2) * u[1]^2) + (-8) * u[4]^2 * ((-2) + u[1]^2)) * u[2]^5 + (-24) * u[3] * u[4] * u[1] * u[2]^6 + (-4) * (4 + 2 * u[3]^2 + 3 * u[4]^2 + (-2) * u[1]^2) * u[2]^7 + (-3) * u[2]^9) * (2 + u[3]^2 + 2 * u[3] * (2^(1 / 2) + (-1) * u[4]) + (-2) * 2^(1 / 2) * u[4] + u[4]^2 + (u[1] + u[2])^2)^(-2)     *(2 + u[3]^2 + 2 * 2^(1 / 2) * u[4] + u[4]^2 + (-2) * u[3] * (2^(1 / 2) + u[4]) + (u[1] + u[2])^2)^(-2) * ((u[3]^2 + u[1]^2) * (2 + (-2) * 2^(1 / 2) * (u[3] + u[4]) + (u[3] + u[4])^2 + (u[1] + (-1) * u[2])^2)^(-1) * (2 + 2 * 2^(1 / 2) * (u[3] + u[4]) + (u[3] + u[4])^2 + (u[1] + ( -1) * u[2])^2)^(-1) * (u[4]^2 + u[2]^2) * (2 + 2 * 2^(1 / 2) * (u[3] + (-1) * u[4]) + ( u[3] + (-1) * u[4])^2 + (u[1] + u[2])^2)^(-1) * (2 + (u[3] + (-1) * u[4])^2 + 2 * 2^(1 / 2) * ((-1) * u[3] + u[4]) + (u[1] + u[2])^2)^(-1))^(-1 / 2);
    du[5] = 0
    du[6] = 0
end

# Hamiltonian_Dimer(q,p,t)=4sqrt(abs( (q[1]^2+p[1]^2)*(q[2]^2+p[2]^2) )/
# (
#         (2-2sqrt(2)*(p[1]+p[2])+(p[1]+p[2])^2+(q[1]-q[2])^2 )*
#         (2+2sqrt(2)*(p[1]+p[2])+(p[1]+p[2])^2+(q[1]-q[2])^2 )*
#         (2-2sqrt(2)*(p[1]-p[2])+(p[1]-p[2])^2+(q[1]+q[2])^2 )*
#         (2+2sqrt(2)*(p[1]-p[2])+(p[1]-p[2])^2+(q[1]+q[2])^2 )

# ));
# function Dimer_Eq_of_M(du, u, p, t)
#     Q=[u[1],u[2]]; P=[u[3],u[4]];
#     Ham=Hamiltonian_Dimer(Q,P,0);
#     du[1]=Ham*(u[3]*(u[3]^2+u[1]^2)^(-1)+(-2)*(u[3]+u[3])*(u[3]^4+4*u[3]^3*u[3]+u[3]^4+ 2*u[3]^2*((-2)+(u[1]+(-1)*u[2])^2)+(2+(u[1]+(-1)*u[2])^2)^2+4*u[3]* u[3]*((-2)+u[3]^2+(u[1]+(-1)*u[2])^2)+2*u[3]^2*((-2)+3*u[3]^2+(u[1]+(-1)*u[2])^2))^(-1)*((-2)+(u[3]+u[3])^2+(u[1]+(-1)*u[2])^2)+(-2)*(u[3]+ (-1)*u[3])*((-2)+(u[3]+(-1)*u[3])^2+(u[1]+u[2])^2)*(u[3]^4+(-4)* u[3]^3*u[3]+u[3]^4+2*u[3]^2*((-2)+(u[1]+u[2])^2)+(2+(u[1]+u[2])^2)^2+(-4) *u[3]*u[3]*((-2)+u[3]^2+(u[1]+u[2])^2)+2*u[3]^2*((-2)+3*u[3]^2+(u[1]+u[2]) ^2))^(-1));
#     du[2]=Ham*((2+u[3]^2+(-2)*2^(1/2)*u[3]+u[3]^2+2*u[3]*((-1)*2^(1/2)+u[3])+(u[1]+(-1)*u[2])^2)^(-1)*(2+u[3]^2+2*2^(1/2)*u[3]+u[3]^2+2*u[3]*(2^(1/2)+u[3])+(u[1]+(-1)*u[2])^2)^(-1)*(u[3]^2+u[2]^2)^(-1)*(2+u[3]^2+  2*u[3]*(2^(1/2)+(-1)*u[3])+(-2)*2^(1/2)*u[3]+u[3]^2+(u[1]+u[2])^2)^(  -1)*(2+u[3]^2+2*2^(1/2)*u[3]+u[3]^2+(-2)*u[3]*(2^(1/2)+u[3])+(u[1]+  u[2])^2)^(-1)*((-3)*u[3]^9+(-24)*u[3]*u[3]^6*u[1]*u[2]+8*u[3]*u[3]^4*u[1]*u[2]*((-4)+2*u[3]^2+(-2)*u[1]^2+(-7)*u[2]^2)+4*u[3]^7*(4+2*u[3]^2+(-2)*u[1]^2+(-3)*u[2]^2)+u[3]*((u[3]^4+2*u[3]^2*((-2)+   u[1]^2)+(2+u[1]^2)^2)^2+8*(u[3]^6+u[3]^4*((-2)+u[1]^2)+(-1)*((-2)+ u[1]^2)*(2+u[1]^2)^2+(-1)*u[3]^2*(4+20*u[1]^2+u[1]^4))*u[2]^2+2*( 20+5*u[3]^4+36*u[1]^2+5*u[1]^4+(-2)*u[3]^2*(18+23*u[1]^2))*  u[2]^4+(-3)*u[2]^8)+(-4)*u[3]^3*u[2]^2*(12+3*u[3]^4+(-4)*u[1]^2+ 3*u[1]^4+2*((-2)+u[1]^2)*u[2]^2+3*u[2]^4+(-2)*u[3]^2*((-2)+u[1]^2+  u[2]^2))+(-8)*u[3]*u[1]*u[2]^3*(u[3]^4+(u[1]^2+(-1)*u[2]^2)^2+2* u[3]^2*(2+u[1]^2+u[2]^2)+(-4)*(3+u[1]^2+u[2]^2))+(-2)*u[3]^5*(12+3* u[3]^4+(-4)*u[1]^2+3*u[1]^4+8*((-2)+u[1]^2)*u[2]^2+9*u[2]^4+(-2)* u[3]^2*((-2)+u[1]^2+4*u[2]^2))+8*u[3]*u[3]^2*u[1]*u[2]*((-12)+u[3]^4+(-4)*u[1]^2+u[1]^4+(-8)*(2+u[1]^2)*u[2]^2+(-5)*u[2]^4+2*u[3]^2*(2+  u[1]^2+4*u[2]^2))));
#     du[3]=Ham*((-1)*u[1]*(u[3]^2+u[1]^2)^(-1)+2*(u[3]^4+4*u[3]^3*u[3]+u[3]^4+2* u[3]^2*((-2)+(u[1]+(-1)*u[2])^2)+(2+(u[1]+(-1)*u[2])^2)^2+4*u[3]*u[3]*  ((-2)+u[3]^2+(u[1]+(-1)*u[2])^2)+2*u[3]^2*((-2)+3*u[3]^2+(u[1]+(-1)*  u[2])^2))^(-1)*(2+(u[3]+u[3])^2+(u[1]+(-1)*u[2])^2)*(u[1]+(-1)*u[2])+2*(u[1]+u[2])*(2+(u[3]+(-1)*u[3])^2+(u[1]+u[2])^2)*(u[3]^4+(-4)*u[3]^3*u[3]+u[3]^4+2*u[3]^2*((-2)+(u[1]+u[2])^2)+(2+(u[1]+u[2])^2)^2+(-4)*u[3]*u[3]*((-2)+u[3]^2+(u[1]+u[2])^2)+2*u[3]^2*((-2)+3*u[3]^2+(u[1]+u[2])^2))^( -1));
#     du[4]=Ham*((2+u[3]^2+(-2)*2^(1/2)*u[3]+u[3]^2+2*u[3]*((-1)*2^(1/2)+u[3])+(u[1]+(-1)*u[2])^2)^(-1)*((-1)*u[1]+u[2])+(2+u[3]^2+2*2^(1/2)*u[3]+u[3]^2+  2*u[3]*(2^(1/2)+u[3])+(u[1]+(-1)*u[2])^2)^(-1)*((-1)*u[1]+u[2])+(-1)*  u[2]*(u[3]^2+u[2]^2)^(-1)+(u[1]+u[2])*(2+u[3]^2+2*u[3]*(2^(1/2)+(-1)* u[3])+(-2)*2^(1/2)*u[3]+u[3]^2+(u[1]+u[2])^2)^(-1)+(u[1]+u[2])*(2+u[3]^2+ 2*2^(1/2)*u[3]+u[3]^2+(-2)*u[3]*(2^(1/2)+u[3])+(u[1]+u[2])^2)^(-1));
#     du[5]=0
#     du[6]=0
# end


text:
export opengene_homedir,
    opengene_datadir,
    parse_config

include("folders.jl")


text:
module steadystate_classical

export OpticalPotential, evaluate_opticalpotential

using ..classical

type OpticalPotential
    N::Int
    k::Vector{Float64}
    ξ::Vector{Float64}
    function OpticalPotential(k::Vector{Float64}, ξ::Vector{Float64})
        @assert length(k)==length(ξ)
        x = new()
        x.N = length(k)
        x.k = k
        x.ξ = ξ
        return x
    end
end

OpticalPotential(s::MultimodeSystem, state::ClassicalState) = OpticalPotential(Float64[mode.k for mode=s.modes], Float64[2*s.modes[n].eta*αn(state,n)[1] for n=1:length(s.modes)])

function evaluate_opticalpotential(U::OpticalPotential, derivative::Int, x)
    f = 0.
    trigfunc = iseven(derivative) ? sin : cos
    sign = iseven(div(derivative, 2)) ? 1. : -1.
    for n=1:U.N
        f += sign*U.ξ[n]*U.k[n]^derivative*trigfunc(U.k[n]*x)
    end
    return f
end

function localextremum(U::OpticalPotential, x0::Float64, eps::Float64=1e-5, dxmax::Float64=1., maxiterations=100)
    function f(x::Float64)
        a = evaluate_opticalpotential(U, 1, x)
        b = evaluate_opticalpotential(U, 2, x)
        abs(b)<0.1 ? NaN : -a/b
    end
    fx = 1
    i = 0
    while abs(fx)>eps
        fx = f(x0)
        i += 1
        if i>maxiterations || isnan(fx)
            return NaN
        end
        x0 = x0 + fx
    end
    return x0
end

function extrema(U::OpticalPotential)
    N = 10*prod(U.k)/(2*pi)^U.N
    eps = 1e-3/N
    solutions = Float64[]
    for x0=linspace(0.,1.,N)
        x = localextremum(U, x0, eps)
        if (!isnan(x)) && (0.<=x<1.) && (!any([abs(y-x)<eps for y=solutions]))
            push!(solutions, x)
        end
    end
    return solutions
end

minima(U::OpticalPotential) = filter(x->evaluate_opticalpotential(U,2,x)>0, extrema(U))
maxima(U::OpticalPotential) = filter(x->evaluate_opticalpotential(U,2,x)<0, extrema(U))
global_minimum(U::OpticalPotential) = (P = extrema(U); length(P)==0 ? NaN : P[indmin(evaluate_opticalpotential(U,0,P))])
global_maximum(U::OpticalPotential) = (P = extrema(U); length(P)==0 ? NaN : P[indmax(evaluate_opticalpotential(U,0,P))])

end #module


text:
using Revise
using DataFrames, BenchmarkTools, LinearAlgebra, Distributions, StatsBase, Statistics, Random
Random.seed!(11)

#set up a siimilar regression situation
function crosssection(N,K; loud::Bool=true)

  #generate our data
  ϵ = rand(Laplace(0,1), N)
  Ξ = hcat((k->rand(Laplace(0,2),N)).(1:K)...)
  þ = 1 ./ (þᵢ->abs(þᵢ) ≤ 10^-5 ? 10.0^-5 * sign(þᵢ) : þᵢ).(rand(Chisq(1), N)./4)
  #þ .*= rand(Bernoulli(),N)
  X = Ξ .* þ
  #abstotX = sum(abs.(X), dims=1)

  function normalizemat!(X)
    X .-= mean(X, dims=1)
    abstot = sum(abs.(X), dims=1)
    #abstot₊ = sum((x->x≥0 ? x : 0.).(X), dims=1)
    #abstot₋ = sum((x->x<0 ? -x : 0.).(X), dims=1)
    #X .= hcat(((t)-> (x->x≥0 ? x/(abstot₊[t[1]]) : x/(abstot₋[t[1]])).(t[2])).(enumerate(eachcol(X)))...)
    X ./= abstot .* 2
  end
  #normalize the unobservable actual
  normalizemat!(X)

  #normalize the observable
  normalizemat!(Ξ)

  #σx = std(X, dims=1)
  #X ./= σx
  loud && display(abstotX)
  #isplay(X)
  X = hcat(ones(N),X)
  β = [50.0; (Float64).(1:K) .* 300  ...;]

  y = X*β  .+ ϵ

  #display(β)
  function err(Xqr,y,b)
    u = y.-X*Xb
    Λ = u .* u
    RinvQp = Xqr.R\I * Xqr.Q'
    Σ = RinvQp * Diagonal(Λ) * (RinvQp')
    return diag(Σ).^0.5
  end

  #errchk(X,y,b) = diag(((X'*X)\I) * X'* Diagonal((y.-X*Xb) .* (y.-X*Xb)) * X * ((X'*X)\I)).^0.5

  loud && @info "Test non-confounded version (β=$β)"
  Xqr = qr(X)
  Xrinv = Xqr.R\I
  Xb =  Xrinv * Xqr.Q' * y
  eX = (y.-X*Xb)
  errX = err(Xqr,y,Xb)#diag(eX' * eX ./ (N-K-1) .* (Xqr.R'*Xqr.R)\I) .^ 0.5
  #errXchk = errchk(X,y,Xb)

  loud && @info "Test observed version"
  Ξ1 =hcat(ones(N),Ξ)
  Ξqr = qr(Ξ1)
  Ξrinv = Ξqr.R\I
  Ξb =  Ξrinv * Ξqr.Q' * y
  errΞ = err(Ξqr,y,Ξb)
  loud && @info "Non-confounded version: $(Xb). ErrX = $errX"

  ΣΞbabs = sum(abs.(Ξ1 .* Ξb'), dims=1)
  loud && @info "Results: Raw b: $Ξb, raw error: $errΞ \n"
  loud && @info "Abs total: $(ΣΞbabs)"

  #=@info "Test confounded version"
  Xqr = qr(X)
  Xrinv = Xqr.R\I
  Xb =  Xrinv * Xqr.Q' * y
  @info "Xb: $(Xb) (vs non-confounded of $(Ξb))"
  ΣXb = vec(sum(, dims=1))=#
  return (Ξb=Ξb, ΣΞbabs=ΣΞbabs)
end

tups = vcat((crosssection(2_000,2, loud=false) for i in 1:1000)...)
#display(tups)
Ξb, As = vcat((t->t.Ξb').(tups)...), vcat((t->t.ΣΞbabs).(tups)...)
display(Ξb)
@info "mean raw: $(mean(Ξb, dims=1))"
@info "mean: $(mean(As, dims=1))"
@info "std: $(std(As, dims=1))"

#crosssection(2_000,2)

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
@testset "Formatters" begin
    rec = DefaultRecord("Logger.example", "info", 20, "blah")

    # NOTE: This check might not be needed as of 0.6 because I can't find a
    # condition where the stacktrace is empty.
    no_lookup = DefaultRecord(
        Memento.Attribute(rec.date), Memento.Attribute(rec.level),
        Memento.Attribute(rec.levelnum), Memento.Attribute(rec.msg),
        Memento.Attribute(rec.name), Memento.Attribute(rec.pid),
        Memento.Attribute(nothing), Memento.Attribute(rec.stacktrace)
    )

    @testset "DefaultFormatter" begin
        fmt = DefaultFormatter("{lookup}|{msg}|{stacktrace}")
        result = Memento.format(fmt, rec)
        parts = split(result, "|")
        @test length(parts) == 3
        @test parts[2] == "blah"
        @test length(parts[1]) > 0
        @test length(parts[3]) > 0
        @test occursin("formatters", parts[3])
        @test !occursin("get_trace", parts[3])
        @test !occursin("DefaultRecord", parts[3])
        @test !occursin("get", parts[3])

        nl_result = Memento.format(fmt, no_lookup)
        @test startswith(nl_result, "<nothing>")

        @testset "Date handling" begin
            withenv("TZ" => "UTC+02") do
                local_date = ZonedDateTime(2012, 1, 1, 3, 1, 1, 123, localzone())
                local_datestr = "2012-01-01 03:01:01"
                utc_date = astimezone(local_date, tz"UTC")
                utc_datestr = Dates.format(utc_date, "yyyy-mm-dd HH:MM:SS")

                fmt = DefaultFormatter("{date}")
                tz_rec = SimpleRecord("info", "blah", local_date)
                @test Memento.format(fmt, tz_rec) == local_datestr

                utc_rec = SimpleRecord("info", "blah", utc_date)
                @test Memento.format(fmt, utc_rec) == local_datestr

                utc_fmt = DefaultFormatter("{date}", tz"UTC")
                @test Memento.format(utc_fmt, tz_rec) == utc_datestr
                @test Memento.format(utc_fmt, utc_rec) == utc_datestr

                date = DateTime(2012, 1, 1, 3, 1, 1)
                notz_rec = SimpleRecord("info", "blah", date)
                @test Memento.format(fmt, notz_rec) == local_datestr

                ts = 1531949014
                ts_rec = SimpleRecord("info", "blah", ts)
                @test Memento.format(fmt, ts_rec) == string(ts)
            end
        end
    end

    @testset "DictFormatter" begin
        fmt = DictFormatter()
        result = Memento.format(fmt, rec)
        for key in [:date, :name, :level, :lookup, :stacktrace, :msg]
            @test occursin(string(key), result)
        end

        @test occursin("blah", result)

        aliases = Dict(
            :logger => :name,
            :level => :level,
            :timestamp => :date,
            :location => :lookup,
            :message => :msg,
            :process_id => :pid,
        )

        fmt2 = DictFormatter(aliases)
        result = Memento.format(fmt2, rec)

        for key in [:location, :message, :timestamp, :process_id]
            @test occursin(string(key), result)
            @test !occursin(string(aliases[key]), result)
        end

        @test occursin("blah", result)

        nl_result = Memento.format(fmt, no_lookup)
        @test occursin("<nothing>", nl_result)
    end
end


text:
# ------------------------------------------------------------------------------------------#
# constants (from bionumbers)       units
# ------------------------------------------------------------------------------------------#
cell_diameter = 12                  # mum
number_of_rnapII = 75000            # copies/cells
number_of_ribosome = 1e6            # copies/cells
mRNA_half_life_TF = 2               # hrs
protein_half_life = 10              # hrs
doubling_time_cell = 19.5           # hrs
max_translation_rate = 5            # aa/sec
max_transcription_rate = 6.0        # nt/sec
average_transcript_length = 15000   # nt
average_protein_length = 5000       # aa
fraction_nucleus = 0.49             # dimensionless
av_number = 6.02e23                 # number/mol
avg_gene_number = 2                 # number of copies of a gene
# ------------------------------------------------------------------------------------------#
# Calculate constants using bionumber values
# ------------------------------------------------------------------------------------------#
# Calculate the volume (convert to L)
V = ((1-fraction_nucleus)*(1/6)*(3.14159)*(cell_diameter)^3)*(1e-15)       # volume

# Calculate the rnapII_concentration and ribosome_concentration
rnapII_concentration = number_of_rnapII*(1/av_number)*(1/V)*1e9                   # nM
ribosome_concentration = number_of_ribosome*(1/av_number)*(1/V)*1e9               # nM

# degrdation rate constants -
degradation_constant_mRNA = -(1/mRNA_half_life_TF)*log(0.5)                       # hr^-1
degradation_constant_protein = -(1/protein_half_life)*log(0.5)                    # hr^-1

# kcats for transcription and translation -
kcat_transcription = max_transcription_rate*(3600/average_transcript_length)      # hr^-1
kcat_translation = max_translation_rate*(3600/average_protein_length)             # hr^-1

# Maximum specific growth rate -
maximum_specific_growth_rate = (1/doubling_time_cell)*log(2)                      # hr^-1

# What is the average gene concentration -
avg_gene_concentration = avg_gene_number*(1/av_number)*(1/V)*1e9                  # nM

# How fast do my cells die?
death_rate_constant = 0.2*maximum_specific_growth_rate                            # hr^-1

# Saturation constants for translation and trascription -
saturation_transcription = 4600*(1/av_number)*(1/V)*1e9                           # nM
saturation_translation = 100000*(1/av_number)*(1/V)*1e9                           # nM
# -------------------------------------------------------------------------------------------#


text:
@inline function get_interpolation(wp::Interpolations.Extrapolation{T, 3}, pt::CartesianPoint{T}, S::Val{:cartesian})::T where {T <: SSDFloat}
    return wp(pt.x, pt.y, pt.z)::T
end
@inline function get_interpolation(wp::Interpolations.Extrapolation{T, 3}, pt::CartesianPoint{T}, S::Val{:cylindrical})::T where {T <: SSDFloat}
    p::CylindricalPoint{T} = CylindricalPoint(pt)
    return wp(p.r, p.φ, p.z)::T
end

function add_signal!(signal::AbstractVector{T}, timestamps::AbstractVector{T}, path::Vector{CartesianPoint{T}}, pathtimestamps::AbstractVector{T}, charge::T,
                        wp::Interpolations.Extrapolation{T, 3}, S::Union{Val{:cylindrical}, Val{:cartesian}})::Nothing where {T <: SSDFloat}
    tmp_signal = Vector{T}(undef, length(pathtimestamps))
    @inbounds for i in eachindex(tmp_signal)
        tmp_signal[i] = get_interpolation(wp, path[i], S)::T * charge
    end
    itp = interpolate( (pathtimestamps,), tmp_signal, Gridded(Linear()))
    t_max::T = last(pathtimestamps)
    i_max::Int = searchsortedlast(timestamps, t_max)
    signal[1:i_max] += itp(timestamps[1:i_max])
    if length(signal) > i_max
        signal[i_max+1:end] .+= tmp_signal[end]
    end
    nothing
end

function add_signal!(signal::AbstractVector{T}, timestamps::AbstractVector{TT}, path::EHDriftPath{T, TT}, charge::T, wp::Interpolations.Extrapolation{T, 3}, S::Union{Val{:cylindrical}, Val{:cartesian}})::Nothing where {T <: SSDFloat, TT}
    add_signal!(signal, timestamps, path.e_path, path.timestamps_e, -charge, wp, S) # electrons induce negative charge
    add_signal!(signal, timestamps, path.h_path, path.timestamps_h,  charge, wp, S)
    nothing
end

function add_signal!(signal::AbstractVector{T}, timestamps::AbstractVector{<:RealQuantity}, paths::Vector{<:EHDriftPath{T}}, charges::Vector{T}, wp::Interpolations.Extrapolation{T, 3}, S::Union{Val{:cylindrical}, Val{:cartesian}})::Nothing where {T <: SSDFloat}
    for ipath in eachindex(paths)
        add_signal!(signal, timestamps, paths[ipath], charges[ipath], wp, S)
    end
    nothing
end

# function get_signal(path::Vector{CartesianPoint{T}}, charge::T, wp::Interpolations.Extrapolation{T, 3}, S::Union{Val{:cylindrical}, Val{:cartesian}})::Vector{T} where {T <: SSDFloat}
#     signal::Vector{T} = zeros(T, length(path))
#     add_signal!(signal, path, charge, wp, S)
#     return signal
# end


text: include("../iml/julia/IML.jl")


text:
using Hecke

function ideals_with_pp_norm(lp::Array{NfOrdIdl, 1}, k::Int)
  l = [degree(x) for x= lp]
#  println("pp with $l and $k")
  #need sum([e[i]*l[i] == k, e[i] >= 0])
  C = cartesian_product_iterator([0:div(k, l[j]) for j = 1:length(l)], inplace = true)
  r = [prod(lp[i]^c[i] for i=1:length(l)) for c = C if sum(c[i]*l[i] for i=1:length(l)) == k]
#  println("r: $r")
  return r
end

function ideals_with_norm(i::fmpz, M::NfOrd)
  @assert M.ismaximal == 1
  if isone(i)
    return [ideal(M, 1)]
  end
  lf = factor(i)
  lp = [ideals_with_pp_norm([x[1] for x = prime_decomposition(M, k)], v) for (k, v) = lf.fac]
#  println(lp)
  return [prod(lp[i][x[i]] for i=1:length(lf.fac)) for x = cartesian_product_iterator([1:length(lp[i]) for i=1:length(lp)], inplace = true)]
end

function orbit_reps(I::Array{NfOrdIdl, 1}, s::Hecke.NfToNfMor)
  O = Set([I[1], Hecke.induce_image(I[1], s)])
  R = [I[1]]
  for i=I
    if i in O
      continue
    end
    push!(R, i)
    push!(O, i)
    push!(O, Hecke.induce_image(i, s))
  end
  return R
end

#Note: this is not optimised, but hopefully correct.
#if you need more, analyse Hasse...

function induce_action(phi::Hecke.NfToNfMor, mR::T) where T <: Map{GrpAbFinGen,
Hecke.FacElemMon{Hecke.NfOrdIdlSet}}
#function induce_action(phi::Hecke.NfToNfMor, mR::Hecke.MapRayClassGrpFacElem{Hecke.GrpAbFinGen})
  lp, x = Hecke.find_gens(
        Hecke.MapFromFunc(x->preimage(mR, FacElem(Dict(x=>1))),
                          base_ring(codomain(mR)),
                          domain(mR)),
        PrimesSet(100, -1))
  return hom(x, GrpAbFinGenElem[preimage(mR, Hecke.induce_image(p, phi)) for p = lp])
end

function s3_with_discriminant(I::NfOrdIdl)
  lI = factor(I)
  kx = PolynomialRing(order(I).nf)[1]
  #we need deccompositions I = df^2
  #and f is squarefree - exccept at 3
  #there can only be wild ramification at primes dividing the degree
  #similarly: d is squarefree outside 2...
  all_poss = Array{Tuple{NfOrdIdl, NfOrdIdl}, 1}()
  l23 = []
  f = ideal(order(I), 1)
  d = ideal(order(I), 1)
  for (i, e) = lI
    if minimum(i) == 2
      l = [(i^e, ideal(order(i), 1))]
      if e>=2
        push!(l, (i^(e-2), i))
      end
      push!(l23, l)
      continue
    end
    if minimum(i) == 3
      if iseven(e)
        f *= i^div(e, 2)
      end
      if isodd(e)
        d *= i
        f *= i^div(e-1, 2)
      end
      continue
    end

    if e==1
      d *= i
    elseif e==2
      if norm(f) % 3 != 1
        println("no 3 in the prime, cannot work")
        return []
      end
      f *= i
    elseif e==3
      if norm(f) % 3 != 1
        println("no 3 in the prime, cannot work")
        return []
      end
      f *= i
      d *= i
    else
      println("expo too large at wrong prime, no S3-disc")
      return []
    end
  end

  res = []
  println("23: $l23")
  C = cartesian_product_iterator([1:length(i) for i in l23], inplace = true)
  #Tommy: use Base.product(l23,...) or similar....
  for c = C
    if length(c) == 0
      D = d
      F = f
    else
      D = d*prod(l23[i][c[i]][1] for i=1:length(c))
      F = f*prod(l23[i][c[i]][2] for i=1:length(c))
    end
    println("need to try $D and $F as conductors")
    #all quadratics with conductor D:
    r, mr = ray_class_group(D, n_quo = 2)

    for s in index_p_subgroups(r, fmpz(2), (A,x) -> quo(A, x)[2])
      a = ray_class_field(mr*pseudo_inv(s))
#      println(a, " with cond ", conductor(a))
      if conductor(a)[1] != D
        continue
      end
      K = number_field(a)
      println("have K: $K")
      @assert length(gens(K)) == 1
      Kr = K
      @assert degree(Kr) == 2
#      @assert Hecke.isradical_extension(Kr)
      conj = Hecke.rel_auto(a.cyc[1])
      Kr = a.cyc[1].A
      k = base_ring(Kr)
      Ka, m1 = absolute_simple_field(Kr)
      m2 = restrict(inv(m1), k)
      sigma = hom(Ka, Ka, m1\(conj(m1(gen(Ka)))))
      #m1: Kr -> Ka, m2: base_ring(Kr) -> Ka
      M = lll(maximal_order(Ka))
      FF = ideal(M, F.gen_one, M(m2(k(F.gen_two))))
      R, mR = ray_class_group(FF, n_quo = 3)
      if order(R) == 1
        println("RCG empty")
        continue
      end

      sigma_R = induce_action(sigma, mR)
#      println(sigma_R)

      for S = Hecke.stable_subgroups(R, [3], [sigma_R], op=quo)
        @assert order(S[1]) == 3
        s, ms = snf(S[1])
        if ms\(S[2](sigma_R(S[2]\(ms(s[1]))))) == s[1]
          #TODO: factor out the part with trivial action
          # ie. kern(sigma_R-I)
          println("action is trivial, no S3")
          continue
        end
        A = ray_class_field(mR*pseudo_inv(S[2]))
#        println(A, " with cond ", conductor(A))
        if conductor(A)[1] != FF
          println("wrong conductor")
          continue
        end
        B = number_field(A)
        Tau = Hecke.extend_aut(A, sigma)
        g = gens(B)[1]
        g += Tau(g)
        if iszero(g)
          g = gens(B)[1]
          g *= Tau(g)
        end
        mg = minpoly(g)
        @assert degree(mg) == 3
        mg = kx([coeff(m1\coeff(mg, i), 0) for i=0:4])
        push!(res, mg)
      end
    end
  end
  return res
end

function Gunter_Qi(r::Range, pref="Qi.new")
  Qt, t = FlintQQ["t"]
  k, a = number_field(t^2+1, "k.1")
  s = Hecke.NfToNfMor(k, k, -a)
  M = maximal_order(k)
  kx, x = k["kx.1"]
  f = open("/tmp/$pref", "a")
  g = open("/tmp/$pref.err", "a")
  err_cnt = 0

  for i = r
    try
      I = ideals_with_norm(fmpz(i), M)
      if length(I)==0
        continue
      end
      I = NfOrdIdl[I[x] for x = linearindices(I)]
      R = orbit_reps(I, s)
      for r in R
        z = s3_with_discriminant(r)
        if length(z)==0
          continue
        end
        for j=z
          println(f, "<$i, $(kx([coeff(j, jj) for jj=0:3]))>, ")
        end
      end
    catch e
      rethrow(e)
      err_cnt += 1
      println("i: $i")
      println(e)
      println(g, "i: $i")
      println(g, e)
      if err_cnt > 10
        rethrow(e)
      end
    end
  end
  close(f)
end


text:
"""
    getdescendant(node, idx)

Obtain a node from a tree by indexing each level of the tree with the elements of `idx`.

This function is defined for all trees regardless of whether they have the [`IndexedChildren`](@ref).
This is because a tree without [`IndexedChildren`](@ref) might have special cases in which all children
are indexable, a prominent example being `Array` which may not have indexable sub-trees (e.g. an
array containing a Dict) but there are common special cases in which array trees are fully indexable
(e.g. a tree in which every non-leaf node is an array).

The elements of `idx` can be any argument to `getindex`, not necessarily integers.  For example,
`getdescendant(Dict("a"=>1), ("a",))` returns `1`.

Note that this is a separate concept from indexed trees which by default do not have `IndexedChildren()`,
see [`IndexNode`](@ref).

## Example
```julia
v = [1, [2, [3, 4]]]

getdescendant(v, (2, 2, 1)) == 3
```
"""
function getdescendant(node, idx)
    n = node
    for j ∈ idx
        n = children(n)[j]
    end
    n
end

"""
    childindices(tree, node_index)

Get the indices of the children of the node of tree `tree` specified by `node_index`.

To be consistent with [`children`](@ref), by default this returns an empty tuple.

**REQUIRED** for indexed trees:  Indexed trees, i.e. trees that do not implement [`children`](@ref) must implement
this function.
"""
childindices(tree, idx) = ()

"""
    parentindex(tree, node_index)

Get the index of the parent of the node of tree `tree` specified by `node_index`.

Nodes that have no parent (i.e. the root node) should return `nothing`.

**OPTIONAL**: Indexed trees with the [`StoredParents`](@ref) trait must implement this.
"""
function parentindex end

"""
    nextsiblingindex(tree, node_index)

Get the index of the next sibling of the node of tree `tree` specified by `node_index`.

Nodes which have no next sibling should return `nothing`.

**OPTIONAL**: Indexed trees with the [`StoredSiblings`](@ref) trait must implement this.
"""
function nextsiblingindex end

"""
    prevsiblingindex(tree, node_index)

Get the index of the previous sibling of the node of tree `tree` specified by `node_index`.

Nodes which have no previous sibling should return `nothing`.

**OPTIONAL**: Indexed trees that have [`StoredSiblings`](@ref) can implement this, but no built-in tree algorithms
require it.
"""
function prevsiblingindex end

"""
    rootindex(tree)

Get the root index of the indexed tree `tree`.

**OPTIONAL**: The single-argument constructor for [`IndexNode`](@ref) requires this, but it is not required for
any built-in tree algorithms.
"""
function rootindex end

"""
    IndexNode{T,I}

The node of a tree which implements the indexed tree interface.  Such a tree consists of an object `tree` from
which nodes can be obtained with the two-argument method of [`nodevalue`](@ref) which by default calls `getindex`.

An `IndexNode` implements the tree interface, and can be thought of an adapter from an object that implements the
indexed tree interface to one that implements the tree interface.

`IndexNode` do not store the value associated with the node but can obtain it by calling [`nodevalue`](@ref).

## Constructors
```julia
IndexNode(tree, node_index)

IndexNode(tree) = IndexNode(tree, rootindex(tree))  # one-argument constructor requires `rootindex`
```

Here `tree` is an object which stores or can obtain information for the entire tree structure, and `node_index`
is the index of the node for which `node_index` is being constructed.
"""
struct IndexNode{T,I}
    tree::T
    index::I
end

IndexNode(tree) = IndexNode(tree, rootindex(tree))

ParentLinks(::Type{<:IndexNode{N,T}}) where {N,T} = ParentLinks(T)
SiblingLinks(::Type{<:IndexNode{N,T}}) where {N,T} = SiblingLinks(T)

nodevalue(idx::IndexNode) = nodevalue(idx.tree, idx.index)

children(idx::IndexNode) = Iterators.map(c -> IndexNode(idx.tree, c), childindices(idx.tree, idx.index))

function parent(idx::IndexNode)
    pidx = parentindex(idx.tree, idx.index)
    isnothing(pidx) ? nothing : IndexNode(idx.tree, pidx)
end

function nextsibling(idx::IndexNode)
    sidx = nextsiblingindex(idx.tree, idx.index)
    isnothing(sidx) ? nothing : IndexNode(idx.tree, sidx)
end

function prevsibling(rn::IndexNode)
    sidx = prevsiblingindex(idx.tree, idx.index)
    isnothing(sidx) ? nothing : IndexNode(idx.tree, sidx)
end


text:
module FrameFunDomains

using DomainSets, GridArrays, StaticArrays

import Base: show, isapprox, in, <, <=, >, >=
import DomainSets: indomain
import GridArrays: boundingbox, indomain_broadcast
include("extensions.jl")
include("fourierdomains.jl")
include("fractals.jl")
include("atomium.jl")
include("characteristic.jl")

include("polardomain.jl")
end


text:
function show(io::IO, emod::ERA5Hourly)
    print(
		io,
		"The ERA5Hourly Module has the following properties:\n",
		"    Module ID      (modID) : ", emod.modID, '\n',
		"    Data Directory (eroot) : ", emod.eroot, '\n',
		"    Date Begin     (dtbeg) : ", emod.dtbeg, '\n',
		"    Date End       (dtend) : ", emod.dtend, '\n',
	)
end

function show(io::IO, emod::ERA5Monthly)
    print(
		io,
		"The ERA5Monthly Module has the following properties:\n",
		"    Module ID      (modID) : ", emod.modID, '\n',
		"    Data Directory (eroot) : ", emod.eroot, '\n',
		"    Date Begin     (dtbeg) : ", emod.dtbeg, '\n',
		"    Date End       (dtend) : ", emod.dtend, '\n',
	)
end

function show(io::IO, emod::ERA5MonthlyHour)
    print(
		io,
		"The ERA5MonthlyHour Module has the following properties:\n",
		"    Module ID      (modID) : ", emod.modID, '\n',
		"    Data Directory (eroot) : ", emod.eroot, '\n',
		"    Date Begin     (dtbeg) : ", emod.dtbeg, '\n',
		"    Date End       (dtend) : ", emod.dtend, '\n',
		"    Hours-of-Day   (hours) : ", emod.hours, '\n',
	)
end


text:
export Instruction, unsafe_delete!,
       metadata, metadata!,
       opcode,
       predicate_int, predicate_real

# forward definition of Instruction in src/core/value/constant.jl
register(Instruction, API.LLVMInstructionValueKind)

const instruction_opcodes = Vector{Type}(fill(Nothing, typemax(API.LLVMOpcode)+1))
function identify(::Type{Instruction}, ref::API.LLVMValueRef)
    opcode = API.LLVMGetInstructionOpcode(ref)
    typ = @inbounds instruction_opcodes[opcode+1]
    typ === Nothing && error("Unknown type opcode $opcode")
    return typ
end
function register(T::Type{<:Instruction}, opcode::API.LLVMOpcode)
    instruction_opcodes[opcode+1] = T
end

function refcheck(::Type{T}, ref::API.LLVMValueRef) where T<:Instruction
    ref==C_NULL && throw(UndefRefError())
    if Base.JLOptions().debug_level >= 2
        T′ = identify(Instruction, ref)
        if T != T′
            error("invalid conversion of $T′ instruction reference to $T")
        end
    end
end

# Construct a concretely typed instruction object from an abstract value ref
function Instruction(ref::API.LLVMValueRef)
    ref == C_NULL && throw(UndefRefError())
    T = identify(Instruction, ref)
    return T(ref)
end

Instruction(inst::Instruction) =
    Instruction(API.LLVMInstructionClone(inst))

unsafe_delete!(::BasicBlock, inst::Instruction) =
    API.LLVMInstructionEraseFromParent(inst)
Base.delete!(::BasicBlock, inst::Instruction) =
    API.LLVMInstructionRemoveFromParent(inst)

parent(inst::Instruction) =
    BasicBlock(API.LLVMGetInstructionParent(inst))

opcode(inst::Instruction) = API.LLVMGetInstructionOpcode(inst)

predicate_int(inst::Instruction) = API.LLVMGetICmpPredicate(inst)
predicate_real(inst::Instruction) = API.LLVMGetFCmpPredicate(inst)

## metadata iteration
# TODO: doesn't actually iterate, since we can't list the available keys

@enum(MD, MD_dbg = 0,
          MD_tbaa = 1,
          MD_prof = 2,
          MD_fpmath = 3,
          MD_range = 4,
          MD_tbaa_struct = 5,
          MD_invariant_load = 6,
          MD_alias_scope = 7,
          MD_noalias = 8,
          MD_nontemporal = 9,
          MD_mem_parallel_loop_access = 10,
          MD_nonnull = 11,
          MD_dereferenceable = 12,
          MD_dereferenceable_or_null = 13,
          MD_make_implicit = 14,
          MD_unpredictable = 15,
          MD_invariant_group = 16,
          MD_align = 17,
          MD_loop = 18,
          MD_type = 19,
          MD_section_prefix = 20,
          MD_absolute_symbol = 21,
          MD_associated = 22)

export InstructionMetadataDict

struct InstructionMetadataDict <: AbstractDict{MD,MetadataAsValue}
    inst::Instruction
end

metadata(inst::Instruction) = InstructionMetadataDict(inst)

Base.isempty(md::InstructionMetadataDict) =
  !convert(Core.Bool, API.LLVMHasMetadata(md.inst))

Base.haskey(md::InstructionMetadataDict, kind::MD) =
  API.LLVMGetMetadata(md.inst, kind) != C_NULL

function Base.getindex(md::InstructionMetadataDict, kind::MD)
    objref = API.LLVMGetMetadata(md.inst, kind)
    objref == C_NULL && throw(KeyError(kind))
    return Metadata(MetadataAsValue(objref))
  end

Base.setindex!(md::InstructionMetadataDict, node::Metadata, kind::MD) =
    API.LLVMSetMetadata(md.inst, kind, Value(node; ctx=context(md.inst)))

Base.delete!(md::InstructionMetadataDict, kind::MD) =
    API.LLVMSetMetadata(md.inst, kind, C_NULL)

## instruction types

const opcodes = [:Ret, :Br, :Switch, :IndirectBr, :Invoke, :Unreachable, :CallBr, :FNeg,
                 :Add, :FAdd, :Sub, :FSub, :Mul, :FMul, :UDiv, :SDiv, :FDiv, :URem, :SRem,
                 :FRem, :Shl, :LShr, :AShr, :And, :Or, :Xor, :Alloca, :Load, :Store,
                 :GetElementPtr, :Trunc, :ZExt, :SExt, :FPToUI, :FPToSI, :UIToFP, :SIToFP,
                 :FPTrunc, :FPExt, :PtrToInt, :IntToPtr, :BitCast, :AddrSpaceCast, :ICmp,
                 :FCmp, :PHI, :Call, :Select, :UserOp1, :UserOp2, :VAArg, :ExtractElement,
                 :InsertElement, :ShuffleVector, :ExtractValue, :InsertValue, :Freeze,
                 :Fence, :AtomicCmpXchg, :AtomicRMW, :Resume, :LandingPad, :CleanupRet,
                 :CatchRet, :CatchPad, :CleanupPad, :CatchSwitch]

for op in opcodes
    typename = Symbol(op, :Inst)
    enum = Symbol(:LLVM, op)
    @eval begin
        @checked struct $typename <: Instruction
            ref::API.LLVMValueRef
        end
        register($typename, API.$enum)
    end
end

## call sites and invocations

# TODO: add this to the actual type hierarchy
const CallBase = Union{CallBrInst, CallInst, InvokeInst}

export callconv, callconv!,
       istailcall, tailcall!,
       called_value, arguments,
       OperandBundleUse, OperandBundleDef, operand_bundles

callconv(inst::CallBase) = API.LLVMGetInstructionCallConv(inst)
callconv!(inst::CallBase, cc) =
    API.LLVMSetInstructionCallConv(inst, cc)

istailcall(inst::CallBase) = convert(Core.Bool, API.LLVMIsTailCall(inst))
tailcall!(inst::CallBase, bool) = API.LLVMSetTailCall(inst, convert(Bool, bool))

called_value(inst::CallBase) = Value(API.LLVMGetCalledValue(inst))

function arguments(inst::CallBase)
    nargs = API.LLVMGetNumArgOperands(inst)
    operands(inst)[1:nargs]
end

# operand bundles

# XXX: these objects are just C structures, whose lifetime isn't tied to LLVM IR structures.
#      that means we need to create a copy when passing them to Julia, necessitating disposal.

@checked mutable struct OperandBundleUse
    ref::API.LLVMOperandBundleUseRef
end
Base.unsafe_convert(::Type{API.LLVMOperandBundleUseRef}, bundle::OperandBundleUse) =
    bundle.ref

struct OperandBundleIterator <: AbstractVector{OperandBundleUse}
    inst::Instruction
end

operand_bundles(inst::Instruction) = OperandBundleIterator(inst)

Base.size(iter::OperandBundleIterator) = (API.LLVMGetNumOperandBundles(iter.inst),)

Base.IndexStyle(::OperandBundleIterator) = IndexLinear()

function Base.getindex(iter::OperandBundleIterator, i::Int)
    @boundscheck 1 <= i <= length(iter) || throw(BoundsError(iter, i))
    bundle = OperandBundleUse(API.LLVMGetOperandBundle(iter.inst, i-1))
    finalizer(bundle) do obj
        API.LLVMDisposeOperandBundleUse(obj)
    end
end

tag_id(bundle::OperandBundleUse) = API.LLVMGetOperandBundleTagID(bundle)

function tag_name(bundle::OperandBundleUse)
    len = Ref{Cuint}()
    data = API.LLVMGetOperandBundleUseTagName(bundle, len)
    unsafe_string(convert(Ptr{Int8}, data), len[])
end

function inputs(bundle::OperandBundleUse)
    nvals = API.LLVMGetOperandBundleUseNumInputs(bundle)
    vals = Vector{API.LLVMValueRef}(undef, nvals)
    API.LLVMGetOperandBundleUseInputs(bundle, vals)
    return [Value(val) for val in vals]
end

@checked mutable struct OperandBundleDef
    ref::API.LLVMOperandBundleDefRef
end
Base.unsafe_convert(::Type{API.LLVMOperandBundleDefRef}, bundle::OperandBundleDef) =
    bundle.ref

function tag_name(bundle::OperandBundleDef)
    len = Ref{Cuint}()
    data = API.LLVMGetOperandBundleDefTag(bundle, len)
    unsafe_string(convert(Ptr{Int8}, data), len[])
end

function OperandBundleDef(bundle_use::OperandBundleUse)
    bundle_def = OperandBundleDef(API.LLVMOperandBundleDefFromUse(bundle_use))
    finalizer(bundle_def) do obj
        API.LLVMDisposeOperandBundleDef(obj)
    end
end

function OperandBundleDef(tag::String, inputs::Vector{<:Value}=Value[])
    bundle = OperandBundleDef(API.LLVMCreateOperandBundleDef(tag, inputs, length(inputs)))
    finalizer(bundle) do obj
        API.LLVMDisposeOperandBundleDef(obj)
    end
end

function inputs(bundle::OperandBundleDef)
    nvals = API.LLVMGetOperandBundleDefNumInputs(bundle)
    vals = Vector{API.LLVMValueRef}(undef, nvals)
    API.LLVMGetOperandBundleDefInputs(bundle, vals)
    return [Value(val) for val in vals]
end

function Base.show(io::IO, bundle::Union{OperandBundleUse,OperandBundleDef})
    # mimic how bundles are rendered in LLVM IR
    print(io, "\"", tag_name(bundle), "\"(")
    join(io, inputs(bundle), ", ")
    print(io, ")")
end

## terminators

export isterminator, isconditional, condition, condition!, default_dest

isterminator(inst::Instruction) = API.LLVMIsATerminatorInst(inst) != C_NULL

isconditional(br::Instruction) = convert(Core.Bool, API.LLVMIsConditional(br))

condition(br::Instruction) =
    Value(API.LLVMGetCondition(br))
condition!(br::Instruction, cond::Value) =
    API.LLVMSetCondition(br, cond)

default_dest(switch::Instruction) =
    BasicBlock(API.LLVMGetSwitchDefaultDest(switch))

# successor iteration

export successors

struct TerminatorSuccessorSet <: AbstractVector{BasicBlock}
    term::Instruction
end

successors(term::Instruction) = TerminatorSuccessorSet(term)

Base.size(iter::TerminatorSuccessorSet) = (API.LLVMGetNumSuccessors(iter.term),)

Base.IndexStyle(::TerminatorSuccessorSet) = IndexLinear()

function Base.getindex(iter::TerminatorSuccessorSet, i::Int)
    @boundscheck 1 <= i <= length(iter) || throw(BoundsError(iter, i))
    return BasicBlock(API.LLVMGetSuccessor(iter.term, i-1))
end

Base.setindex!(iter::TerminatorSuccessorSet, bb::BasicBlock, i::Int) =
    API.LLVMSetSuccessor(iter.term, i-1, bb)

## phi nodes

# incoming iteration

export PhiIncomingSet

struct PhiIncomingSet <: AbstractVector{Tuple{Value,BasicBlock}}
    phi::Instruction
end

incoming(phi::Instruction) = PhiIncomingSet(phi)

Base.size(iter::PhiIncomingSet) = (API.LLVMCountIncoming(iter.phi),)

Base.IndexStyle(::PhiIncomingSet) = IndexLinear()

function Base.getindex(iter::PhiIncomingSet, i::Int)
    @boundscheck 1 <= i <= length(iter) || throw(BoundsError(iter, i))
    return tuple(Value(API.LLVMGetIncomingValue(iter.phi, i-1)),
                       BasicBlock(API.LLVMGetIncomingBlock(iter.phi, i-1)))
end

function Base.append!(iter::PhiIncomingSet, args::Vector{Tuple{V, BasicBlock}} where V <: Value)
    vals, blocks = zip(args...)
    API.LLVMAddIncoming(iter.phi, collect(vals), collect(blocks), length(args))
end

Base.push!(iter::PhiIncomingSet, args::Tuple{<:Value, BasicBlock}) = append!(iter, [args])


text:
include("euler/euler.jl")
using .Calculus: len_collatz_chain

function compute(n::Int)::Int
    result, max_chain, hashmap = 0, 0, Dict(1 => 1)
    for i ∈ 2:n-1
        if len_collatz_chain(i, hashmap) > max_chain result, max_chain = i, hashmap[i] end
    end
    return result
end

max_stars_repo_name: Dynortice/Project-Euler

max_issues_repo_name: Dynortice/Project-Euler

max_forks_repo_name: Dynortice/Project-Euler


text:
# Binary heap

immutable MutableBinaryHeapNode{T}
    value::T
    handle::Int
end

#################################################
#
#   core implementation
#
#################################################

function _heap_bubble_up!{Comp, T}(comp::Comp,
    nodes::Vector{MutableBinaryHeapNode{T}}, nodemap::Vector{Int}, nd_id::Int)

    @inbounds nd = nodes[nd_id]
    v::T = nd.value

    swapped = true  # whether swap happens at last step
    i = nd_id

    while swapped && i > 1  # nd is not root
        p = i >> 1
        @inbounds nd_p = nodes[p]

        if compare(comp, v, nd_p.value)
            # move parent downward
            @inbounds nodes[i] = nd_p
            @inbounds nodemap[nd_p.handle] = i
            i = p
        else
            swapped = false
        end
    end

    if i != nd_id
        nodes[i] = nd
        nodemap[nd.handle] = i
    end
end

function _heap_bubble_down!{Comp, T}(comp::Comp,
    nodes::Vector{MutableBinaryHeapNode{T}}, nodemap::Vector{Int}, nd_id::Int)

    @inbounds nd = nodes[nd_id]
    v::T = nd.value

    n = length(nodes)
    last_parent = n >> 1

    swapped = true
    i = nd_id

    while swapped && i <= last_parent
        il = i << 1

        if il < n   # contains both left and right children
            ir = il + 1

            # determine the better child
            @inbounds nd_l = nodes[il]
            @inbounds nd_r = nodes[ir]

            if compare(comp, nd_r.value, nd_l.value)
                # consider right child
                if compare(comp, nd_r.value, v)
                    @inbounds nodes[i] = nd_r
                    @inbounds nodemap[nd_r.handle] = i
                    i = ir
                else
                    swapped = false
                end
            else
                # consider left child
                if compare(comp, nd_l.value, v)
                    @inbounds nodes[i] = nd_l
                    @inbounds nodemap[nd_l.handle] = i
                    i = il
                else
                    swapped = false
                end
            end

        else  # contains only left child
            nd_l = nodes[il]
            if compare(comp, nd_l.value, v)
                @inbounds nodes[i] = nd_l
                @inbounds nodemap[nd_l.handle] = i
                i = il
            else
                swapped = false
            end
        end
    end

    if i != nd_id
        @inbounds nodes[i] = nd
        @inbounds nodemap[nd.handle] = i
    end
end

function _binary_heap_pop!{Comp,T}(comp::Comp,
    nodes::Vector{MutableBinaryHeapNode{T}}, nodemap::Vector{Int})

    # extract root node
    rt = nodes[1]
    v = rt.value
    @inbounds nodemap[rt.handle] = 0

    if length(nodes) == 1
        # clear
        empty!(nodes)
    else
        # place last node to root
        @inbounds nodes[1] = new_rt = pop!(nodes)
        @inbounds nodemap[new_rt.handle] = 1

        if length(nodes) > 1
            _heap_bubble_down!(comp, nodes, nodemap, 1)
        end
    end
    v
end

function _make_mutable_binary_heap{Comp,T}(comp::Comp, ty::Type{T}, values)
    # make a static binary index tree from a list of values

    n = length(values)
    nodes = Array(MutableBinaryHeapNode{T}, n)
    nodemap = Array(Int, n)

    i::Int = 0
    for v in values
        i += 1
        @inbounds nodes[i] = MutableBinaryHeapNode{T}(v, i)
        @inbounds nodemap[i] = i
    end

    for i = 1 : n
        _heap_bubble_up!(comp, nodes, nodemap, i)
    end
    return nodes, nodemap
end

#################################################
#
#   Binary Heap type and constructors
#
#################################################

type MutableBinaryHeap{VT, Comp} <: AbstractMutableHeap{VT,Int}
    comparer::Comp
    nodes::Vector{MutableBinaryHeapNode{VT}}
    node_map::Vector{Int}

    function MutableBinaryHeap(comp::Comp)
        nodes = Array(MutableBinaryHeapNode{VT}, 0)
        node_map = Array(Int, 0)
        new(comp, nodes, node_map)
    end

    function MutableBinaryHeap(comp::Comp, xs)  # xs is an iterable collection of values
        nodes, node_map = _make_mutable_binary_heap(comp, VT, xs)
        new(comp, nodes, node_map)
    end
end

mutable_binary_minheap{T}(ty::Type{T}) = MutableBinaryHeap{T,LessThan}(LessThan())
mutable_binary_maxheap{T}(ty::Type{T}) = MutableBinaryHeap{T,GreaterThan}(GreaterThan())

mutable_binary_minheap{T}(xs::AbstractVector{T}) = MutableBinaryHeap{T,LessThan}(LessThan(), xs)
mutable_binary_maxheap{T}(xs::AbstractVector{T}) = MutableBinaryHeap{T,GreaterThan}(GreaterThan(), xs)

function show(io::IO, h::MutableBinaryHeap)
    print(io, "MutableBinaryHeap(")
    nodes = h.nodes
    n = length(nodes)
    if n > 0
        print(io, string(nodes[1].value))
        for i = 2 : n
            print(io, ", $(nodes[i].value)")
        end
    end
    print(io, ")")
end

#################################################
#
#   interfaces
#
#################################################

length(h::MutableBinaryHeap) = length(h.nodes)

isempty(h::MutableBinaryHeap) = isempty(h.nodes)

function push!{T}(h::MutableBinaryHeap{T}, v::T)
    nodes = h.nodes
    nodemap = h.node_map
    i = length(nodemap) + 1
    nd_id = length(nodes) + 1
    push!(nodes, MutableBinaryHeapNode(v, i))
    push!(nodemap, nd_id)
    _heap_bubble_up!(h.comparer, nodes, nodemap, nd_id)
    i
end

top(h::MutableBinaryHeap) = h.nodes[1].value

pop!{T}(h::MutableBinaryHeap{T}) = _binary_heap_pop!(h.comparer, h.nodes, h.node_map)

function update!{T}(h::MutableBinaryHeap{T}, i::Int, v::T)
    nodes = h.nodes
    nodemap = h.node_map
    comp = h.comparer

    nd_id = nodemap[i]
    v0 = nodes[nd_id].value
    nodes[nd_id] = MutableBinaryHeapNode(v, i)
    if compare(comp, v, v0)
        _heap_bubble_up!(comp, nodes, nodemap, nd_id)
    else
        _heap_bubble_down!(comp, nodes, nodemap, nd_id)
    end
end

max_stars_repo_name: ambujone/DIDebug

max_issues_repo_name: ambujone/DIDebug

max_forks_repo_name: ambujone/DIDebug