text:
#
# Trailer path planning library
#
# author: Atsushi Sakai(@Atsushi_twi)
#

module trailerlib

using PyPlot
using NearestNeighbors

# Vehicle parameter
const WB = 3.7  #[m] wheel base: rear to front steer
const LT = 8.0 #[m] rear to trailer wheel
const W = 2.6 #[m] width of vehicle
const LF = 4.5 #[m] distance from rear to vehicle front end of vehicle
const LB = 1.0 #[m] distance from rear to vehicle back end of vehicle
const LTF = 1.0 #[m] distance from rear to vehicle front end of trailer
const LTB = 9.0 #[m] distance from rear to vehicle back end of trailer
const MAX_STEER = 0.6 #[rad] maximum steering angle
const TR = 0.5 # Tyre radius [m] for plot
const TW = 1.0 # Tyre width [m] for plot

# for collision check
const WBUBBLE_DIST = 3.5 #distance from rear and the center of whole bubble
const WBUBBLE_R = 10.0 # whole bubble radius
const B = 4.45 # distance from rear to vehicle back end
const C = 11.54 # distance from rear to vehicle front end
const I = 8.55 # width of vehicle
const VRX = [C, C, -B, -B, C ]
const VRY = [-I/2.0, I/2.0, I/2.0, -I/2.0, -I/2.0]

function check_collision(x::Array{Float64},
                         y::Array{Float64},
                         yaw::Array{Float64},
                         kdtree::KDTree,
                         ox::Array{Float64},
                         oy::Array{Float64},
                         wbd::Float64,
                         wbr::Float64,
                         vrx::Array{Float64},
                         vry::Array{Float64})::Bool

    for (ix, iy, iyaw) in zip(x, y, yaw)
        cx = ix + wbd*cos(iyaw)
        cy = iy + wbd*sin(iyaw)

        # Whole bubble check
        ids = inrange(kdtree, [cx, cy], wbr, true)
        # println(length(ids))
        if length(ids) == 0 continue end

        if !rect_check(ix, iy, iyaw, ox[ids], oy[ids], vrx, vry)
            # println("collision")
            return false #collision
        end

        # println(ids)
    end

    return true #OK

end

function rect_check(ix::Float64, iy::Float64, iyaw::Float64,
                    ox::Array{Float64}, oy::Array{Float64},
                    vrx::Array{Float64}, vry::Array{Float64}
                   )

    c = cos(-iyaw)
    s = sin(-iyaw)

    for (iox, ioy) in zip(ox, oy)
        tx = iox - ix
        ty = ioy - iy
        lx = (c*tx - s*ty)
        ly = (s*tx + c*ty)

        sumangle = 0.0
        for i in 1:length(vrx)-1
            x1 = vrx[i] - lx
            y1 = vry[i] - ly
            x2 = vrx[i+1] - lx
            y2 = vry[i+1] - ly
            d1 = hypot(x1,y1)
            d2 = hypot(x2,y2)
            theta1 = atan(y1,x1)
            tty = (-sin(theta1)*x2 + cos(theta1)*y2)
            tmp = (x1*x2+y1*y2)/(d1*d2)

            if tmp >= 1.0
                tmp = 1.0
            elseif tmp <= 0.0
                tmp = 0.0
            end

            if tty >= 0.0
                sumangle += acos(tmp)
            else
                sumangle -= acos(tmp)
            end
        end

        if sumangle >= pi
            return false #collision
        end
    end

    return true #OK
end

function calc_trailer_yaw_from_xyyaw(
                   x::Array{Float64},
                   y::Array{Float64},
                   yaw::Array{Float64},
                   init_tyaw::Float64,
                   steps::Array{Float64})::Array{Float64}
    """
    calc trailer yaw from x y yaw lists
    """

    tyaw = fill(0.0, length(x))
    tyaw[1] = init_tyaw

    ## for trailer truck
    # for i in 2:length(x)
    #     tyaw[i] += tyaw[i-1] + steps[i-1]/LT*sin(yaw[i-1] - tyaw[i-1])
    # end

    for i in 2:length(x)
        tyaw[i] = yaw[i]# tyaw[i-1] + steps[i-1]/LT*sin(yaw[i-1] - tyaw[i-1])
    end
    return tyaw
end

function trailer_motion_model(x, y, yaw0, yaw1, D, d, L, delta)
    """
    Motion model for trailer
    see:
    http://planning.cs.uiuc.edu/node661.html#77556
    """
    x += D*cos(yaw0)
    y += D*sin(yaw0)
    yaw0 += D/L*tan(delta)
    yaw1 += D/d*sin(yaw0 - yaw1)

    return x, y, yaw0, yaw1
end

function check_trailer_collision(
                    ox::Array{Float64},
                    oy::Array{Float64},
                    x::Array{Float64},
                    y::Array{Float64},
                    yaw0::Array{Float64},
                    yaw1::Array{Float64};
                    kdtree = nothing
                   )
    """
    collision check function for trailer

    """

    if kdtree == nothing
		kdtree = KDTree([ox'; oy'])
    end

    vrxt = [LTF, LTF, -LTB, -LTB, LTF]
    vryt = [-W/2.0, W/2.0, W/2.0, -W/2.0, -W/2.0]

    # bubble parameter
    DT = (LTF + LTB)/2.0 - LTB
    DTR = (LTF + LTB)/2.0 + 0.3

    # check trailer
    if !check_collision(x, y, yaw1, kdtree, ox, oy, DT, DTR, vrxt, vryt)
        return false
    end

    vrxf = [LF, LF, -LB, -LB, LF]
    vryf = [-W/2.0, W/2.0, W/2.0, -W/2.0, -W/2.0]

    # bubble parameter
    DF = (LF + LB)/2.0 - LB
    DFR = (LF + LB)/2.0 + 0.3

    # check front trailer
    if !check_collision(x, y, yaw0, kdtree, ox, oy, DF, DFR, vrxf, vryf)
        return false
    end

    return true #OK
end

function plot_trailer(x::Float64,
                      y::Float64,
                      yaw::Float64,
                      yaw1::Float64,
                      steer::Float64)

    truckcolor = "-k"

    LENGTH = LB+LF
    LENGTHt = LTB+LTF

    truckOutLine = [-LB (LENGTH - LB) (LENGTH - LB) (-LB) (-LB);
                    W/2 W/2 -W/2 -W/2 W / 2]
    trailerOutLine = [-LTB (LENGTHt - LTB) (LENGTHt - LTB) (-LTB) (-LTB);
                    W/2 W/2 -W/2 -W/2 W / 2]

    rr_wheel = [TR -TR -TR TR TR;
                -W/12.0+TW  -W/12.0+TW W/12.0+TW W/12.0+TW -W/12.0+TW]

    rl_wheel = [TR -TR -TR TR TR;
                -W/12.0-TW  -W/12.0-TW W/12.0-TW W/12.0-TW -W/12.0-TW]

    fr_wheel = [TR -TR -TR TR TR
                -W/12.0+TW  -W/12.0+TW W/12.0+TW W/12.0+TW -W/12.0+TW]

    fl_wheel = [TR -TR -TR TR TR;
                -W/12.0-TW  -W/12.0-TW W/12.0-TW W/12.0-TW -W/12.0-TW]
    tr_wheel = [TR -TR -TR TR TR
                -W/12.0+TW  -W/12.0+TW W/12.0+TW W/12.0+TW -W/12.0+TW]

    tl_wheel = [TR -TR -TR TR TR;
                -W/12.0-TW  -W/12.0-TW W/12.0-TW W/12.0-TW -W/12.0-TW]

    Rot1 = [cos(yaw) sin(yaw);
           -sin(yaw) cos(yaw)]
    Rot2 = [cos(steer) sin(steer);
           -sin(steer) cos(steer)]
    Rot3 = [cos(yaw1) sin(yaw1);
           -sin(yaw1) cos(yaw1)]

    fr_wheel = (fr_wheel' * Rot2)'
    fl_wheel = (fl_wheel' * Rot2)'
    fr_wheel[1,:] .+= WB
    fl_wheel[1,:] .+= WB
    fr_wheel = (fr_wheel' * Rot1)'
    fl_wheel = (fl_wheel' * Rot1)'

    tr_wheel[1,:] .-= LT
    tl_wheel[1,:] .-= LT
    tr_wheel = (tr_wheel' * Rot3)'
    tl_wheel = (tl_wheel' * Rot3)'

    truckOutLine = (truckOutLine' * Rot1)'
    trailerOutLine = (trailerOutLine' * Rot3)'
    rr_wheel = (rr_wheel' * Rot1)'
    rl_wheel = (rl_wheel' * Rot1)'

    truckOutLine[1, :] .+= x
    truckOutLine[2, :] .+= y
    trailerOutLine[1, :] .+= x
    trailerOutLine[2, :] .+= y
    fr_wheel[1, :] .+= x
    fr_wheel[2, :] .+= y
    rr_wheel[1, :] .+= x
    rr_wheel[2, :] .+= y
    fl_wheel[1, :] .+= x
    fl_wheel[2, :] .+= y
    rl_wheel[1, :] .+= x
    rl_wheel[2, :] .+= y

    tr_wheel[1, :] .+= x
    tr_wheel[2, :] .+= y
    tl_wheel[1, :] .+= x
    tl_wheel[2, :] .+= y

    plot(truckOutLine[1, :], truckOutLine[2, :], truckcolor)
    plot(trailerOutLine[1, :], trailerOutLine[2, :], truckcolor)
    plot(fr_wheel[1, :], fr_wheel[2, :], truckcolor)
    plot(rr_wheel[1, :], rr_wheel[2, :], truckcolor)
    plot(fl_wheel[1, :], fl_wheel[2, :], truckcolor)
    plot(rl_wheel[1, :], rl_wheel[2, :], truckcolor)

    plot(tr_wheel[1, :], tr_wheel[2, :], truckcolor)
    plot(tl_wheel[1, :], tl_wheel[2, :], truckcolor)
    plot(x, y, "*")
end

function main()
    x = 0.0
    y = 0.0
    yaw0 = deg2rad(-90.0)
    yaw1 = deg2rad(-90.0)

    plot_trailer(x, y, yaw0, yaw1, 0.0)

    DF = (LF + LB)/2.0 - LB
    DFR = (LF + LB)/2.0 + 0.3

    DT = (LTF + LTB)/2.0 - LTB
    DTR = (LTF + LTB)/2.0 + 0.3

    axis("equal")

    show()

end

# @time main()
if length(PROGRAM_FILE)!=0 &&
	occursin(PROGRAM_FILE, @__FILE__)
    @time main()
end

end

max_stars_repo_name: sunwenjun/HybridAStarTrailer

max_issues_repo_name: sunwenjun/HybridAStarTrailer

max_forks_repo_name: sunwenjun/HybridAStarTrailer


text:
mutable struct Data{T}
    """
    Data structure for the solution of
        min         ½x'Px + q'x
        subject to  Ax ≤ b
    with an active set algorithm.

    The most important element is F, which holds
    - F.QR:      an updatable QR factorization of the working constraints
    - F.Z:       a view on an orthonormal matrix spanning the nullspace of the working constraints
    - F.P:       the hessian of the problem
    - F.U & F.D: the ldl factorization of the projected hessian, i.e. F.U'*F.D*F.D = F.Z'*F.P*F.Z

    The Data structure also keeps matrices of the constraints not in the working set (A_ignored and b_ignored)
    stored in a continuous manner. This is done to increase use of BLAS.
    """
    x::Vector{T}

    n::Int
    m::Int

    F::NullspaceHessianLDL{T}
    q::Vector{T}
    A::Matrix{T}
    b::Vector{T}
    working_set::Vector{Int}
    ignored_set::Vector{Int}
    λ::Vector{T}
    residual::T

    iteration::Int
    done::Bool

    e::Vector{T} # Just an auxiliary vector used in computing step directions

    A_ignored::SubArray{T, 2, Matrix{T}, Tuple{UnitRange{Int}, Base.Slice{Base.OneTo{Int}}}, false}
    b_ignored::SubArray{T, 1, Vector{T}, Tuple{UnitRange{Int}}, true}
    A_shuffled::Matrix{T}
    b_shuffled::Vector{T}

    verbosity::Int
    printing_interval::Int
    r_min::T
    r_max::T

    function Data(P::Matrix{T}, q::Vector{T}, A::Matrix{T}, b::Vector{T},
        x::Vector{T}; r_min::T=zero(T), r_max::T=Inf, verbosity=1, printing_interval=50) where T
        if r_min < 1e-9; r_min = -one(T); end

        m, n = size(A)
        working_set = zeros(Int, 0)
        @assert(maximum(A*x - b) < 1e-9, "The initial point is infeasible!")
        @assert(norm(x) - r_min > -1e-9, "The initial point is infeasible!")
        @assert(norm(x) - r_max < 1e-9, "The initial point is infeasible!")

        ignored_set = setdiff(1:m, working_set)

        F = NullspaceHessianLDL(P, Matrix(view(A, working_set, :)'))
        if F.m == 0 # To many artificial constraints...
            remove_constraint!(F, 0)
        end
        A_shuffled = zeros(T, m, n)
        l = length(ignored_set)
        A_shuffled[end-l+1:end, :] .= view(A, ignored_set, :)
        b_shuffled = zeros(T, m)
        b_shuffled[end-l+1:end] .= view(b, ignored_set)

        e = zeros(T, n);
        λ = zeros(T, m);

        new{T}(x, n, m, F, q, A, b, working_set, ignored_set, λ,
            NaN, 0, false, e,
            view(A_shuffled, m-l+1:m, :),
            view(b_shuffled, m-l+1:m),
            A_shuffled, b_shuffled,
            verbosity, printing_interval,
            r_min, r_max)
    end
end

function solve(P::Matrix{T}, q::Vector{T}, A::Matrix{T}, b::Vector{T},
    x::Vector{T}; max_iter::Int=5000, kwargs...) where T

    data = Data(P, q, A, b, x; kwargs...)

    if data.verbosity > 0
        print_header(data)
        print_info(data)
    end

    while !data.done && data.iteration <= max_iter && norm(data.x) <= data.r_max - 1e-10 && norm(data.x) >= data.r_min + 1e-10
        iterate!(data)

        if data.verbosity > 0
            mod(data.iteration, 10*data.printing_interval) == 0 && print_header(data)
            mod(data.iteration, data.printing_interval) == 0 && print_info(data)
        end
    end
    data.verbosity > 0 && print_info(data)

    return x
end

function iterate!(data::Data{T}) where{T}
    direction, stepsize, new_constraints = calculate_step(data)
    data.x .+= stepsize*direction

    if !isempty(new_constraints)
        add_constraint!(data, new_constraints[1])
    end
    #ToDo: break next condition into simpler ones or a more representative flag
    if (isempty(new_constraints) || data.F.m == 0) && norm(data.x) <= data.r_max - 1e-10
        if data.F.artificial_constraints > 0
            remove_constraint!(data.F, 0)
        else
            idx = check_kkt!(data)
            !data.done && remove_constraint!(data, idx)
        end
    end
    data.iteration += 1
end

function calculate_step(data)
    gradient = data.F.P*data.x + data.q
    if data.F.D[end] >= data.F.indefinite_tolerance
        qw = data.F.Z'*(gradient)
        if norm(qw) <= 1e-10 # We are alread on the optimizer of the current subproblem
            return zeros(data.n), 0, []
        end
        # Gill & Murray (1978) calculate the direction/stepsize as:
        # direction = data.F.Z*reverse(data.F.U\e)
        # direction .*= -sign(dot(direction, gradient))
        # α_min = -dot(gradient, direction)/data.F.D[end]
        # This assumes that we start from a feasible vertex or a feasible stationary point.
        # This assumption is suitable e.g. when we use a Phase I algorithm to find a feasible point
        # and it allows for a faster/simpler calculation of the direction.
        # Nevertheless, for purposes of generality, we calculate the direction in a more general way.
        direction = -data.F.Z*reverse(data.F.U\(data.F.D\(data.F.U'\reverse(qw))))
        α_min = 1
    else
        e = view(data.e, 1:data.F.m); e[end] = 1
        #=
        if norm(data.F.U[:, end]) <= 1e-11
            data.F.U[end] = 1e-11
        end
        =#
        direction = data.F.Z*reverse(data.F.U\e)
        if dot(direction, gradient) >= 0
            direction .= -direction
        end
        e .= 0
        α_min = Inf
    end

    Α_times_direction = data.A_ignored*direction
    ratios = abs.(data.b_ignored - data.A_ignored*data.x)./(Α_times_direction)
    ratios[Α_times_direction .<= 1e-11] .= Inf

    α_constraint = minimum(ratios)
    idx = findlast(ratios .== α_constraint)
    # @show α_constraint, data.ignored_set[idx]

    α = min(α_min, α_constraint)
    if α == Inf
        # variable "direction" holds the unbounded ray
        @info "GeneralQP.jl has detected the problem to be unbounded (unbounded ray found)."
    end

    α_max = Inf
    if isfinite(data.r_max)
        # Calculate the maximum allowable step α_max so that r_min ≤ norm(x) ≤ r_max
        roots_rmax = roots(Poly([norm(data.x)^2 - data.r_max^2, 2*dot(direction, data.x), norm(direction)^2]))
        if data.r_min > 0
            roots_rmin = roots(Poly([norm(data.x)^2 - data.r_min^2, 2*dot(direction, data.x), norm(direction)^2]))
            roots_all = [roots_rmin; roots_rmax]
        else
            roots_all = roots_rmax
        end
        # Discard complex and negative steps
        roots_all = real.(roots_all[isreal.(roots_all)])
        roots_all = roots_all[roots_all .>= 0]
        if length(roots_all) > 0
            α_max = minimum(roots_all)
        end
    end
    stepsize = min(α, α_max)
    @assert(isfinite(stepsize), "Set the keyword argument r_max to a finite value if you want to get bounded solutions.")
    if α_constraint == stepsize
        new_constraints = [idx]
    else
        new_constraints = []
    end

    return direction, stepsize, new_constraints
end

function check_kkt!(data)
    grad = data.F.P*data.x + data.q
    λ = -data.F.QR.R1\data.F.QR.Q1'*grad
    data.λ .= 0.0
    data.λ[data.working_set] .= λ
    data.residual = norm(data.F.Z'*grad)
    # data.residual = norm(grad + data.A[data.working_set, :]'*λ)

    idx = NaN
    if all(λ .>= -1e-8)
        data.done = true
    else
        data.done = false
        idx = argmin(λ)
    end

    return idx
end

ext: jl


text:
@rlibrary ggforce

"""
    Plot landmark positions of guided of forward simulated path
    X: sample path containing landmark evolution
    Xᵒ: guided sample path containing landmark evolution
    n: number of landmarks
    model: either :ms or :ahs
    v0: starting configuration (positions)
    vT: ending configuration (positions)
    nfs: contains info over noisefields
    db: domainbound (plot on [-db,db] x [-db,db])
"""
function plotlandmarkpositions(X,Xᵒ,P,v0,vT;db=5)
    n = P.n
    if isa(P,Landmarks)
        nfs = P.nfs
    elseif isa(P,MarslandShardlow)
        nfs = 0.0
    end
    # construct df for X
    out = [Any[X.tt[i], [X.yy[i][CartesianIndex(c, k)][l] for l in 1:d, c in 1:2]..., "point$k"] for k in 1:n, i in eachindex(X.tt) ][:]
    df = DataFrame(time=extractcomp(out,1),pos1=extractcomp(out,2),pos2=extractcomp(out,3),mom1=extractcomp(out,4),mom2=extractcomp(out,5),pointID=extractcomp(out,6))
    # construct df for Xᵒ
    outg = [Any[Xᵒ.tt[i], [Xᵒ.yy[i][CartesianIndex(c, k)][l] for l in 1:d, c in 1:2]..., "point$k"] for k in 1:n, i in eachindex(Xᵒ.tt) ][:]
    dfg = DataFrame(time=extractcomp(outg,1),pos1=extractcomp(outg,2),pos2=extractcomp(outg,3),mom1=extractcomp(outg,4),mom2=extractcomp(outg,5),pointID=extractcomp(outg,6))
    # construct df for noisefields
    if isa(P,Landmarks)
        nfscales = [nfs[j].γ for j in eachindex(nfs)]
        nfsloc = [nfs[j].δ for j in eachindex(nfs)]
        df_nfs = DataFrame(locx =  extractcomp(nfsloc,1), locy =  extractcomp(nfsloc,2),
        lambda1=  extractcomp(nfscales,1), lambda2=extractcomp(nfscales,2), nfstd=fill(nfstd,length(nfs)))
    elseif isa(P,MarslandShardlow)
        df_nfs =DataFrame(locx=Int64[], locy=Int64[], lambda1=Int64[], lambda2=Int64[],nfstd=Int64[])
    end

    # construct df for initial and final configs
    # (need to repeat to connect also beginning and end)
    dfT = DataFrame(pos1=repeat(extractcomp(vT,1),2), pos2=repeat(extractcomp(vT,2),2))
    df0= DataFrame(pos1=repeat(extractcomp(v0,1),2), pos2=repeat(extractcomp(v0,2),2))

    if isa(P,MarslandShardlow)
        titel = "Marsland-Shardlow model, "
    elseif isa(P,Landmarks)
        titel = "Arnaudon-Holm-Sommer model, "
    end
    titel = titel * string(n)*" landmarks"

    @rput titel
    @rput df
    @rput dfg
    @rput df_nfs
    @rput dfT
    @rput df0
    @rput db

    R"""
    library(tidyverse)
    library(ggforce) # for circles
    dfpaths <- rbind(df,dfg)
    dfpaths$fg <- rep(c("forward","guided"),each=nrow(df))

    g <- ggplot() +
    geom_segment(data=df_nfs,aes(x=locx, y=locy, xend=locx, yend=locy+lambda2), arrow = arrow(length = unit(.1, "cm")),color="Grey")+
    geom_segment(data=df_nfs,aes(x=locx, y=locy, xend=locx+lambda1, yend=locy), arrow = arrow(length = unit(.1, "cm")),color="Grey")+
     geom_circle(aes(x0 = locx, y0 = locy, r = nfstd), data = df_nfs,color="Grey",linetype="dashed") +
     geom_path(data=df0, mapping=aes(x=pos1,y=pos2),colour='black',size=0.7) +
     geom_path(data=dfT, mapping=aes(x=pos1,y=pos2),colour='orange',size=0.7) +
     geom_point(data=dfT, mapping=aes(x=pos1,y=pos2),colour='orange',size=0.7) +
     geom_point(data=df0, mapping=aes(x=pos1,y=pos2),colour='black',size=0.7)+
     geom_path(data=dfpaths, mapping=aes(x=pos1,pos2,group=pointID,colour=time)) +
          facet_wrap(~fg) + geom_point(data=df_nfs, aes(x=locx,y=locy),size=2,color="Grey")+
              theme_light() + xlab('horizontal position') +
               coord_cartesian(xlim = c(-db,db), ylim = c(-db,db)) +
          ylab('vertical position') + ggtitle(titel)
    show(g)
    """
end

"""
    Plot landmark positions of both guided and forward simulated paths
    X: sample path containing landmark evolution
    n: number of landmarks
    model: either :ms or :ahs
    v0: starting configuration (positions)
    vT: ending configuration (positions)
    nfs: contains info over noisefields
    db: domainbound (plot on [-db,db] x [-db,db])
"""
function plotlandmarkpositions(X,P,v0,vT;db=5)
    n = P.n
    if isa(P,Landmarks)
        nfs = P.nfs
    elseif isa(P,MarslandShardlow)
        nfs = 0.0
    end
    # construct df for X
    out = [Any[X.tt[i], [X.yy[i][CartesianIndex(c, k)][l] for l in 1:d, c in 1:2]..., "point$k"] for k in 1:n, i in eachindex(X.tt) ][:]
    df = DataFrame(time=extractcomp(out,1),pos1=extractcomp(out,2),pos2=extractcomp(out,3),mom1=extractcomp(out,4),mom2=extractcomp(out,5),pointID=extractcomp(out,6))
    # construct df for noisefields
    if isa(P,Landmarks)
        nfscales = [nfs[j].γ for j in eachindex(nfs)]
        nflocs = [nfs[j].δ for j in eachindex(nfs)]
        df_nfs = DataFrame(locx =  extractcomp(nflocs,1), locy =  extractcomp(nflocs,2),
        lambda1=  extractcomp(nfscales,1), lambda2=extractcomp(nfscales,2), nfstd=fill(nfstd,length(nfs)))
    elseif isa(P,MarslandShardlow)
        df_nfs =DataFrame(locx=Int64[], locy=Int64[], lambda1=Int64[], lambda2=Int64[],nfstd=Int64[])
    end

    # construct df for initial and final configs
    # (need to repeat to connect also beginning and end)
    dfT = DataFrame(pos1=repeat(extractcomp(vT,1),2), pos2=repeat(extractcomp(vT,2),2))
    df0= DataFrame(pos1=repeat(extractcomp(v0,1),2), pos2=repeat(extractcomp(v0,2),2))

    if isa(P,MarslandShardlow)
        titel = "Marsland-Shardlow model, "
    elseif isa(P,Landmarks)
        titel = "Arnaudon-Holm-Sommer model, "
    end
    titel = titel * string(n)*" landmarks"

    @rput titel
    @rput df

    @rput df_nfs
    @rput dfT
    @rput df0
    @rput db

    R"""
    library(tidyverse)
    library(ggforce) # for circles
    dfpaths <- df

    g <- ggplot() +
    geom_segment(data=df_nfs,aes(x=locx, y=locy, xend=locx, yend=locy+lambda2), arrow = arrow(length = unit(.1, "cm")),color="Grey")+
    geom_segment(data=df_nfs,aes(x=locx, y=locy, xend=locx+lambda1, yend=locy), arrow = arrow(length = unit(.1, "cm")),color="Grey")+
     geom_circle(aes(x0 = locx, y0 = locy, r = nfstd), data = df_nfs,color="Grey",linetype="dashed") +
     geom_path(data=df0, mapping=aes(x=pos1,y=pos2),colour='black',size=0.7) +
     geom_path(data=dfT, mapping=aes(x=pos1,y=pos2),colour='orange',size=0.7) +
     geom_point(data=dfT, mapping=aes(x=pos1,y=pos2),colour='orange',size=0.7) +
     geom_point(data=df0, mapping=aes(x=pos1,y=pos2),colour='black',size=0.7)+
     geom_path(data=dfpaths, mapping=aes(x=pos1,pos2,group=pointID,colour=time)) +
              theme_light() + xlab('horizontal position') +
              coord_cartesian(xlim = c(-db,db), ylim = c(-db,db))+
          ylab('vertical position') + ggtitle(titel)
    show(g)
    """
end

"""
    On the square grid provided by (-db:δ:db) x (-db:δ:db) compute
    σq(q) and make heatmap
"""
function plotσq(db,nfs ;δ=0.1)
    rfine = -db:δ:db
    nfloc_fineres = Point.(collect(product(rfine, rfine)))[:]
    σq_fineres = σq(nfs).(nfloc_fineres)
    df_sigmaq = DataFrame(x=extractcomp(nfloc_fineres,1), y=extractcomp(nfloc_fineres,2), horiz =map(x->x[1,1], σq_fineres), vertic = map(x->x[2,2], σq_fineres) )
    if model ==:ahs
        nfscales = [nfs[j].γ for j in eachindex(nfs)]
        nfsloc = [nfs[j].δ for j in eachindex(nfs)]
        df_nfs = DataFrame(locx =  extractcomp(nfsloc,1), locy =  extractcomp(nfsloc,2),
        lambda1=  extractcomp(nfscales,1), lambda2=extractcomp(nfscales,2), nfstd=fill(nfstd,length(nfs)))
    else
        df_nfs =DataFrame(locx=Int64[], locy=Int64[], lambda1=Int64[], lambda2=Int64[],nfstd=Int64[])
    end
    @rput df_sigmaq
    @rput df_nfs
    R"""
    library(tidyverse)
    library(gridExtra)
    phor <- ggplot() + geom_point(data=df_sigmaq, aes(x=x, y=y,colour=horiz)) +
        ggtitle("Horizontal noise field on position")#+ geom_point(data=df_nfs, aes(x=locx,y=locy),size=2,color="Grey")
    pvert <- ggplot() + geom_point(data=df_sigmaq, aes(x=x, y=y,colour=vertic)) +
        ggtitle("Vertical noise field on position")#+ geom_point(data=df_nfs, aes(x=locx,y=locy),size=2,color="Grey")
    grid.arrange(phor,pvert)
    """
end

# R"""
# dfg %>% dplyr::select(time,mom1,mom2,pointID) %>%
#     filter(pointID %in% c("point1","point2","point3","point40","point41","point42")) %>%
#      ggplot(aes(x=mom1,mom2,group=pointID,colour=time)) + geom_path() +
#      facet_wrap(~pointID,scales="free") + theme_light()
# """

if false
    R"""
    dfg %>% dplyr::select(time,mom1,mom2,pointID) %>%
    filter(pointID %in% c("point1","point2","point3","point40","point41","point42")) %>%
     ggplot(aes(x=mom1,mom2,group=pointID,colour=time)) + geom_path() +
     facet_wrap(~pointID) + theme_light()
    """
    R"""
    library(tidyverse)
    dfsub <- df %>% dplyr::select(time,pos1,pos2,pointID) %>% filter(time<0.99)
    dfsubg <- dfg %>% dplyr::select(time,pos1,pos2,pointID) %>% filter(time<0.99)

    sub <- rbind(dfsub,dfsubg)
    sub$fg <- rep(c("forward","guided"),each=nrow(dfsub))

    ggplot() +
          geom_path(data=sub, mapping=aes(x=pos1,pos2,group=pointID,colour=fg)) +
          geom_point(data=dfT, mapping=aes(x=pos1,y=pos2),colour='orange',size=0.7) +
          geom_point(data=df0, mapping=aes(x=pos1,y=pos2),colour='black',size=0.7)+
          #facet_wrap(~fg) +
          theme_light() + xlab('horizontal position') +
          ylab('vertical position') + ggtitle(titel)
    """
end
if false
    R"""
    df %>% dplyr::select(time,mom1,mom2,pointID) %>%
        filter(pointID %in% c("point1","point2","point3","point40","point41","point42")) %>%
         ggplot(aes(x=mom1,mom2,group=pointID,colour=time)) + geom_path() +
         facet_wrap(~pointID,scales="free") + theme_light()
    """
end


text:
module Junet

# Use-only imports
import Random: randexp  # for generate/random.jl
import GZip             # for io/edgelist.jl
using Cairo             # for draw/plot.jl
import Colors: RGB      # for draw/plot.jl

# Imports to extend Base
import Base:
    length, size, getindex, setindex!,  # common
    @_inline_meta, @_propagate_inbounds_meta,
    isless, ==,             # primitives.jl
    reverse, adjoint,       # graph.jl
    iterate, lastindex, eltype,  # graph_iter.jl
    IteratorEltype, HasEltype, IteratorSize, HasLength,
    deleteat!,              # attributes.jl
    push!, pop!, rand,      # util/rand_dist.jl
    Matrix,                 # io/matrix.jl
    summary, show           # io/show.jl

include("params.jl")                        # getting to a minimally usable state
include("primitives.jl")
include("graph.jl")
include("attributes.jl")
include("graph_info.jl")
include("graph_iter.jl")
include("graph_operations.jl")

include("util/discrete_sampler.jl")         # commonly used methods
include("util/relabel.jl")

include("generate/classic.jl")              # graph generators
include("generate/random.jl")
include("generate/small_world.jl")
include("generate/preferential.jl")

include("measure/centrality/degree.jl")     # graph measures
include("measure/centrality/subtree.jl")
include("measure/centrality/power_iter.jl")
include("measure/centrality/closeness.jl")
include("measure/dyadic.jl")
include("measure/triadic/clustering.jl")
include("measure/global/components.jl")
include("measure/global/k-core.jl")
include("community/label_propagation.jl")   # community detection

include("io/show.jl")                       # input/output
include("io/matrix.jl")
include("io/formats/edgelist.jl")
include("io/formats/pajek.jl")
include("io/formats/gml.jl")
include("io/formats/ngraph.jl")
include("io/formats/ucinet.jl")
include("io/file_io.jl")

include("draw/layout.jl")                   # visualization
include("draw/node_paths.jl")
include("draw/edge_paths.jl")
include("draw/styles.jl")
include("draw/plot.jl")

export
    # graph.jl
    Graph, DirectedGraph, UndirectedGraph, MultiGraph, SimpleGraph, LightGraph,
    isdirected, ismultigraph, directed, undirected,

    # graph_info.jl
    memsize,
    nodes, nodecount, edgecount, maxedgecount,

    # graph_operations.jl
    neighbors, outneighbors, inneighbors,
    addnode!, addnodes!, remnode!, remnodes!,
    edges, outedges, inedges,
    addedge!, addedges!, hasedge, remedge!, remedges!,
    getindex, setindex!,
    addnodeattr!, addedgeattr!, addnodeattrs!, addedgeattrs!,

    # generate/
    graph_path, graph_cycle, graph_star, graph_wheel, graph_complete,
    graph_grid, graph_web, graph_tree,
    graph_gilbert, graph_erdos_renyi, graph_random, graph_erdosrenyi, erdosrenyi,
    graph_small_world, graph_smallworld, smallworld,
    graph_preferential,

    # measure/centrality/
    degree, indegree, outdegree,        # degree.jl
    isolates, isolatecount,
    leaves, leafcount,
    sources, sourcecount,
    sinks, sinkcount,

    subtree_size,                       # subtree.jl

    eigenvector, pagerank,              # power_iter.jl

    closeness, harmonic_centrality,     # closeness.jl

    # measure/dyadic.jl
    density, mutuality, reciprocity, dyadcensus,
    selfloopcount, selfloopnodes,

    # measure/triadic/
    maxtrianglecount, triangles, clustering,

    # measure/global/
    components, kcores,

    # io/
    readgraph,

    # community/
    label_propagation,

    # draw/
    layout_random, layout_circle, layout_line,
    layout_fruchterman_reingold,
    plot

end


text:
function solve(r::Rheology,p::Params,σᵢⱼ_i,ϵᵢⱼ_i,D_i,ϵ̇11,Δt)
    ps = p.solver
    D = D_i
    # get first guess of the unknowns with an elastic solve
    ϵᵢⱼnext = insert_into(ϵᵢⱼ_i, (ϵᵢⱼ_i[1,1] + ϵ̇11*Δt), (1,1))
    σᵢⱼnext = σᵢⱼ_i#σᵢⱼnext = compute_σij(r,D,ϵᵢⱼnext)
    #σᵢⱼnext = insert_into(σᵢⱼnext, -1e6, (2,2))
    u = Vec(σᵢⱼnext[1,1], σᵢⱼnext[3,3], ϵᵢⱼnext[2,2])
    #@debug "u_i = $u"
    for i in 1:ps.newton_maxiter
        # get residual and its gradient with respect to u
        ∇res , res = Tensors.gradient(u -> residual(r,D,ϵᵢⱼ_i,ϵᵢⱼnext,σᵢⱼ_i,σᵢⱼnext,Δt,u), u, :all)
        #@debug "norm res = $(norm(res))"
        #@debug "∇res = $∇res"

        # update u with Newton algo
        δu = - ∇res\res
        u = u + δu
        #@debug "δu = $δu"
        #@debug "typeof(u) = $(typeof(u))"

        (norm(res) <= ps.newton_abstol) && break

        (i == ps.newton_maxiter) && @debug("ending norm res = $(norm(res))")#("max newton iteration reached ($i), residual still higher than abstol with $(norm(res))")
    end
    # update ϵᵢⱼnext and σᵢⱼnext with converged u
    ϵᵢⱼnext = insert_into(ϵᵢⱼnext, u[3], (2,2))
    σᵢⱼnext = insert_into(σᵢⱼnext, (u[1], u[2]), ((1,1),(3,3)))
    _ , Ḋ = compute_ϵ̇ij(r,D,σᵢⱼ_i,σᵢⱼnext,Δt)
    D = D + Ḋ*Δt
    return σᵢⱼnext, ϵᵢⱼnext, D, u
    end

    function residual(r,D,ϵij,ϵijnext,σij,σijnext,Δt,u)
    σ11next, σ33next, ϵ22next = u
    σijnext = insert_into(σijnext, (σ11next, σ33next), ((1,1),(3,3)))
    ϵijnext = insert_into(ϵijnext, ϵ22next, (2,2))
    ϵ̇ij_analytical, _ = compute_ϵ̇ij(r,D,σij,σijnext,Δt)
    ϵ̇ij = (ϵijnext - ϵij)/Δt
    Δϵ̇ij = ϵ̇ij_analytical - ϵ̇ij
    #@debug " KI = $(compute_KI(r,σijnext, D))"
    return Vec(Δϵ̇ij[1,1],Δϵ̇ij[2,2],Δϵ̇ij[3,3])
end

function adaptative_solve_1_point(r,p,σᵢⱼ_i,ϵᵢⱼ_i,D_i,ϵ̇11,θ,Δt)
    e₀ = p.solver.e₀

    σᵢⱼnext1, ϵᵢⱼnext1, Dnext1, u1 = solve_1_point(r,p,σᵢⱼ_i,ϵᵢⱼ_i,D_i,ϵ̇11,θ,Δt)
    σᵢⱼmid, ϵᵢⱼmid, Dmid, umid = solve_1_point(r,p,σᵢⱼ_i,ϵᵢⱼ_i,D_i,ϵ̇11,θ,Δt/2)
    σᵢⱼnext2, ϵᵢⱼnext2, Dnext2, u2 = solve_1_point(r,p,σᵢⱼmid,ϵᵢⱼmid,Dmid,ϵ̇11,θ,Δt/2)

    # compute errors for each unknowns physical quantity and ponderate
    eD = Dnext2-Dnext1
    eσ = norm((u2-u1)[1:2])
    eϵ = norm((u2-u1)[3])
    if e₀ isa Real
        e₀ref = e₀
        e_normalized = (eD, eσ/r.G, eϵ)
        acceptable_error_flag = all(e_normalized.<e₀)
        e = maximum(e_normalized)
    elseif e₀ isa NamedTuple
        e₀ref = e₀.D
        acceptable_error_flag = (eD<e₀.D) && (eσ<(e₀.σ)) && (eϵ<e₀.ϵ)
        e_normalized = (eD, eσ*(e₀ref/e₀.σ), eϵ*(e₀ref/e₀.ϵ))
        e,ind = findmax(e_normalized)
        ok_flag || @debug("maximum error comes from indice $(ind) of (D,σ,ϵ)")
    end

    if acceptable_error_flag
        # increse timestep
        Δt_next = min(Δt*abs(e₀ref/e),Δt*2)
        # keep best solution
        return σᵢⱼnext2, ϵᵢⱼnext2, Dnext2, Δt, Δt_next
    else
        # recursively run with decreased timestep
        adaptative_solve_1_point(r,p,σᵢⱼ_i,ϵᵢⱼ_i,D_i,ϵ̇11,θ,Δt*abs(e₀ref/e))
        #initialy Δt*abs(e₀ref/e)^2 but without the square seems to generaly require less iterations.
    end
end

function solve_1_point(r::Rheology,p::Params,σᵢⱼ_i,ϵᵢⱼ_i,D_i,ϵ̇11,θ,Δt)
    # this function is similar to solve but takes and returns tensor in the shear band coordinates instead of the principal stresses coordinates
    ps = p.solver
    D = D_i

    ϵᵢⱼ_i = principal_coords(ϵᵢⱼ_i,θ)
    σᵢⱼ_i = principal_coords(σᵢⱼ_i,θ)
    #println(eps()*norm(σᵢⱼ_i))
    @assert isapprox(σᵢⱼ_i[1,2],0.0,atol=1e-12*norm(σᵢⱼ_i)) # put this in tests at some point

    # get first guess of the unknowns with an elastic solve
    ϵᵢⱼnext = insert_into(ϵᵢⱼ_i, (ϵᵢⱼ_i[1,1] + ϵ̇11*Δt), (1,1))
    σᵢⱼnext = σᵢⱼ_i#σᵢⱼnext = compute_σij(r,D,ϵᵢⱼnext)
    #σᵢⱼnext = insert_into(σᵢⱼnext, -1e6, (2,2))
    u = Vec(σᵢⱼnext[1,1], σᵢⱼnext[3,3], ϵᵢⱼnext[2,2])
    #@debug "u_i = $u"
    for i in 1:ps.newton_maxiter
        # get residual and its gradient with respect to u
        ∇res , res = Tensors.gradient(u -> residual(r,D,ϵᵢⱼ_i,ϵᵢⱼnext,σᵢⱼ_i,σᵢⱼnext,Δt,u), u, :all)
        #@debug "norm res = $(norm(res))"
        #@debug "∇res = $∇res"

        # update u with Newton algo
        δu = - ∇res\res
        u = u + δu
        #@debug "δu = $δu"
        #@debug "typeof(u) = $(typeof(u))"

        (norm(res) <= ps.newton_abstol) && break

        (i == ps.newton_maxiter) && @debug("ending norm res = $(norm(res))")#("max newton iteration reached ($i), residual still higher than abstol with $(norm(res))")
    end
    # update ϵᵢⱼnext and σᵢⱼnext with converged u
    ϵᵢⱼnext = insert_into(ϵᵢⱼnext, u[3], (2,2))
    σᵢⱼnext = insert_into(σᵢⱼnext, (u[1], u[2]), ((1,1),(3,3)))
    #_ , Ḋ = compute_ϵ̇ij(r,D,σᵢⱼ_i,σᵢⱼnext,Δt) # replaced by two next lines
    KI = compute_KI(r,(σᵢⱼnext+σᵢⱼ_i)/2,D)
    Ḋ = compute_subcrit_damage_rate(r,KI,D)
    Dnext = D + Ḋ*Δt
    return band_coords(σᵢⱼnext,θ), band_coords(ϵᵢⱼnext,θ), Dnext, u
end

function adaptative_solve_2_points(r,p,S,σ₃,σⁱᵢⱼ,σᵒᵢⱼ,ϵⁱᵢⱼ,Dᵒ,Dⁱ,ϵ̇ⁱξη,θ,irecursions,Δt ; damage_growth_out=true)
    # unpacking
    e₀ = p.solver.e₀
    flags = p.flags
    maxrecursions = p.solver.adaptative_maxrecursions
    irecursions += 1
    # solving procedures
    # integration over Δt
    Snext1, σⁱᵢⱼnext1, σᵒᵢⱼnext1, ϵⁱᵢⱼnext1, Dⁱnext1, Dᵒnext1, u1 = solve_2_points(r,p,S,σ₃,σⁱᵢⱼ,σᵒᵢⱼ,ϵⁱᵢⱼ,Dᵒ,Dⁱ,ϵ̇ⁱξη,θ,Δt ; damage_growth_out)
    flags.nan && (flags.nan1 = true)

    # integration over Δt/2 twice
    Snext2, σⁱᵢⱼnext2, σᵒᵢⱼnext2, ϵⁱᵢⱼnext2, Dⁱnext2, Dᵒnext2, u2 = Snext1, σⁱᵢⱼnext1, σᵒᵢⱼnext1, ϵⁱᵢⱼnext1, Dⁱnext1, Dᵒnext1, u1
    if !flags.nan
        Smid, σⁱᵢⱼmid, σᵒᵢⱼmid, ϵⁱᵢⱼmid, Dⁱmid, Dᵒmid, umid = solve_2_points(r,p,S,σ₃,σⁱᵢⱼ,σᵒᵢⱼ,ϵⁱᵢⱼ,Dᵒ,Dⁱ,ϵ̇ⁱξη,θ,Δt/2 ; damage_growth_out)
        flags.nan && (flags.nan2 = true)
        # run second part of the integration if first on did not produces NaNs
        Snext2, σⁱᵢⱼnext2, σᵒᵢⱼnext2, ϵⁱᵢⱼnext2, Dⁱnext2, Dᵒnext2, u2 = solve_2_points(r,p,Smid,σ₃,σⁱᵢⱼmid,σᵒᵢⱼmid,ϵⁱᵢⱼmid,Dᵒ,Dⁱmid,ϵ̇ⁱξη,θ,Δt/2 ; damage_growth_out)
        flags.nan && (flags.nan2 = true)
    end

    print_nans_error(p,Δt) # prints an error message and set nan_flag to true if nan1 or nan2 flags are true
    if flags.nan
        if irecursions == maxrecursions
            @warn "NaNs and maximum number of recursives calls to solve. Aborting the simulation"
            return S, σⁱᵢⱼ, σᵒᵢⱼ, ϵⁱᵢⱼ, Dⁱ, Dᵒ, Δt, Δt # return input values to save a last working state
        end
        factor = 0.2
        adaptative_solve_2_points(r,p,S,σ₃,σⁱᵢⱼ,σᵒᵢⱼ,ϵⁱᵢⱼ,Dᵒ,Dⁱ,ϵ̇ⁱξη,θ,irecursions, Δt*factor ; damage_growth_out)
    end

    # if nans
    # compute errors for each unknowns physical quantity and ponderate
    eD = max(abs((Dⁱnext2-Dⁱnext1)/Dⁱnext2),abs((Dᵒnext2-Dᵒnext1)/Dᵒnext2))
    eS = abs((Snext2-Snext1)/Snext2)
    eσⁱ = maximum(filter(!isnan,abs.((σⁱᵢⱼnext2 - σⁱᵢⱼnext1)./σⁱᵢⱼnext2)))
    eσᵒ = maximum(filter(!isnan,abs.((σᵒᵢⱼnext2 - σᵒᵢⱼnext1)./σᵒᵢⱼnext2)))
    eσ = max(eσⁱ,eσᵒ)
    eϵ = maximum(filter(!isnan,abs.((ϵⁱᵢⱼnext2 - ϵⁱᵢⱼnext1)./ϵⁱᵢⱼnext2)))
    e_vec = Vec(eD,eS,eσ,eϵ)
    e_norm_max = 0
    id_e_norm_max = 0
    if e₀ isa Real
        e_normalized = e_vec./e₀
        e_norm_max, id_e_norm_max = findmax(e_normalized)
    elseif e₀ isa NamedTuple
        e_normalized = e_vec ./ Vec(e₀.D,e₀.S,e₀.σ,e₀.ϵ)
        e_norm_max, id_e_norm_max = findmax(e_normalized)
    end
    flags.acceptable_error = (e_norm_max < 1)

    if irecursions == maxrecursions
        flags.acceptable_error = true
        components_tuple = ("D","S","σ","ϵ")
        @warn "50 recursive call to lower timestep, Δt = $(Δt), continuing with max relative error emax of $(e_norm_max) on $(components_tuple[id_e_norm_max])..."
    end

    if flags.acceptable_error
        # increse timestep
        # println("e₀ref = ", e₀ref)
        # println("emax = ", emax)
        Δt_next::Float64 = min(Δt/e_norm_max,Δt*2)
        # keep best solution
        return Snext2, σⁱᵢⱼnext2, σᵒᵢⱼnext2, ϵⁱᵢⱼnext2, Dⁱnext2, Dᵒnext2, Δt, Δt_next
    else
        # recursively run with decreased timestep
        # println("e₀ref = ", e₀ref)
        # println("emax = ", emax)
        factor::Float64 = 0.2
        adaptative_solve_2_points(r,p,S,σ₃,σⁱᵢⱼ,σᵒᵢⱼ,ϵⁱᵢⱼ,Dᵒ,Dⁱ,ϵ̇ⁱξη,θ,irecursions, Δt*factor ; damage_growth_out)
        #initialy Δt*abs(e₀ref/e)^2 but without the square seems to generaly require less iterations.
    end
end

function solve_2_points(r::Rheology,p::Params,S,σ₃,σⁱᵢⱼ,σᵒᵢⱼ,ϵⁱᵢⱼ,Dᵒ,Dⁱ,ϵ̇ⁱξη,θ,Δt ; damage_growth_out=true)
    # unpacking :
    ps = p.solver
    flags=p.flags

    # get first guess of the unknowns with an elastic solve
    u = Vec(S, σⁱᵢⱼ[1,1], σⁱᵢⱼ[3,3], ϵⁱᵢⱼ[2,2]) # Snext, σⁱξξnext, σⁱoopnext, ϵⁱηηnext
    #@debug "u_i = $u"
    for i in 1:ps.newton_maxiter
        # get residual and its gradient with respect to u
        #∇res , res = Tensors.gradient(u -> residual_2_points(r,S,σ₃,Dⁱ,Dᵒ,ϵⁱᵢⱼ,σⁱᵢⱼ,σᵒᵢⱼ,ϵ̇ⁱξη,Δt,u), u, :all)
        result = DiffResults.JacobianResult(u)
        try
            ForwardDiff.jacobian!(result, u -> residual_2_points(r,S,σ₃,Dⁱ,Dᵒ,ϵⁱᵢⱼ,σⁱᵢⱼ,σᵒᵢⱼ,ϵ̇ⁱξη,Δt,θ,u;damage_growth_out), u)
        catch
            @warn "Jacobian calculation generates nans.\n Exiting newton iter $(i) \n"
            flags.nan = true
            break
        end
        ∇res = DiffResults.jacobian(result)
        res  = DiffResults.value(result)

        # update u with Newton algo
        δu = - ∇res\res
        if any(isnan.(δu))
            @warn "Solution vector update δu contains nans : $(δu).\n Exiting newton iter $(i) with residual norm of $(norm(res)) \n"
            flags.nan = true
            break
        else
            u = u + δu
        end
        #@debug "δu = $δu"
        #@debug "typeof(u) = $(typeof(u))"

        (norm(res) <= ps.newton_abstol) && (@debug("Newton iter $i ending norm res = $(norm(res))") ; break)

        (i == ps.newton_maxiter) && @debug("Newton maxiter ending norm res = $(norm(res))")#("max newton iteration reached ($i), residual still higher than abstol with $(norm(res))")
    end
    # update strains and stresses with converged u
    # any(isnan.((S,Dⁱ,Dᵒ))) && println("(S,Dⁱ,Dᵒ) : ", (S,Dⁱ,Dᵒ))
    # any(isnan.(ϵⁱᵢⱼ)) && println("ϵⁱᵢⱼ : ",ϵⁱᵢⱼ)
    # any(isnan.(σⁱᵢⱼ)) && println("σⁱᵢⱼ : ",σⁱᵢⱼ)
    # any(isnan.(σᵒᵢⱼ)) && println("σᵒᵢⱼ : ",σᵒᵢⱼ)
    # any(isnan.(u)) && println("u : ",u)

    flags.nan && (return S, σⁱᵢⱼ, σᵒᵢⱼ, ϵⁱᵢⱼ, Dⁱ, Dᵒ, u)

    Ṡ, σ̇ᵒᵢⱼ, σ̇ⁱᵢⱼ, ϵ̇ⁱᵢⱼ = compute_stress_strain_derivatives_from_u(r,S,σ₃,Dⁱ,Dᵒ,ϵⁱᵢⱼ,σⁱᵢⱼ,σᵒᵢⱼ,ϵ̇ⁱξη,Δt,θ,u)
    Snext = S + Ṡ*Δt
    σᵒᵢⱼnext = σᵒᵢⱼ + σ̇ᵒᵢⱼ*Δt
    σⁱᵢⱼnext = σⁱᵢⱼ + σ̇ⁱᵢⱼ*Δt
    ϵⁱᵢⱼnext = ϵⁱᵢⱼ + ϵ̇ⁱᵢⱼ*Δt

    # damage growth out
    if damage_growth_out
        KIᵒ = compute_KI(r,(σᵒᵢⱼnext+σᵒᵢⱼ)/2,Dᵒ)
        Ḋᵒ = compute_subcrit_damage_rate(r,KIᵒ,Dᵒ)
        Dᵒnext = Dᵒ + Ḋᵒ*Δt
    else
        Dᵒnext = Dᵒ
    end
    # damage growth in
    KIⁱ = compute_KI(r,(σⁱᵢⱼnext+σⁱᵢⱼ)/2,Dⁱ) # compute KI at intermediate stress : TO TEST.
    Ḋⁱ = compute_subcrit_damage_rate(r,KIⁱ,Dⁱ)
    Dⁱnext = Dⁱ + Ḋⁱ*Δt
    #_ , Ḋⁱ2 = compute_ϵ̇ij(r,Dⁱ,σⁱᵢⱼ,σⁱᵢⱼnext,Δt) # compute damage growth inside the shear band
    #@assert Ḋⁱ == Ḋⁱ2 # shouldn't error out, then remove preceeding line

    return Snext, σⁱᵢⱼnext, σᵒᵢⱼnext, ϵⁱᵢⱼnext, Dⁱnext, Dᵒnext, u

end

function two_points_update_du(u::SVector{8,T},p,t) where T<:Real
    # unpack
    S, σⁱξξ, σⁱoop, ϵⁱηη, ϵⁱηη_plas, ϵⁱξη_plas, Dᵒ, Dⁱ = u
    r = p.r
    mp = p.mp
    Ṡ, σ̇ⁱξξ, σ̇ⁱoop, ϵ̇ⁱηη, Ḋᵒ, Ḋⁱ = p.du
    σ₃, ϵ̇ⁱξη, θ = p.scalars
    flags = mp.flags

    # debug
    (Dᵒ < r.D₀) && @show(Ḋᵒ, Dᵒ)

    # use static arrays in constructors to avoid allocation
    σᵒᵢⱼ_p = build_principal_stress_tensor(r,S,σ₃,Dᵒ ; abstol=1e-15)
    σᵒᵢⱼ = band_coords(σᵒᵢⱼ_p,θ)
    σⁱᵢⱼ = SymmetricTensor{2,3}(SA[σⁱξξ σᵒᵢⱼ[1,2] zero(σⁱξξ) ; σᵒᵢⱼ[1,2] σᵒᵢⱼ[2,2] zero(σⁱξξ) ; zero(σⁱξξ) zero(σⁱξξ) σⁱoop])

    # compute damage rates :
    Ḋᵒ = compute_subcrit_damage_rate(r, compute_KI(r,σᵒᵢⱼ,Dᵒ), Dᵒ)
    Ḋⁱ = compute_subcrit_damage_rate(r, compute_KI(r,σⁱᵢⱼ,Dⁱ), Dⁱ)

    #((Ḋᵒ < 0) | (Ḋⁱ < 0)) && (p.terminate_flag = true)
    #@show Dᵒ Dⁱ
    #@show Ḋᵒ Ḋⁱ
    #println()

    # vector to be solved by non linear iterations (scaled in order to homogenize componant magnitude)
    du_nl = SA_F64[Ṡ*σ₃, σ̇ⁱξξ, σ̇ⁱoop, ϵ̇ⁱηη*r.G]
    ∇res = du_nl*du_nl'
    for i in 1:mp.solver.newton_maxiter
        # get residual and its gradient with respect to du_nl

        # try
        #     ∇res = ForwardDiff.jacobian(du_nl -> residual_2_points_2(du_nl, σᵒᵢⱼ, σⁱᵢⱼ, Dᵒ, Ḋᵒ, Dⁱ, Ḋⁱ, p), du_nl)
        # catch e
        #     @warn "Jacobian calculation generates nans.\n Exiting newton iter $(i) \n"
        #     throw(e) # comment this line to continue without erroring
        #     flags.nan = true
        #     break
        # end
        ∇res = ForwardDiff.jacobian(du_nl -> residual_2_points_2(du_nl, σᵒᵢⱼ, σⁱᵢⱼ, Dᵒ, Ḋᵒ, Dⁱ, Ḋⁱ, p), du_nl)
        res  = residual_2_points_2(du_nl, σᵒᵢⱼ, σⁱᵢⱼ, Dⁱ, Ḋⁱ, Dᵒ, Ḋᵒ, p)

        # newton correction term
        δdu_nl = - ∇res\res
        if any(isnan.(δdu_nl))
            @warn "Solution vector update δu contains nans : $(δu).\n Exiting newton iter $(i) with residual norm of $(norm(res)) \n"
            flags.nan = true
            break
        else # update if no NaNs
            du_nl = du_nl + δdu_nl
        end
        #@debug "δu = $δu"
        #@debug "typeof(u) = $(typeof(u))"

        (norm(res) <= mp.solver.newton_abstol) && (@debug("Newton iter $i ending norm res = $(norm(res))") ; break)

        if i == mp.solver.newton_maxiter
            #@debug("Newton maxiter ending norm res = $(norm(res))")#("max newton iteration reached ($i), residual still higher than abstol with $(norm(res))")
            #@debug @show(du_nl)
        end
    end
    # rescale appropriate components
    Ṡ = du_nl[1] / σ₃
    σ̇ⁱξξ = du_nl[2]
    σ̇ⁱoop = du_nl[3]
    ϵ̇ⁱηη = du_nl[4]/r.G

    # compute plastic strain
    σ̇ᵒᵢⱼ = compute_rotated_stress_rate_from_band_coords_2(r,Ṡ,σ₃,σᵒᵢⱼ,Dᵒ,Ḋᵒ,θ; damaged_allowed=p.allow_Ḋᵒ)
    σ̇ⁱᵢⱼ = SymmetricTensor{2,3}(SA[σ̇ⁱξξ σ̇ᵒᵢⱼ[1,2] zero(σ̇ⁱξξ) ; σ̇ᵒᵢⱼ[1,2] σ̇ᵒᵢⱼ[2,2] zero(σ̇ⁱξξ) ; zero(σ̇ⁱξξ) zero(σ̇ⁱξξ) σ̇ⁱoop])
    ϵⁱᵢⱼ = compute_ϵ̇ij_2(r,Dⁱ,Ḋⁱ,σⁱᵢⱼ,σ̇ⁱᵢⱼ; damaged_allowed=true)
    ϵⁱᵢⱼ_elast = compute_ϵ̇ij_2(r,Dⁱ,Ḋⁱ,σⁱᵢⱼ,σ̇ⁱᵢⱼ; damaged_allowed=false)
    ϵⁱᵢⱼ_plas = ϵⁱᵢⱼ - ϵⁱᵢⱼ_elast

    # build full derivatives vector
    du = SA_F64[Ṡ, σ̇ⁱξξ, σ̇ⁱoop, ϵ̇ⁱηη, ϵⁱᵢⱼ_plas[2,2], ϵⁱᵢⱼ_plas[1,2], Ḋᵒ, Ḋⁱ]

    (p.du isa Vector) && (p.du .= du)
    (p.du isa SVector) && (p.du = du)
    return du
end

function two_points_update_du(u::SVector{6,T},p,t) where T<:Real
    # unpack
    S, σⁱξξ, σⁱoop, ϵⁱηη, Dᵒ, Dⁱ = u
    r = p.r
    mp = p.mp
    Ṡ, σ̇ⁱξξ, σ̇ⁱoop, ϵ̇ⁱηη, Ḋᵒ, Ḋⁱ = p.du
    σ₃, ϵ̇ⁱξη, θ = p.scalars
    flags = mp.flags

    # debug
    (Dᵒ < r.D₀) && @show(Ḋᵒ, Dᵒ)

    # use static arrays in constructors to avoid allocation
    σᵒᵢⱼ_p = build_principal_stress_tensor(r,S,σ₃,Dᵒ ; abstol=1e-15)
    σᵒᵢⱼ = band_coords(σᵒᵢⱼ_p,θ)
    σⁱᵢⱼ = SymmetricTensor{2,3}(SA[σⁱξξ σᵒᵢⱼ[1,2] zero(σⁱξξ) ; σᵒᵢⱼ[1,2] σᵒᵢⱼ[2,2] zero(σⁱξξ) ; zero(σⁱξξ) zero(σⁱξξ) σⁱoop])

    # compute damage rates :
    Ḋᵒ = compute_subcrit_damage_rate(r, compute_KI(r,σᵒᵢⱼ,Dᵒ), Dᵒ)
    Ḋⁱ = compute_subcrit_damage_rate(r, compute_KI(r,σⁱᵢⱼ,Dⁱ), Dⁱ)

    ((Ḋᵒ < 0) | (Ḋⁱ < 0)) && (p.terminate_flag = true)
    #@show Dᵒ Dⁱ
    #@show Ḋᵒ Ḋⁱ
    #println()

    # vector to be solved by non linear iterations (scaled in order to homogenize componant magnitude)
    du_nl = SA_F64[Ṡ*σ₃, σ̇ⁱξξ, σ̇ⁱoop, ϵ̇ⁱηη*r.G]
    ∇res = du_nl*du_nl'
    for i in 1:mp.solver.newton_maxiter
        # get residual and its gradient with respect to du_nl

        # try
        #     ∇res = ForwardDiff.jacobian(du_nl -> residual_2_points_2(du_nl, σᵒᵢⱼ, σⁱᵢⱼ, Dᵒ, Ḋᵒ, Dⁱ, Ḋⁱ, p), du_nl)
        # catch e
        #     @warn "Jacobian calculation generates nans.\n Exiting newton iter $(i) \n"
        #     throw(e) # comment this line to continue without erroring
        #     flags.nan = true
        #     break
        # end
        ∇res = ForwardDiff.jacobian(du_nl -> residual_2_points_2(du_nl, σᵒᵢⱼ, σⁱᵢⱼ, Dᵒ, Ḋᵒ, Dⁱ, Ḋⁱ, p), du_nl)
        res  = residual_2_points_2(du_nl, σᵒᵢⱼ, σⁱᵢⱼ, Dⁱ, Ḋⁱ, Dᵒ, Ḋᵒ, p)

        # newton correction term
        δdu_nl = - ∇res\res
        if any(isnan.(δdu_nl))
            @warn "Solution vector update δu contains nans : $(δu).\n Exiting newton iter $(i) with residual norm of $(norm(res)) \n"
            flags.nan = true
            break
        else # update if no NaNs
            du_nl = du_nl + δdu_nl
        end
        #@debug "δu = $δu"
        #@debug "typeof(u) = $(typeof(u))"

        (norm(res) <= mp.solver.newton_abstol) && (@debug("Newton iter $i ending norm res = $(norm(res))") ; break)

        if i == mp.solver.newton_maxiter
            #@debug("Newton maxiter ending norm res = $(norm(res))")#("max newton iteration reached ($i), residual still higher than abstol with $(norm(res))")
            #@debug @show(du_nl)
        end
    end
    du = SA_F64[du_nl[1]/σ₃, du_nl[2], du_nl[3], du_nl[4]/r.G, Ḋᵒ, Ḋⁱ]

    (p.du isa Vector) && (p.du .= du)
    (p.du isa SVector) && (p.du = du)
    return du
end

function residual_2_points_2(u, σᵒᵢⱼ, σⁱᵢⱼ, Dᵒ, Ḋᵒ, Dⁱ, Ḋⁱ, p)
    Ṡ, σ̇ⁱξξ, σ̇ⁱoop, ϵ̇ⁱηη = u
    σ₃, ϵ̇ⁱξη, θ = p.scalars
    r = p.r

    # rescale appropriate components
    Ṡ = Ṡ/σ₃
    ϵ̇ⁱηη = ϵ̇ⁱηη/r.G

    # state out
    σ̇ᵒᵢⱼ = compute_rotated_stress_rate_from_band_coords_2(r,Ṡ,σ₃,σᵒᵢⱼ,Dᵒ,Ḋᵒ,θ; damaged_allowed=p.allow_Ḋᵒ)
    ϵ̇ᵒᵢⱼ = compute_ϵ̇ij_2(r,Dᵒ,Ḋᵒ,σᵒᵢⱼ,σ̇ᵒᵢⱼ; damaged_allowed=p.allow_Ḋᵒ)

    # state in
    σ̇ⁱᵢⱼ = SymmetricTensor{2,3}(SA[σ̇ⁱξξ σ̇ᵒᵢⱼ[1,2] zero(σ̇ⁱξξ) ; σ̇ᵒᵢⱼ[1,2] σ̇ᵒᵢⱼ[2,2] zero(σ̇ⁱξξ) ; zero(σ̇ⁱξξ) zero(σ̇ⁱξξ) σ̇ⁱoop])
    #σ̇ⁱᵢⱼ = SymmetricTensor{typeof(σ̇ⁱξξ),2,3}(SA[σ̇ⁱξξ σ̇ᵒᵢⱼ[1,2] 0 ; σ̇ᵒᵢⱼ[1,2] σ̇ᵒᵢⱼ[2,2] 0 ; 0 0 σ̇ⁱoop])
    ϵ̇ⁱᵢⱼ = compute_ϵ̇ij_2(r,Dⁱ,Ḋⁱ,σⁱᵢⱼ,σ̇ⁱᵢⱼ; damaged_allowed=true)

    res = SA[ϵ̇ᵒᵢⱼ[1,1] - ϵ̇ⁱᵢⱼ[1,1],
             ϵ̇ⁱηη - ϵ̇ⁱᵢⱼ[2,2],
             -ϵ̇ⁱᵢⱼ[3,3],
             ϵ̇ⁱξη - ϵ̇ⁱᵢⱼ[1,2]] / ϵ̇ⁱξη

    return res
end

function residual_2_points(r,S,σ₃,Dⁱ,Dᵒ,ϵⁱᵢⱼ,σⁱᵢⱼ,σᵒᵢⱼ,ϵ̇ⁱξη,Δt,θ,u ; damage_growth_out=true)

    Ṡ, σ̇ᵒᵢⱼ, σ̇ⁱᵢⱼ, ϵ̇ⁱᵢⱼ = compute_stress_strain_derivatives_from_u(r,S,σ₃,Dⁱ,Dᵒ,ϵⁱᵢⱼ,σⁱᵢⱼ,σᵒᵢⱼ,ϵ̇ⁱξη,Δt,θ,u; damage_growth_out)
    σⁱᵢⱼnext = σⁱᵢⱼ + σ̇ⁱᵢⱼ*Δt

    ϵ̇ⁱᵢⱼ_analytical, _ = compute_ϵ̇ij(r,Dⁱ,σⁱᵢⱼ,σⁱᵢⱼnext,Δt ; damaged_allowed=true)
    Δϵ̇ⁱᵢⱼ = ϵ̇ⁱᵢⱼ_analytical - ϵ̇ⁱᵢⱼ
    res = Vec(Δϵ̇ⁱᵢⱼ[1,1]* 2r.G, Δϵ̇ⁱᵢⱼ[2,2]* 2r.G, Δϵ̇ⁱᵢⱼ[3,3]* 2r.G, Δϵ̇ⁱᵢⱼ[1,2]* 2r.G)
    return res
end

function compute_stress_strain_derivatives_from_u(r,S,σ₃,Dⁱ,Dᵒ,ϵⁱᵢⱼ,σⁱᵢⱼ,σᵒᵢⱼ,ϵ̇ⁱξη,Δt,θ,u ; damage_growth_out=true)
    Snext, σⁱξξnext, σⁱoopnext, ϵⁱηηnext = u
    σⁱξξ = σⁱᵢⱼ[1,1]
    σⁱoop = σⁱᵢⱼ[3,3]
    ϵⁱηη = ϵⁱᵢⱼ[2,2]
    ϵⁱξη = ϵⁱᵢⱼ[1,2]
    Ṡ = (Snext - S) / Δt
    ϵⁱξηnext = ϵⁱξη + ϵ̇ⁱξη*Δt

    if isnan(Ṡ)
        println("Snext : ",Snext.value)
        println("S : ",S)
        println("Δt : ",Δt)
    end

    σ̇ᵒᵢⱼ = compute_rotated_stress_rate_from_band_coords(r,Ṡ,σ₃,σᵒᵢⱼ,Dᵒ,Δt,θ ; damaged_allowed=damage_growth_out) # no damage outside the band
    σ̇ⁱᵢⱼ = get_σ̇ⁱij_from_primitive_vars(σ̇ᵒᵢⱼ,σⁱξξ,σⁱoop,σⁱξξnext,σⁱoopnext,Δt)
    #σ̇ⁱᵢⱼ = set_plane_strain_oop_stress_rate(σⁱᵢⱼ,σ̇ⁱᵢⱼ,r,Dⁱ,Δt ; damaged_allowed=true)
    ϵ̇ⁱᵢⱼ = get_ϵ̇ⁱᵢⱼ_from_primitive_vars(r,σᵒᵢⱼ,σ̇ᵒᵢⱼ,ϵⁱηη,ϵⁱξη,Dⁱ,ϵⁱηηnext,ϵⁱξηnext,Δt ; damage_growth_out)
    return Ṡ, σ̇ᵒᵢⱼ, σ̇ⁱᵢⱼ, ϵ̇ⁱᵢⱼ
end

function get_σ̇ⁱij_from_primitive_vars(σ̇ᵒᵢⱼ_band,σⁱξξ,σⁱoop,σⁱξξnext,σⁱoopnext,Δt)
    σ̇ⁱξξ = (σⁱξξnext - σⁱξξ)/Δt
    σ̇ⁱoop = (σⁱoopnext - σⁱoop)/Δt
    σ̇ⁱij_band = SymmetricTensor{2,3}([    σ̇ⁱξξ       σ̇ᵒᵢⱼ_band[1,2]     0    ;
                                        σ̇ᵒᵢⱼ_band[2,1] σ̇ᵒᵢⱼ_band[2,2]     0    ;
                                            0              0          σ̇ⁱoop   ])
    return σ̇ⁱij_band
end

function get_ϵ̇ⁱᵢⱼ_from_primitive_vars(r,σᵒᵢⱼ,σ̇ᵒᵢⱼ,ϵⁱηη,ϵⁱξη,D,ϵⁱηηnext,ϵⁱξηnext,Δt ; damage_growth_out=true)
    σᵒᵢⱼnext = σᵒᵢⱼ + σ̇ᵒᵢⱼ*Δt
    ϵ̇ⁱηη = (ϵⁱηηnext - ϵⁱηη) /Δt
    ϵ̇ⁱξη = (ϵⁱξηnext - ϵⁱξη) /Δt
    ϵ̇ᵒᵢⱼ, _ = compute_ϵ̇ij(r,D,σᵒᵢⱼ,σᵒᵢⱼnext,Δt ; damaged_allowed=damage_growth_out)
    ϵ̇ⁱᵢⱼ = SymmetricTensor{2,3}([ ϵ̇ᵒᵢⱼ[1,1]   ϵ̇ⁱξη   0  ;
                                    ϵ̇ⁱξη      ϵ̇ⁱηη   0  ;
                                        0         0     0  ]) # ϵ̇ᵒᵢⱼ[3,3] put zero on last term if needed to force plain strain
    return ϵ̇ⁱᵢⱼ
end

function compute_rotated_stress_rate_from_principal_coords(r,Ṡ,σ₃,σᵢⱼ_principal,D,Δt,θ ; damaged_allowed=true)
    σ̇ᵢⱼ_principal = SymmetricTensor{2,3}([Ṡ*σ₃ 0 0 ; 0 0 0 ; 0 0 r.ν*Ṡ*σ₃])
    σ̇ᵢⱼ_principal = set_plane_strain_oop_stress_rate(σᵢⱼ_principal,σ̇ᵢⱼ_principal,r,D,Δt ; abstol=1e-16, damaged_allowed)
    σ̇ᵢⱼ_band = band_coords(σ̇ᵢⱼ_principal,θ)
    return σ̇ᵢⱼ_band
end
function compute_rotated_stress_rate_from_band_coords(r,Ṡ,σ₃,σᵢⱼ_band,D,Δt,θ; damaged_allowed=true)
    σ̇ᵢⱼ_band = compute_rotated_stress_rate_guess(r,Ṡ,σ₃,θ)
    #any(isnan.(σ̇ᵢⱼ_band)) && @error("isnan here !")
    σ̇ᵢⱼ_band = set_plane_strain_oop_stress_rate(σᵢⱼ_band,σ̇ᵢⱼ_band,r,D,Δt ; abstol=1e-16, damaged_allowed)
    return σ̇ᵢⱼ_band
end

function compute_rotated_stress_rate_from_band_coords_2(r,Ṡ,σ₃,σᵢⱼ_band,D,Ḋ,θ; damaged_allowed=true)
    σ̇ᵢⱼ_band = compute_rotated_stress_rate_guess(r,Ṡ,σ₃,θ)
    #any(isnan.(σ̇ᵢⱼ_band)) && @error("isnan here !")
    σ̇ᵢⱼ_band = set_plane_strain_oop_stress_rate_2(r,σᵢⱼ_band,σ̇ᵢⱼ_band,D,Ḋ ; abstol=1e-16, maxiter=100, damaged_allowed)
    return σ̇ᵢⱼ_band
end

function compute_rotated_stress_rate_guess(r,Ṡ,σ₃,θ)
    if isnan(Ṡ)
        throw(@error("isnan here !"))
    end
    σ̇ᵢⱼ_principal = SymmetricTensor{2,3}(SA[Ṡ*σ₃ 0 0 ; 0 0 0 ; 0 0 r.ν*Ṡ*σ₃])
    σ̇ᵢⱼ_band = band_coords(σ̇ᵢⱼ_principal,θ)
    return σ̇ᵢⱼ_band
end


text:
"""
    qitf_hamiltonian(lx, g, J)

generate the quantum Ising in a transverse field Haimltonian that is H
= J ∑_i Z_i Z_{i+1} + g ∑_i X_i where `lx` is the size of the
system, `J` ising exchange amplitude and `g` is the strength of the
transverse field. The phase transition happens at `g=J` which is the
default value.

"""
function qitf_hamiltonian(lx::Int, g::Float64=1.0, J::Float64=-1.0;
                          envelope::Vector{Float64} = ones(lx),
                          boundary::Symbol=:OBC)
    sz, sp ,sm = spinoperators(1/2)
    X = sp + sm
    Z = 2 * sz
    bond_term = J * (Z ⊗ Z) + g/2. * (X ⊗ I(2) + I(2) ⊗ X)
    Hmat = envelope[1] * bond_term
    if lx > 2
        for i=2:lx-1
            Hmat = Hmat ⊗ I(2) + envelope[i] * (I(2^(i-1)) ⊗ bond_term)
        end
    end

    # Should the below line be added only for periodic?! Nahhh...
    Hmat += envelope[1] * g/2 * (X ⊗ I(2^(lx-1)) + I(2^(lx-1)) ⊗ X)

    if boundary == :OBC
        return Hmat
    elseif boundary == :PBC
        Hmat += envelope[lx] * J * (Z ⊗ I(2^(lx-2)) ⊗ Z)
        return Hmat
    else
        throw("unrecognized boundary condition : $boundary")
    end
end

function qitf_bondtensor(g::Float64=1.0, J::Float64=-1.0)
    sz, sp ,sm = spinoperators(1/2)
    X = sp + sm
    Z = 2 * sz
    bond_term = J * (Z ⊗ Z) + g/2. * (X ⊗ I(2) + I(2) ⊗ X)
    reshape(Matrix(bond_term), 2, 2, 2, 2)
end

function qitf_energy_exact(g::Float64=1.0, J::Float64=-1.0)
    e_exact, e_exact_err = quadgk(k -> -sqrt(1+g^2-2*g*cos(k))/pi, 0, pi)
    e_exact
end

max_stars_repo_name: amiragha/mps

max_issues_repo_name: amiragha/mps

max_forks_repo_name: amiragha/mps


text:
module TestErrorHandling

using Test
using ThreadPools

function testfunc(err)
    fn(x) = x == err && error("Error")
    return fn
end

@testset "exception handling" begin

    @testset "foreach functions" begin
        for i in 1:16
            @test_throws TaskFailedException tforeach(testfunc(i), 1:16)
            @test_throws TaskFailedException bforeach(testfunc(i), 1:16)
            @test_throws TaskFailedException qforeach(testfunc(i), 1:16)
            @test_throws TaskFailedException qbforeach(testfunc(i), 1:16)
        end
    end

    @testset "map functions" begin
        for i in 1:16
            @test_throws TaskFailedException tmap(testfunc(i), 1:16)
            @test_throws TaskFailedException bmap(testfunc(i), 1:16)
            @test_throws TaskFailedException qmap(testfunc(i), 1:16)
            @test_throws TaskFailedException qbmap(testfunc(i), 1:16)
        end
    end

    @testset "macros" begin
        for i in 1:16
            @test_throws TaskFailedException begin
                @bthreads for j in 1:16
                    testfunc(i)(j)
                end
            end
            @test_throws TaskFailedException begin
                @qthreads for j in 1:16
                    testfunc(i)(j)
                end
            end
            @test_throws TaskFailedException begin
                @qbthreads for j in 1:16
                    testfunc(i)(j)
                end
            end
        end
    end
end

@testset "stop on error" begin
    for fn in [tforeach, bforeach, qforeach, qbforeach, tmap, bmap, qmap, qbmap]
        let
            count = Threads.Atomic{Int}(0)
            function testfn(x)
                if x == 2
                    sleep(0.001)
                    error("Error")
                else
                    sleep(0.02)
                    Threads.atomic_add!(count, 1)
                end
            end
            @test_throws TaskFailedException fn(testfn, 1:16)
            @test count[] <= Threads.nthreads()
        end
    end
end

end


text:
using DashBootstrapComponents

dropdown = dbc_dropdownmenu(
    label = "Menu",
    menu_variant = "dark",
    children = [
        dbc_dropdownmenuitem("Item 1"),
        dbc_dropdownmenuitem("Item 2"),
        dbc_dropdownmenuitem("Item 3"),
    ],
)

max_stars_repo_name: glsdown/dash-bootstrap-components

max_issues_repo_name: glsdown/dash-bootstrap-components

max_forks_repo_name: glsdown/dash-bootstrap-components


text:
@testset "$(format("Loss: Poisson[construction]"))" begin
    @test typeof(Loss{Poisson}(Poisson())) == Loss{Poisson}
    @test typeof(Loss(Poisson())) == Loss{Poisson}
    @test typeof(Loss(:Poisson)) == Loss{Poisson}
end

@testset "$(format("Optimizer: Poisson Loss [Small][Forgiving][AutoGrad]"))"  begin
    @timeit to "Optimizer: Poisson Loss [Small][Forgiving][AutoGrad]"  begin
        let
            is_admissable(result)  = begin
                if result["relative-error[#within-radius(1e-5)]"] < 0.1 || result["relative-error[#within-radius(1)]"] <0.1
                    return true
                end
                @warn @sprintf("expected %f, got %f",0.1,result["relative-error[#within-radius(1e-5)]"])
                return false
            end
            @test !isnothing(provide(Loss(Poisson())))
            @test !isnothing(provide(Loss(:Poisson)))
            @test is_admissable(unit_test_train_subloss(gradient_eval      = provide(Loss(Poisson())),
                                                        input_distribution = Distributions.Poisson(5),
                                                        input_size         = 1000,
                                                        ρ                  = 0,
                                                        step_size          = 0.1,
                                                        max_iter           = 100))
            @test is_admissable(unit_test_train_subloss(gradient_eval      = provide(Loss(Poisson())),
                                                        input_distribution = Distributions.Poisson(5),
                                                        input_size         = 3000,
                                                        ρ                  = 0,
                                                        step_size          = 0.1,
                                                        max_iter          = 100))
            @test is_admissable(unit_test_train_subloss(gradient_eval      = provide(Loss(Poisson())),
                                                        input_distribution = Distributions.Poisson(5),
                                                        input_size         = 5000,
                                                        ρ                  = 0,
                                                        step_size          = 0.1,
                                                        max_iter           = 100))
            @test is_admissable(unit_test_train_subloss(gradient_eval      = provide(Loss(Poisson())),
                                                        input_distribution = Distributions.Poisson(5),
                                                        input_size         = 1000,
                                                        ρ                  = 0.1,
                                                        step_size          = 0.1,
                                                        max_iter           = 100))
            @test is_admissable(unit_test_train_subloss(gradient_eval      = provide(Loss(Poisson())),
                                                        input_distribution = Distributions.Poisson(5),
                                                        input_size         = 3000,
                                                        ρ                  = 0.1,
                                                        step_size          = 0.1,
                                                        max_iter = 100))
            @test is_admissable(unit_test_train_subloss(gradient_eval      = provide(Loss(Poisson())),
                                                        input_distribution = Distributions.Poisson(5),
                                                        input_size         = 5000,
                                                        ρ                  = 0.2,
                                                        step_size          = 0.1,
                                                        max_iter           = 100))
            @test is_admissable(unit_test_train_subloss(gradient_eval      = provide(Loss(:Poisson)),
                                                        input_distribution = Distributions.Poisson(5),
                                                        input_size         = 1000,
                                                        ρ                  = 0,
                                                        step_size          = 0.1,
                                                        max_iter           = 100))
            @test is_admissable(unit_test_train_subloss(gradient_eval      = provide(Loss(:Poisson)),
                                                        input_distribution = Distributions.Poisson(5),
                                                        input_size         = 3000,
                                                        ρ                  = 0,
                                                        step_size          = 0.1,
                                                        max_iter          = 100))
            @test is_admissable(unit_test_train_subloss(gradient_eval      = provide(Loss(:Poisson)),
                                                        input_distribution = Distributions.Poisson(5),
                                                        input_size         = 5000,
                                                        ρ                  = 0,
                                                        step_size          = 0.1,
                                                        max_iter           = 100))
            @test is_admissable(unit_test_train_subloss(gradient_eval      = provide(Loss(:Poisson)),
                                                        input_distribution = Distributions.Poisson(5),
                                                        input_size         = 1000,
                                                        ρ                  = 0.1,
                                                        step_size          = 0.1,
                                                        max_iter           = 100))
            @test is_admissable(unit_test_train_subloss(gradient_eval      = provide(Loss(:Poisson)),
                                                        input_distribution = Distributions.Poisson(5),
                                                        input_size         = 3000,
                                                        ρ                  = 0.1,
                                                        step_size          = 0.1,
                                                        max_iter = 100))
            @test is_admissable(unit_test_train_subloss(gradient_eval      = provide(Loss(:Poisson)),
                                                        input_distribution = Distributions.Poisson(5),
                                                        input_size         = 5000,
                                                        ρ                  = 0.2,
                                                        step_size          = 0.1,
                                                        max_iter           = 100))
        end
    end
end

@testset "$(format("Optimizer: Poisson Loss [Small][Forgiving][Native]"))" begin
    @timeit to "Optimizer: Poisson Loss [Small][Forgiving][Native]" begin
        is_admissable(result)  = begin
                if result["relative-error[#within-radius(1e-5)]"] < 0.1 || result["relative-error[#within-radius(1)]"] <0.1
                    return true
                end
                @warn @sprintf("expected %f, got %f",0.1,result["relative-error[#within-radius(1e-5)]"])
                return false
            end
        @test is_admissable(unit_test_train_subloss(gradient_eval      = Loss(Poisson()),
                                                    input_distribution = Distributions.Poisson(5),
                                                    input_size         = 1000,
                                                    ρ                  = 0,
                                                    step_size          = 0.1,
                                                    max_iter           = 100))
        @test is_admissable(unit_test_train_subloss(gradient_eval      = Loss(Poisson()),
                                                    input_distribution = Distributions.Poisson(5),
                                                    input_size         = 3000,
                                                    ρ                  = 0,
                                                    step_size          = 0.1,
                                                    max_iter          = 100))
        @test is_admissable(unit_test_train_subloss(gradient_eval      = Loss(Poisson()),
                                                    input_distribution = Distributions.Poisson(5),
                                                    input_size         = 5000,
                                                    ρ                  = 0,
                                                    step_size          = 0.1,
                                                    max_iter           = 100))
        @test is_admissable(unit_test_train_subloss(gradient_eval      = Loss(Poisson()),
                                                    input_distribution = Distributions.Poisson(5),
                                                    input_size         = 1000,
                                                    ρ                  = 0.2,
                                                    step_size          = 0.1,
                                                    max_iter           = 100))
        @test is_admissable(unit_test_train_subloss(gradient_eval      = Loss(Poisson()),
                                                    input_distribution = Distributions.Poisson(5),
                                                    input_size         = 3000,
                                                    ρ                  = 0.2,
                                                    step_size          = 0.1,
                                                    max_iter = 100))
        @test is_admissable(unit_test_train_subloss(gradient_eval      = Loss(Poisson()),
                                                    input_distribution = Distributions.Poisson(5),
                                                    input_size         = 5000,
                                                    ρ                  = 0.2,
                                                    step_size          = 0.1,
                                                    max_iter           = 100))
        @test is_admissable(unit_test_train_subloss(gradient_eval      = Loss(:Poisson),
                                                    input_distribution = Distributions.Poisson(5),
                                                    input_size         = 1000,
                                                    ρ                  = 0,
                                                    step_size          = 0.1,
                                                    max_iter           = 100))
        @test is_admissable(unit_test_train_subloss(gradient_eval      = Loss(:Poisson),
                                                    input_distribution = Distributions.Poisson(5),
                                                    input_size         = 3000,
                                                    ρ                  = 0,
                                                    step_size          = 0.1,
                                                    max_iter          = 100))
        @test is_admissable(unit_test_train_subloss(gradient_eval      = Loss(:Poisson),
                                                    input_distribution = Distributions.Poisson(5),
                                                    input_size         = 5000,
                                                    ρ                  = 0,
                                                    step_size          = 0.1,
                                                    max_iter           = 100))
        @test is_admissable(unit_test_train_subloss(gradient_eval      = Loss(:Poisson),
                                                    input_distribution = Distributions.Poisson(5),
                                                    input_size         = 1000,
                                                    ρ                  = 0.2,
                                                    step_size          = 0.1,
                                                    max_iter           = 100))
        @test is_admissable(unit_test_train_subloss(gradient_eval      = Loss(:Poisson),
                                                    input_distribution = Distributions.Poisson(5),
                                                    input_size         = 3000,
                                                    ρ                  = 0.2,
                                                    step_size          = 0.1,
                                                    max_iter = 100))
        @test is_admissable(unit_test_train_subloss(gradient_eval      = Loss(:Poisson),
                                                    input_distribution = Distributions.Poisson(5),
                                                    input_size         = 5000,
                                                    ρ                  = 0.2,
                                                    step_size          = 0.1,
                                                    max_iter           = 100))
    end
end

#include("test_impl_poissonloss.jl")
#include("test_impl_gammaloss.jl")

@label END_OF_TEST


text:
eval(quote
    import .Polyhedra: polyhedron
    export polyhedron
    using .Polyhedra: HRep, VRep,
                      removehredundancy!, removevredundancy!
    import JuMP, GLPK

    function default_polyhedra_backend(P, N::Type{<:AbstractFloat})
        return Polyhedra.default_library(LazySets.dim(P), N)
    end

    function default_polyhedra_backend(P, N::Type{<:Rational})
        return Polyhedra.default_library(LazySets.dim(P), N)
    end

    # NOTE: exists in parallel to `default_lp_solver` because we use different
    # interfaces (see #1493)
    function default_lp_solver_polyhedra(N::Type{<:AbstractFloat})
        return JuMP.with_optimizer(GLPK.Optimizer)
    end

    # NOTE: exists in parallel to `default_lp_solver` because we use different
    # interfaces (see #1493)
    function default_lp_solver_polyhedra(N::Type{<:Rational})
        return JuMP.with_optimizer(GLPK.Optimizer, method=GLPK.EXACT)
    end
end)

eval(load_polyhedra_hpolytope())
eval(load_polyhedra_hpolyhedron())
eval(load_polyhedra_vpolytope())

eval(initialize_mesh())


text:
"""
	ArrayFuse{AT, T, P} <: AbstractArray{T, 1}

GPU Friendly type to wrap around two arrays - `visible` and `hidden`, for which when we `setindex!` some value `v` at index `i`
we get

```
visible[i] = p[1] * visible[i] + p[2] * v
hidden[i] = hidden[i] + p[3] * visible[i]
```

where p is a parameter tuple of size 3.
"""
struct ArrayFuse{AT, T, P} <: AbstractArray{T, 1}
	visible::AT
	hidden::AT
	p::P
end

ArrayFuse(visible::AT, hidden::AT, p) where {AT} = ArrayFuse{AT, eltype(visible), typeof(p)}(visible, hidden, p)

@inline function Base.copyto!(af::OrdinaryDiffEq.ArrayFuse{AT, T, P}, src::Base.Broadcast.Broadcasted) where {AT, T, P}
	@. af.visible = af.p[1] * af.visible + af.p[2] * src
	@. af.hidden = af.hidden + af.p[3] * af.visible
end

@inline function Base.copyto!(af::OrdinaryDiffEq.ArrayFuse{AT, T, P}, src::Base.Broadcast.Broadcasted{F1, Axes, F, Args}) where {AT, T, P, F1<:Base.Broadcast.AbstractArrayStyle{0}, Axes, F, Args<:Tuple}
	@. af.visible = af.p[1] * af.visible + af.p[2] * src
	@. af.hidden = af.hidden + af.p[3] * af.visible
end

# not recommended but good to have
@inline function Base.getindex(af::ArrayFuse, index)
	return af.visible[index]
end

@inline function Base.setindex!(af::ArrayFuse, value, index)
	af.visible[index] = af.p[1] * af.visible[index] + af.p[2] * value
	af.hidden[index] =  muladd(af.p[3], af.visible[index], af.hidden[index])
end

@inline Base.size(af::ArrayFuse) = length(af.visible)
@inline Base.axes(af::ArrayFuse) = axes(af.visible)

Adapt.adapt_structure(to, af::ArrayFuse{AT, T, P}) where {AT, T, P} =
	ArrayFuse(adapt(to, af.visible), adapt(to, af.hidden), af.p)


text:
# div σ = f
# div u = g
# σ = -p I + 2 η ε(u)


text:
#"""
#"""
#function writevtk(
#  trian::Grid, filebase; celldata=Dict(), nodaldata=Dict())
#  write_vtk_file(trian,filebase,celldata=celldata,nodaldata=nodaldata)
#end

"""
    writevtk(reffe::LagrangianRefFE,filebase)
"""
function writevtk(reffe::LagrangianRefFE,filebase)

  p = get_polytope(reffe)
  writevtk(p,filebase)

  node_coords = get_node_coordinates(reffe)
  node_comp_to_dof = get_node_and_comp_to_dof(reffe)
  nodaldata = [
    "dof" => node_comp_to_dof,
    "node" => collect(1:num_nodes(reffe))]
  writevtk(node_coords, "$(filebase)_nodes"; nodaldata = nodaldata)

end

"""
    writevtk(x::AbstractVector{<:Point}, filebase; kwargs...)
"""
function writevtk(x::AbstractVector{<:Point}, filebase; kwargs...)
  vtkfile = createvtk(x,filebase; kwargs...)
  outfiles = vtk_save(vtkfile)
end

"""
    createvtk(x::AbstractVector{<:Point}, filebase; kwargs...)
"""
function createvtk(x::AbstractVector{<:Point}, filebase; kwargs...)
  grid = UnstructuredGrid(x)
  create_vtk_file(grid,filebase; kwargs...)
end

"""
    writevtk(p::Polytope,filebase)
"""
function writevtk(p::Polytope,filebase)
  for d in 0:(num_dims(p)-1)
    grid = Grid(ReferenceFE{d},p)
    write_vtk_file(grid,"$(filebase)_$d")
  end
end

"""
    writevtk(model::DiscreteModel,filebase)
    writevtk(model::DiscreteModel, labels::FaceLabeling, filebase)
"""
function writevtk(model::DiscreteModel,filebase)
  labels = get_face_labeling(model)
  writevtk(model,labels,filebase)
end

function writevtk(model::DiscreteModel, labels::FaceLabeling, filebase)
  for d in 0:num_cell_dims(model)
    grid = Grid(ReferenceFE{d},model)
    cdat = _prepare_cdata_model(labels,d)
    write_vtk_file(grid,"$(filebase)_$d";celldata=cdat)
  end
end

function _prepare_cdata_model(labels,d)
  dface_to_entity = get_face_entity(labels,d)
  cdat = []
  for tag in 1:num_tags(labels)
    dface_to_isontag = zeros(Int,num_faces(labels,d))
    for entity in get_tag_entities(labels,tag)
      _set_entity!(dface_to_isontag,dface_to_entity,entity,tag)
    end
    name = get_tag_name(labels,tag)
    push!(cdat, name => dface_to_isontag )
  end
  push!(cdat,"entity" => dface_to_entity)
  cdat
end

function _set_entity!(dface_to_isontag,dface_to_entity,entity,tag)
  for i in 1:length(dface_to_entity)
    if dface_to_entity[i] == entity
      dface_to_isontag[i] = tag
    end
  end
end

"""

    write_vtk_file(
      trian::Grid,
      filebase;
      celldata=Dict(),
      nodaldata=Dict())

Low level entry point to vtk. Other vtk-related routines in Gridap eventually call this one.

"""
function write_vtk_file(
  trian::Grid, filebase; celldata=Dict(), nodaldata=Dict())
  vtkfile = create_vtk_file(trian, filebase, celldata=celldata, nodaldata=nodaldata)
  outfiles = vtk_save(vtkfile)
end

"""

    create_vtk_file(
      trian::Grid,
      filebase;
      celldata=Dict(),
      nodaldata=Dict())

Low level entry point to vtk. Other vtk-related routines in Gridap eventually call this one.
This function only creates the vtkFile, without writing to disk.

"""
function create_vtk_file(
  trian::Grid, filebase; celldata=Dict(), nodaldata=Dict())

  points = _vtkpoints(trian)
  cells = _vtkcells(trian)
  vtkfile = vtk_grid(filebase, points, cells, compress=false)

  for (k,v) in celldata
    vtk_cell_data(vtkfile, _prepare_data(v), k)
  end
  for (k,v) in nodaldata
    vtk_point_data(vtkfile, _prepare_data(v), k)
  end

  return vtkfile
end

function _vtkpoints(trian)
  D = num_point_dims(trian)
  x = get_node_coordinates(trian)
  xflat = collect(x)
  reshape(reinterpret(Float64,xflat),(D,length(x)))
end

function _vtkcells(trian)

  type_to_reffe = get_reffes(trian)
  cell_to_type = get_cell_type(trian)
  type_to_vtkid = map(get_vtkid, type_to_reffe)
  type_to_vtknodes = map(get_vtknodes, type_to_reffe)
  cell_to_nodes = get_cell_nodes(trian)
  cache = array_cache(cell_to_nodes)

  _generate_vtk_cells(
    cache,
    cell_to_nodes,
    cell_to_type,
    type_to_vtkid,
    type_to_vtknodes)
end

function _generate_vtk_cells(
  cache,
  cell_to_nodes,
  cell_to_type,
  type_to_vtkid,
  type_to_vtknodes)

  V = eltype(cell_to_nodes)
  meshcells = MeshCell{WriteVTK.VTKCellTypes.VTKCellType,V}[]

  d = _vtkcelltypedict()

  cells = 1:length(cell_to_type)
  for cell in cells

    t = cell_to_type[cell]
    vtkid = type_to_vtkid[t]
    vtknodes = type_to_vtknodes[t]

    nodes = getindex!(cache,cell_to_nodes,cell)
    meshcell = MeshCell(d[vtkid], nodes[vtknodes])

    push!(meshcells,meshcell)

  end

  meshcells

end

_prepare_data(v) = v

function _prepare_data(v::AbstractArray{<:MultiValue})
  a = collect(v)
  reinterpret(a)
end

function _prepare_data(v::AbstractArray{<:VectorValue{2,T}}) where T
  a = collect(v)
  b = reshape(reinterpret(T,a),(2,length(a)))
  z = zeros((1,size(b,2)))
  vcat(b,z)
end

"""
"""
function get_vtkid(reffe::LagrangianRefFE)
  basis = get_prebasis(reffe)
  p = get_polytope(reffe)
  get_vtkid(p,basis)
end

"""
"""
function get_vtknodes(reffe::LagrangianRefFE)
  basis = get_prebasis(reffe)
  p = get_polytope(reffe)
  get_vtknodes(p,basis)
end

function get_vtkid(p::ExtrusionPolytope, basis::MonomialBasis)
  exponents = get_exponents(basis)
  vtkid, _ = _vtkinfo_extrusion_polytope(p,exponents)
  vtkid
end

function get_vtknodes(p::ExtrusionPolytope, basis::MonomialBasis)
  exponents = get_exponents(basis)
  _, vtknodes = _vtkinfo_extrusion_polytope(p,exponents)
  vtknodes
end

function get_vtkid(p::SerendipityPolytope,basis::MonomialBasis)
  get_vtkid(p.hex,basis)
end

function get_vtknodes(p::SerendipityPolytope,basis::MonomialBasis)
  get_vtknodes(p.hex,basis)
end

function _vtkinfo_extrusion_polytope(p,exponents)

  # Taken from the vtk specification
  # https://vtk.org/wp-content/uploads/2015/04/file-formats.pdf

  n_nodes = length(exponents)

  if p == VERTEX
    if n_nodes == 1
      vtkid = 1
      vtknodes = [1,]
    else
      @notimplemented
    end

  elseif p == SEGMENT
    if n_nodes == 2
      vtkid = 3
      vtknodes = [1,2]
    elseif n_nodes == 3
      vtkid = 21
      vtknodes = [1,2,3]
    else
      @notimplemented
    end

  elseif p == TRI
    if n_nodes == 3
      vtkid = 5
      vtknodes = [1,2,3]
    elseif n_nodes == 6
      vtkid = 22
      vtknodes = [1,2,3,4,6,5]
    elseif n_nodes == 10
      vtkid = 69
      vtknodes = [1,2,3,4,5,8,9,7,6,10]
    elseif n_nodes == 15
      vtkid = 69
      vtknodes = [1,2,3,4,5,6,10,11,12,9,8,7,13,14,15]
    elseif n_nodes == 21
      vtkid = 69
      vtknodes = [1,2,3,4,5,6,7,12,13,14,15,11,10,9,8,16,18,21,17,20,19]
    else
      @notimplemented
    end

  elseif p == QUAD
    if n_nodes == 4
      vtkid = 9
      vtknodes = [1,2,4,3]
    elseif n_nodes == 8
      vtkid = 23
      vtknodes = [1,2,4,3,5,8,6,7]
    elseif n_nodes == 9
      vtkid = 28
      vtknodes = [1,2,4,3,5,8,6,7,9]
    elseif n_nodes == 16
      vtkid = 70
      vtknodes = [1,2,4,3,5,6,11,12,7,8,9,10,13,14,15,16]
    elseif n_nodes == 25
      vtkid = 70
      vtknodes = [1,2,4,3,5,6,7,14,15,16,8,9,10,11,12,13,17,18,19,20,21,22,23,24,25]
    else
      @notimplemented
    end

  elseif p == TET
    if n_nodes == 4
      vtkid = 10
      vtknodes = [1,2,3,4]
    elseif n_nodes == 10
      vtkid = 24
                 #0,1,2,3,4,5,6,7,8,9
      vtknodes = [1,2,3,4,5,7,6,8,9,10]
    elseif n_nodes == 20
      vtkid = 71
      vtknodes = [1,2,3,4,5,6,9,10,7,8,11,12,13,14,15,16,18,20,19,17]
    else
      @notimplemented
    end

  elseif p == HEX
    if n_nodes == 8
      vtkid = 12
      vtknodes = [1,2,4,3,5,6,8,7]
    elseif n_nodes == 27
      vtkid = 72
      vtknodes = [1,2,4,3,5,6,8,7,9,14,10,13,11,16,12,15,17,18,19,20,25,26,23,24,21,22,27]
    elseif n_nodes == 64
      vtkid = 72
      vtknodes = [1,2,4,3,5,6,8,7,9,10,19,20,11,12,17,18,13,14,23,24,15,16,21,22,25,26,27,
                  28,29,30,31,32,49,50,51,52,53,54,55,56,41,42,43,44,45,46,47,48,33,34,35,
                  36,37,38,39,40,57,58,59,60,61,62,63,64]
    else
      @notimplemented
    end

  elseif p == WEDGE
    if n_nodes == 6
      vtkid = 13
      vtknodes = [1,3,2,4,6,5]
    else
      @notimplemented
    end

  elseif p == PYRAMID
    if n_nodes == 5
      vtkid = 14
      vtknodes = [1,2,4,3,5]
    else
      @notimplemented
    end

  else
    @notimplemented "vtkid not implemented for given ExtrusionPolytope"
  end

  (vtkid, vtknodes)
end

function _vtkcelltypedict()
  d = Dict{Int,WriteVTK.VTKCellTypes.VTKCellType}()
  d[VTK_VERTEX.vtk_id] = VTK_VERTEX
  d[VTK_LINE.vtk_id] = VTK_LINE
  d[VTK_LINE.vtk_id] = VTK_LINE
  d[VTK_TRIANGLE.vtk_id] = VTK_TRIANGLE
  d[VTK_QUAD.vtk_id] = VTK_QUAD
  d[VTK_TETRA.vtk_id] = VTK_TETRA
  d[VTK_HEXAHEDRON.vtk_id] = VTK_HEXAHEDRON
  d[VTK_WEDGE.vtk_id] = VTK_WEDGE
  d[VTK_QUADRATIC_QUAD.vtk_id] = VTK_QUADRATIC_QUAD
  d[VTK_BIQUADRATIC_QUAD.vtk_id] = VTK_BIQUADRATIC_QUAD
  d[VTK_QUADRATIC_TRIANGLE.vtk_id] = VTK_QUADRATIC_TRIANGLE
  d[VTK_QUADRATIC_TETRA.vtk_id] = VTK_QUADRATIC_TETRA
  d[VTK_QUADRATIC_EDGE.vtk_id] = VTK_QUADRATIC_EDGE
  d[VTK_QUADRATIC_HEXAHEDRON.vtk_id] = VTK_QUADRATIC_HEXAHEDRON
  d[VTK_PYRAMID.vtk_id] = VTK_PYRAMID
  d[VTK_LAGRANGE_TRIANGLE.vtk_id] = VTK_LAGRANGE_TRIANGLE
  d[VTK_LAGRANGE_QUADRILATERAL.vtk_id] = VTK_LAGRANGE_QUADRILATERAL
  d[VTK_LAGRANGE_TETRAHEDRON.vtk_id] = VTK_LAGRANGE_TETRAHEDRON
  d[VTK_LAGRANGE_HEXAHEDRON.vtk_id] = VTK_LAGRANGE_HEXAHEDRON
  #d[VTK_BIQUADRATIC_HEXAHEDRON.vtk_id] = VTK_BIQUADRATIC_HEXAHEDRON
  d
end

"""
"""
function writevtk(trian::Triangulation, filebase; order=-1, nsubcells=-1, celldata=Dict(), cellfields=Dict())
  vtkfile = createvtk(trian,filebase, order=order, nsubcells=nsubcells, celldata=celldata, cellfields=cellfields)
  outfiles = vtk_save(vtkfile)
end

"""
"""
function createvtk(trian::Triangulation, filebase; order=-1, nsubcells=-1, celldata=Dict(), cellfields=Dict())
  visdata = visualization_data(trian; order=order,
    nsubcells=nsubcells, celldata=celldata, cellfields=cellfields)
  create_vtk_file(visdata.grid,filebase,celldata=visdata.celldata,nodaldata=visdata.nodaldata)
end

"""
    writevtk(
      cell_to_points::AbstractArray{<:AbstractArray{<:Point}},
      filename;
      celldata=Dict(),
      nodaldata=Dict())
"""
function writevtk(
  cell_to_points::AbstractArray{<:AbstractArray{<:Point}},
  filename; celldata=Dict(), nodaldata=Dict())
  vtkfile = createvtk(cell_to_points, filename, celldata=celldata, nodaldata=nodaldata)
  outfiles = vtk_save(vtkfile)
end

"""
    createvtk(
      cell_to_points::AbstractArray{<:AbstractArray{<:Point}},
      filename;
      celldata=Dict(),
      nodaldata=Dict())
"""
function createvtk(
  cell_to_points::AbstractArray{<:AbstractArray{<:Point}},
  filename; celldata=Dict(), nodaldata=Dict())

  node_to_point, cell_to_offset = _prepare_node_to_coords(cell_to_points)
  nnodes = length(node_to_point)
  cell_to_nodes = identity_table(Int,Int32,nnodes)
  cell_to_ctype = fill(Int8(1),nnodes)
  ctype_to_reffe = [VERTEX1,]

  node_to_cell = _prepare_node_to_cell(cell_to_offset,nnodes)

  grid = UnstructuredGrid(
    node_to_point,
    cell_to_nodes,
    ctype_to_reffe,
    cell_to_ctype,
    Val{true}())

  cdata = _prepare_cdata(celldata,node_to_cell)
  pdata = _prepare_pdata_for_cell_points(nodaldata)

  create_vtk_file(grid,filename,celldata=cdata,nodaldata=pdata)

end

function _prepare_node_to_cell(cell_to_offset,nnodes)
  node_to_cell = zeros(Int,nnodes)
  ncells = length(cell_to_offset)
  for cell in 1:ncells
    nini = cell_to_offset[cell]+1
    if cell < ncells
      nend = cell_to_offset[cell+1]
    else
      nend = nnodes
    end
    for node in nini:nend
      node_to_cell[node] = cell
    end
  end
  node_to_cell
end

function _prepare_pdata_for_cell_points(nodaldata)
  pdata = Dict()
  for (k,v) in nodaldata
    pdata[k], = _prepare_node_to_coords(v)
  end
  pdata
end


text:
using ImagineInterface
using Test, IntervalSets, Unitful

import ImagineInterface: check_max_speed, ValidationState, check_piezo, check_piezos
import ImagineInterface: check_pulse_duration, check_interpulse_duration, check_pulse_interval
import ImagineInterface: check_camera, check_cameras, check_laser, check_lasers

#check_max_speed
a = ones(Int16, 10)
a[3] = 3
@test check_max_speed(a, 5, 5) == true
@test check_max_speed(a, 4, 5) == true
@test check_max_speed(a, 3, 5) == false
a[3] = 4
@test check_max_speed(a, 5, 5) == false
a[3] = 1
a[5] = 4
@test check_max_speed(a, 5, 5) == false
a[5] = 1
a[10] = 4
@test check_max_speed(a, 5, 5) == true
a[10] = 1
a[9] = 4
@test check_max_speed(a, 5, 5) == false

#check_piezo
rig = "ocpi-2"
ocpi2 = rigtemplate(rig; sample_rate = 1000*inv(Unitful.s))
pos = getpositioners(ocpi2)[1]
b = zeros(rawtype(pos), 2000)
append!(pos, "b", b)

speed_lim = ImagineInterface.PIEZO_MAX_SPEED[rig]
win_sz = 10
thresh_um = (win_sz/samprate(pos)) * speed_lim
#thresh_raw is the allowable cumulative raw change in piezo signal within the sliding window
thresh_raw = ceil(Int, (thresh_um/width(interval_world(pos))) * width(interval_raw(pos)))

@test isa(check_piezo(pos; window_sz = win_sz), ValidationState)

b[1] = round(eltype(b), thresh_raw + 1)
replace!(pos, "b", b)
@test_throws Exception check_piezo(pos; window_sz = win_sz)

b[1] = 0
b[end] = round(eltype(b), thresh_raw + 1)
replace!(pos, "b", b)
@test_throws Exception check_piezo(pos; window_sz = win_sz)

b[end] = 0
b[1001] = round(eltype(b), ceil(Int, thresh_raw/2) + 1)
replace!(pos, "b", b)
@test_throws Exception check_piezo(pos; window_sz = win_sz)

b[1001] = round(eltype(b), ceil(Int, thresh_raw/2) -2)
replace!(pos, "b", b)
@test isa(check_piezo(pos; window_sz = win_sz), ValidationState)

#check_piezos
@test isa(check_piezos(getpositioners(ocpi2); window_sz = win_sz), ValidationState)

#test pulse validation code camera
srate = 1000000*inv(Unitful.s)
rig = "ocpi-2"
ocpirig = rigtemplate(rig; sample_rate = srate)
function test_seq1(first_i, samps_on, start2start)
    start_is = Int[first_i; first_i+start2start; first_i + start2start*2; first_i + start2start*3]
    stop_is = Int[first_i + samps_on-1; first_i+start2start+samps_on-1; first_i+start2start*2+samps_on-1; first_i+start2start*3+samps_on-1]
    return start_is, stop_is
end

function test_seq2(first_i, samps_off, start2start)
    samps_on = start2start - samps_off
    test_seq1(first_i, samps_on, start2start)
end

function gen_samp_seq(start_is, stop_is, tot_length)
    @assert length(start_is) == length(stop_is)
    seq = falses(tot_length)
    for i = 1:length(start_is)
        @assert start_is[i] <= stop_is[i]
        seq[start_is[i]:stop_is[i]] .= true
    end
    return seq
end

sigs = getcameras(ocpirig)
ton = ImagineInterface.CAMERA_ON_TIME[rig]
toff = ImagineInterface.CAMERA_OFF_TIME[rig]
samps_on_tol = ceil(Int, uconvert(inv(unit(srate)), ton) * srate)
samps_off_tol = ceil(Int, uconvert(inv(unit(srate)), toff) * srate)
frate = max_framerate(rig, chip_size(rig)...)
start2start_tol = ceil(Int, uconvert(inv(unit(srate)), 1/frate) * srate)
first_i = 2

start_is, stop_is = test_seq1(first_i, samps_on_tol, start2start_tol)
@test check_pulse_duration(start_is, stop_is, samps_on_tol, samps_off_tol) == 0
@test check_interpulse_duration(start_is, stop_is, samps_on_tol, samps_off_tol) == 0
@test check_pulse_interval(start_is, start2start_tol) == 0
samps = gen_samp_seq(start_is, stop_is, stop_is[end]+1)
append!(sigs[1], "test", samps)
check_camera(sigs[1])
check_cameras(sigs)
append!(sigs[2], "test2", samps)
check_cameras(sigs)

#check_pulse_duration fails because on time is too short
start_is, stop_is = test_seq1(first_i, samps_on_tol-1, start2start_tol)
@test check_pulse_duration(start_is, stop_is, samps_on_tol, samps_off_tol) > 0
@test check_interpulse_duration(start_is, stop_is, samps_on_tol, samps_off_tol) == 0

samps = gen_samp_seq(start_is, stop_is, stop_is[end]+1)
replace!(sigs[1], "test", samps)
@test_throws Exception check_camera(sigs[1])

#check_pulse_interval fails because start2start is to short
start_is, stop_is = test_seq1(first_i, samps_on_tol, start2start_tol-1)
@test check_pulse_duration(start_is, stop_is, samps_on_tol, samps_off_tol) == 0
@test check_interpulse_duration(start_is, stop_is, samps_on_tol, samps_off_tol) == 0

@test check_pulse_interval(start_is, start2start_tol) > 0
samps = gen_samp_seq(start_is, stop_is, stop_is[end]+1)
replace!(sigs[1], "test", samps)
@test_throws Exception check_camera(sigs[1])

start_is, stop_is = test_seq2(first_i, samps_off_tol, start2start_tol)
@test check_pulse_duration(start_is, stop_is, samps_on_tol, samps_off_tol) == 0
@test check_interpulse_duration(start_is, stop_is, samps_on_tol, samps_off_tol) == 0

@test check_pulse_interval(start_is, start2start_tol) == 0
samps = gen_samp_seq(start_is, stop_is, stop_is[end]+1)
replace!(sigs[1], "test", samps)
check_camera(sigs[1])

#check_interpulse_duration fails because off time is too short
start_is, stop_is = test_seq2(first_i, samps_off_tol-1, start2start_tol)
@test check_pulse_duration(start_is, stop_is, samps_on_tol, samps_off_tol) == 0
@test check_interpulse_duration(start_is, stop_is, samps_on_tol, samps_off_tol) > 0

@test check_pulse_interval(start_is, start2start_tol) == 0
samps = gen_samp_seq(start_is, stop_is, stop_is[end]+1)
replace!(sigs[1], "test", samps)
@test_throws Exception check_camera(sigs[1])

#check_pulse_interval fails because start2start is to short
start_is, stop_is = test_seq2(first_i, samps_off_tol, start2start_tol-1)
@test check_pulse_duration(start_is, stop_is, samps_on_tol, samps_off_tol) == 0
@test check_interpulse_duration(start_is, stop_is, samps_on_tol, samps_off_tol) == 0

@test check_pulse_interval(start_is, start2start_tol) > 0
samps = gen_samp_seq(start_is, stop_is, stop_is[end]+1)
replace!(sigs[1], "test", samps)
@test_throws Exception check_camera(sigs[1])

############################ Similar tests for lasers (should think about simplifying this)
sigs = getlasers(ocpirig)
ton = ImagineInterface.LASER_ON_TIME[rig]
toff = ImagineInterface.LASER_OFF_TIME[rig]
samps_on_tol = ceil(Int, uconvert(inv(unit(srate)), ton) * srate)
samps_off_tol = ceil(Int, uconvert(inv(unit(srate)), toff) * srate)
start2start_tol = samps_on_tol+samps_on_tol+2 #currently we don't limit laser pulse rate (except indirectly through ton and toff) so just setting this to something safe
first_i = 2

start_is, stop_is = test_seq1(first_i, samps_on_tol, start2start_tol)
@test check_pulse_duration(start_is, stop_is, samps_on_tol, samps_off_tol) == 0
@test check_interpulse_duration(start_is, stop_is, samps_on_tol, samps_off_tol) == 0

@test check_pulse_interval(start_is, start2start_tol) == 0
samps = gen_samp_seq(start_is, stop_is, stop_is[end]+1)
append!(sigs[1], "test3", samps)
check_laser(sigs[1])
check_lasers(sigs)
append!(sigs[2], "test4", samps)
check_lasers(sigs)
append!(sigs[3], "test5", trues(length(samps))) #laser always on

#check_pulse_duration fails because on time is too short
start_is, stop_is = test_seq1(first_i, samps_on_tol-1, start2start_tol)
@test check_pulse_duration(start_is, stop_is, samps_on_tol, samps_off_tol) > 0
@test check_interpulse_duration(start_is, stop_is, samps_on_tol, samps_off_tol) == 0

@test check_pulse_interval(start_is, start2start_tol) == 0
samps = gen_samp_seq(start_is, stop_is, stop_is[end]+1)
replace!(sigs[1], "test3", samps)
@test_throws Exception check_laser(sigs[1])

start_is, stop_is = test_seq2(first_i, samps_off_tol, start2start_tol)
@test check_pulse_duration(start_is, stop_is, samps_on_tol, samps_off_tol) == 0
@test check_interpulse_duration(start_is, stop_is, samps_on_tol, samps_off_tol) == 0

@test check_pulse_interval(start_is, start2start_tol) == 0
samps = gen_samp_seq(start_is, stop_is, stop_is[end]+1)
replace!(sigs[1], "test3", samps)
check_laser(sigs[1])

#check_interpulse_duration fails because off time is too short
start_is, stop_is = test_seq2(first_i, samps_off_tol-1, start2start_tol)
@test check_pulse_duration(start_is, stop_is, samps_on_tol, samps_off_tol) == 0
@test check_interpulse_duration(start_is, stop_is, samps_on_tol, samps_off_tol) > 0

@test check_pulse_interval(start_is, start2start_tol) == 0
samps = gen_samp_seq(start_is, stop_is, stop_is[end]+1)
replace!(sigs[1], "test3", samps)
@test_throws Exception check_laser(sigs[1])

@testset "Reduced frame size" begin
    μm, s = Unitful.μm, Unitful.s
    sample_rate = 50000s^-1 #analog output samples per second
    rig = "realm"
    realm = rigtemplate(rig; sample_rate = sample_rate)
    positioners = getpositioners(realm)
    pos = positioners[1]
    ############STACK PARAMETERS#################
    pmin = 0.0*μm #Piezo start position in microns
    pmax = 200.0*μm #Piezo stop position in microns
    stack_img_time = 0.26s #Time to complete the imaging sweep with the piezo (remember we also need to reset it to its starting position)
    reset_time = 0.001s #Time to reset piezo to starting position.  This time plus "stack_img_time" determines how long it takes to complete an entire stack and be ready to start a new stack
    z_spacing = 3.1μm #The space between slices in the z-stack.
    z_pad = 1.0μm #Set this greater than 0 if you want to ignore the edges of the positioner sweep (only take slices in a central region)

    hmax, vmax = chip_size(rig)
    h = 1000 #doesn't matter with current cameras
    v = 400
    mx_f = max_framerate(rig, h,v)
    mn_exp = 1/mx_f #This is the minimum possible exposure time
    mx_f = max_framerate(rig, hmax, vmax)
    mx_exp = 1/mx_f
    exp_time = 0.002s  #Exposure time of the camera. Make sure this is greater than mn_exp and less than mx_exp above
    flash_frac = 1.1 #fraction of time to keep laser on during exposure.  If you set this greater than 1 then the laser will stay on constantly during the imaging sweep

    d = gen_bidirectional_stack(pmin, pmax, z_spacing, stack_img_time, exp_time, sample_rate, flash_frac; z_pad = z_pad, rig = rig)
    nframes = d["nframes"] #note that there are twice as many frames as in the unidirectional example because one cycle of the positioner includes two stacks

    #Let's append our newly-created sample vectors to their respective commands in the template
    append!(pos, "bi_stack_pos", d["positioner"])
    lasers = getlasers(realm)
    las1 = lasers[1] #You could also append to any of the other lasers
    append!(las1, "bi_stack_las1", d["laser"])
    cams = getcameras(realm)
    cam1 = cams[1] #You could also append to the other camera
    append!(cam1, "bi_stack_cam1", d["camera"])

    nstacks = 1
    fn = tempname() * ".json"
    write_commands(fn, realm, nstacks, nframes, exp_time; isbidi = true, vertical_lines=v)
    rm(fn)
end

max_forks_repo_name: HolyLab/ImagineInterface


text:
#################### Model Initialization ####################

function setinits!(m::Model, inits::Dict{Symbol})
  m.hasinputs || throw(ArgumentError("inputs must be set before inits"))
  m.iter = 0
  for key in keys(m, :dependent)
    node = m[key]
    if isa(node, AbstractStochastic)
      haskey(inits, key) ||
        throw(ArgumentError("missing initial value for node : $key"))
      setinits!(node, m, inits[key])
    else
      setinits!(node, m)
    end
  end
  m.hasinits = true
  m
end

function setinits!(m::Model, inits::Vector{V} where V<:Dict{Symbol})
  n = length(inits)
  m.states = Array{ModelState}(undef, n)
  for i in n:-1:1
    setinits!(m, inits[i])
    m.states[i] = ModelState(unlist(m), deepcopy(gettune(m)))
  end
  m
end

function setinputs!(m::Model, inputs::Dict{Symbol})
  for key in keys(m, :input)
    haskey(inputs, key) ||
      throw(ArgumentError("missing inputs for node : $key"))
    isa(inputs[key], AbstractDependent) &&
      throw(ArgumentError("inputs cannot be Dependent types"))
    m.nodes[key] = deepcopy(inputs[key])
  end
  m.hasinputs = true
  m
end

function setsamplers!(m::Model, samplers::Vector{V} where V<:Sampler)
  m.samplers = deepcopy(samplers)
  for sampler in m.samplers
    sampler.targets = keys(m, :target, sampler.params)
  end
  m
end


text:
"""
    walk(T) -> trajectory

模拟随机行走
"""

function walk(T)
    x = 0
    trajectory = [x] #轨迹, 数组

    for t in 1:T
        if rand()< 0.5
            x += 1
        else
            x -= 1
        end
        push!(trajectory, x)
    end

    trajectory
end

# using GR,Interact
# T = 100
# plot(0:T,walk(T))

# @manioulate for T in 10:100
#     plot(0:T, walk(T))
# end

max_stars_repo_name: liushooter/OneDayOneCommit

max_issues_repo_name: liushooter/OneDayOneCommit

max_forks_repo_name: liushooter/OneDayOneCommit


text:
using StatsBase, JuMP, GLPK, Test

function optimal_price_probabilities(
    price, demand, capacity; verbose = false)
    revenue = price .* demand
    model = Model(with_optimizer(GLPK.Optimizer))
    @variable(model, x[1:length(price)], lower_bound=0)
    # Objective: maximize profit
    @objective(model, Max, revenue' * x)
    # Constraint: can carry all
    @constraint(model, sum(x) <= 1)
    #=
    for i in 1:length(prices)
        @constraint(model, x[i] <= demand[i])
    end
    =#
    # Solve problem using MIP solver
    JuMP.optimize!(model)
    if verbose
        println("Objective is: ", JuMP.objective_value(model))
        println("Solution is:")
        for i in 1:length(price)
            print("x[$i] = ", JuMP.value(x[i]))
            println(", revenue = ", price[i] * JuMP.value(x[i]))
        end
    end
    @test JuMP.termination_status(model) == MOI.OPTIMAL
    @test JuMP.primal_status(model) == MOI.FEASIBLE_POINT
    Weights(JuMP.value.(x))
end


text:
using CvxCompress, Blosc, LinearAlgebra

Blosc.set_num_threads(4)

nz, nx = 2048, 4096
bz, bx = 32, 32
scl = 1e0

x = rand(Float32, nz, nx)

# CvxCompress
y = zeros(UInt32, nz*nx)
clength_cvx = CvxCompress.compress!(y, CvxCompressor((bz,bx),scl), x)
t_cvx_compress = @elapsed clength_cvx = CvxCompress.compress!(y, CvxCompressor((bz,bx),scl), x)
xx = similar(x)
CvxCompress.decompress!(xx, CvxCompressor((bz,bx),scl), y, clength_cvx)
t_cvx_decompress = @elapsed CvxCompress.decompress!(xx, CvxCompressor((bz,bx),scl), y, clength_cvx)
snr_cvx = 10*log10(norm(x)^2 / norm(x-xx)^2)

# Blosc (no quantization)
y = zeros(UInt8, nz*nx*4+512)
clength_blosc = Blosc.compress!(y, x)
t_blosc_compress = @elapsed clength_blosc = Blosc.compress!(y, x)
xx=zeros(Float32,nz*nx)
Blosc.decompress!(xx, y)
t_blosc_decompress = @elapsed Blosc.decompress!(xx, y)
xx = reshape(xx,nz,nx)
snr_blosc = 10*log10(norm(x)^2 / (norm(x-xx)^2))

# Blosc (quantization)
function blosc_quant(nz,nx,x,y)
    xx = zeros(UInt16, nz, nx)
    mn = minimum(x)
    d = maximum(x) - mn + eps(Float32)
    sc = d > eps(Float32) ? typemax(UInt16)/d : 1.0f0
    for i = 1:length(x)
        @inbounds xx[i] = unsafe_trunc(UInt16, (x[i]-mn)*sc)
    end
    clength_blosc_quant = Blosc.compress!(y, xx)
    clength_blosc_quant, mn, sc
end

function blosc_dequant(nz,nx,x,y,mn,sc)
    xx = zeros(UInt16, nz*nx)
    Blosc.decompress!(xx,y)
    for i = 1:length(x)
        x[i] = xx[i]/sc + mn
    end
end

clength_blosc_quant, mn, sc = blosc_quant(nz,nx,x,y)
t_bloscquant_compress = @elapsed blosc_quant(nz,nx,x,y)
blosc_dequant(nz,nx,x,y,mn,sc)
t_bloscquant_decompress = @elapsed blosc_dequant(nz,nx,x,y,mn,sc)
snr_bloscquant = 10*log10(norm(x)^2 / (norm(x-xx)^2))

#readme = open("README.md", "w")
readme = stdout
for io in (stdout, readme)
    write(io, "# times (compress, de-compress):\n")
    write(io, "* cvx:         $(t_cvx_compress),$(t_cvx_decompress) seconds\n")
    write(io, "* blosc:       $(t_blosc_compress),$(t_blosc_decompress) seconds\n")
    write(io, "* blosc-quant: $(t_bloscquant_compress),$(t_bloscquant_decompress) seconds\n")

    write(io, "\n")

    write(io, "# compression:\n")
    write(io, "* cvx:         $((clength_cvx*4)/(nz*nx*4))\n")
    write(io, "* blosc:       $((clength_blosc)/(nz*nx*4))\n")
    write(io, "* blosc-quant: $((clength_blosc_quant)/(nz*nx*4))\n")

    write(io, "\n")

    write(io, "# signal-to-noise:\n")
    write(io, "* cvx:         $(snr_cvx)\n")
    write(io, "* blosc:       $(snr_blosc)\n")
    write(io, "* blosc-quant: $(snr_bloscquant)\n")
end
#close(readme)