text:
# Example 6.3 of
# R. Jungers, A. Cicone and N. Guglielmi,
# "Lifted polytope methods for computing the joint spectral radius."
# SIAM Journal on Matrix Analysis and Applications, SIAM, 2014, 35, 391-410.
# The JSR is 1.618033988749895

@testset "[JCG14] Example 6.3" begin
    expected_lb = [1.1441, 1.3606]
    expected_ub = [1.6181, 1.6181]
    s = discreteswitchedsystem([[ 0  1  1;
                                  1  0  0;
                                  0 -1  0],
                                [ 0  1  0;
                                 -1  0  1;
                                 -1  0  0]])
    smp = periodicswitching(s, [1, 2])
    for factory in sdp_factories
        sosdata(s).lb = 0
        tol = 1e-4
        for d in 1:2
            lb, ub = soslyapb(s, d, factory=factory, tol=tol)
            @test isapprox(lb, expected_lb[d], rtol=tol)
            @test isapprox(ub, expected_ub[d], rtol=tol)
            seq = sosbuildsequence(s, d, p_0=:Primal)
            psw = findsmp(seq)
            @test psw !== nothing
            @test psw == smp
        end
    end
end

ext: jl


text:
"""
    Valuation

A data structure for approximating the valuation of a solution path ``x(t)``.
"""
mutable struct Valuation
    val_x::Vector{Float64}
    val_tẋ::Vector{Float64}
    Δval_x::Vector{Float64}
    Δval_tẋ::Vector{Float64}

    # Data for computing the val_x and val_tẋ and ν̈ by finite differences
    val_x_data::NTuple{2,Vector{Float64}}
    val_ẋ_data::NTuple{2,Vector{Float64}}
    logx_data::NTuple{2,Vector{Float64}}
    logẋ_data::NTuple{2,Vector{Float64}}
    logt_data::NTuple{2,Float64}
end

Valuation(x::AbstractVector) = Valuation(length(x))
Valuation(n::Integer) =
    Valuation((zeros(n) for i = 1:4)..., ((zeros(n), zeros(n)) for i = 1:4)..., (NaN, NaN))

function Base.show(io::IO, val::Valuation)
    print(io, "Valuation :")
    for field in [:val_x, :val_tẋ, :Δval_x, :Δval_tẋ]
        vs = [Printf.@sprintf("%#.4g", v) for v in getfield(val, field)]
        print(io, "\n • ", field, " → ", "[", join(vs, ", "), "]")
    end
end
Base.show(io::IO, ::MIME"application/prs.juno.inline", v::Valuation) = v

"""
    init!(val::Valuation)

Initialize the valuation estimator to start from scratch.
"""
function init!(val::Valuation)
    val.val_x .= 0.0
    val.val_tẋ .= 0.0
    val.Δval_x .= 0.0
    val.Δval_tẋ .= 0.0
    val.val_x_data[1] .= val.val_x_data[2] .= 0.0
    val.val_ẋ_data[1] .= val.val_ẋ_data[2] .= 0.0
    val.logx_data[1] .= val.logx_data[2] .= 0.0
    val.logẋ_data[1] .= val.logẋ_data[2] .= 0.0
    val.logt_data = (NaN, NaN)
    val
end

function ν(x, ẋ, t)
    u, v = reim(x)
    u¹, v¹ = reim(ẋ)
    xx = abs2(x)

    μ = u * u¹ + v * v¹
    l = μ / xx

    t * l
end

function ν_ν¹(x, ẋ, x², t)
    u, v = reim(x)
    u¹, v¹ = reim(ẋ)
    u², v² = reim(x²)

    xx = abs2(x)

    μ = u * u¹ + v * v¹
    l = μ / xx

    μ¹ = u * u² + u¹^2 + v * v² + v¹^2
    l¹ = μ¹ / xx - 2 * l^2

    t * l, t * l¹ + l
end

logabs(x) = log(fast_abs(x))

function update!(val::Valuation, pred::Predictor, t::Real)
    logt = log(t)
    ν₂, ν₁ = val.val_x_data
    logx₂, logx₁, = val.logx_data
    logẋ₂, logẋ₁, = val.logẋ_data
    logt₂, logt₁ = val.logt_data
    val_ẋ₂, val_ẋ₁ = val.val_ẋ_data

    # We use finite difference when the predictor is in hermite mode since then we only
    # have ẋ available
    diff = pred.winding_number > 1 && !isnan(logt₂)
    for (i, (xᵢ, ẋᵢ, x²ᵢ, x³ᵢ)) in enumerate(pred.tx³)
        logxᵢ = logabs(xᵢ)
        logẋᵢ = logabs(ẋᵢ)

        if diff
            νᵢ = ν(xᵢ, ẋᵢ, t)
            Δνᵢ = finite_diff(νᵢ, logt, ν₂[i], logt₂, ν₁[i], logt₁)
            val.val_x[i] = νᵢ
            val.Δval_x[i] = Δνᵢ

            val_ẋᵢ = finite_diff(logẋᵢ, logt, logẋ₂[i], logt₂, logẋ₁[i], logt₁)
            Δval_ẋᵢ = finite_diff(val_ẋᵢ, logt, val_ẋ₂[i], logt₂, val_ẋ₁[i], logt₁)
            val.val_tẋ[i] = val_ẋᵢ + 1
            val.Δval_tẋ[i] = Δval_ẋᵢ
        else
            νᵢ, ν¹ = ν_ν¹(xᵢ, ẋᵢ, 2x²ᵢ, t)
            val.val_x[i] = νᵢ
            val.Δval_x[i] = t * ν¹
            ν₂[i], ν₁[i] = νᵢ, ν₂[i]

            val_ẋᵢ, ν¹ = ν_ν¹(ẋᵢ, 2x²ᵢ, 6x³ᵢ, t)
            val.val_tẋ[i] = val_ẋᵢ + 1
            val.Δval_tẋ[i] = t * ν¹
        end
        logx₂[i], logx₁[i] = logxᵢ, logx₂[i]
        logẋ₂[i], logẋ₁[i] = logẋᵢ, logẋ₂[i]
    end

    val.logt_data = (logt, logt₂)

    val
end

"""
    finite_diff(ν₃, s₃, ν₂, s₂, ν₁, s₁)

Compute the derivative `ν̇` of a function ``ν(s)`` at `s₃` by using a finite
difference scheme. This uses the values `νᵢ = ν(sᵢ)` for ``i=1,…,3``.
Since we have a non-uniform grid, we need a more elaborate difference scheme.
The implementation follows the formulas derived in [^BS05].

[^BS05]: Bowen, M. K., and Ronald Smith. "Derivative formulae and errors for non-uniformly
  spaced points." Proceedings of the Royal Society A: Mathematical, Physical and Engineering
  Sciences 461.2059 (2005): 1975-1997.
"""
function finite_diff(ν, s, ν₂, s₂, ν₁, s₁)
    Δ₁, Δ₂, Δ₁₂ = s - s₁, s - s₂, s₁ - s₂
    (Δ₂ * ν₁) / (Δ₁₂ * Δ₁) - ((Δ₁₂ + Δ₂) * ν₂) / (Δ₁₂ * Δ₂) - (Δ₁₂ * ν) / (Δ₁ * Δ₂)
end

function at_infinity_tol!(
    at_infinity_tols,
    val::Valuation;
    finite_tol::Float64,
    zero_is_finite::Bool,
)
    @unpack val_x, val_tẋ, Δval_x, Δval_tẋ = val

    for (i, val_xᵢ) in enumerate(val_x)
        ε∞ = max(
            abs(1.0 - val_tẋ[i] / val_xᵢ),
            abs(Δval_x[i] / val_xᵢ),
            abs(Δval_tẋ[i] / val_tẋ[i]),
        )
        if isnan(ε∞)
            at_infinity_tols[i] = Inf
            continue
        end

        # ∞ check
        if val_xᵢ + ε∞ < -finite_tol
            at_infinity_tols[i] = ε∞
            # 0 check
        elseif !zero_is_finite && val_xᵢ - ε∞ > -finite_tol
            at_infinity_tols[i] = ε∞
        else
            at_infinity_tols[i] = Inf
        end
    end
    at_infinity_tols
end

function is_finite(
    val::Valuation;
    finite_tol::Float64,
    zero_is_finite::Bool,
    max_winding_number::Int,
)
    @unpack val_x, val_tẋ, Δval_x, Δval_tẋ = val

    δ = inv(max_winding_number)
    for (i, val_xᵢ) in enumerate(val_x)
        # Case a: val(x) = 0
        if abs(val_xᵢ) < finite_tol
            if !(abs(Δval_x[i]) < finite_tol) || val_tẋ[i] < 0.5δ
                return false
            end
            # Case b: val(x) = val(tẋ) > 0
        elseif zero_is_finite && val_xᵢ > (δ - finite_tol)
            ε∞ = max(
                abs(1.0 - val_tẋ[i] / val_xᵢ),
                abs(Δval_x[i] / val_xᵢ),
                abs(Δval_tẋ[i] / val_tẋ[i]),
            )
            if !(ε∞ < finite_tol)
                return false
            end
        else
            return false
        end
    end
    return true
end

function estimate_winding_number(val::Valuation; max_winding_number::Int)
    m = 1
    min_err = Inf
    for k = 1:max_winding_number
        err = check_winding_number(val, k)
        if err < min_err
            m = k
            min_err = err
        end
    end
    m, min_err
end

function check_winding_number(val::Valuation, m::Int)
    err = 0.0
    for vᵢ in val.val_tẋ
        mvᵢ = m * vᵢ
        errᵢ = abs(round(mvᵢ) - mvᵢ)
        err = ifelse(errᵢ > err, errᵢ, err)
    end
    err
end


text:
#utilities for QUtilities
#some of these are candidates for Finometrics
const CRSP_MISSINGS = [
  "A", "B", "C", "D", "E", "S", "T", "P",
  -55.0, -66.0, -77.0, -88.0, -99.0,
  -55, -66, -77, -88, -99,
  "-55.0", "-66.0", "-77.0", "-88.0", "-99.0"]

const CRSP_PARSED_MISSINGS = (-55.0, -66.0, -77.0, -88.0, -99.0)

const FRET = :ret
const CQUART_FVALS = [:mknegcash, :bknegcash, :mkflev, :bkflev, :mkliab, :bkliab]
const CQUART_FCONTROLS = :LLret12m
const CQUART_FWEIGHTONS = nothing
const CQUART_FNAMES = (s::Symbol->Symbol("P_",s)).(CQUART_FVALS)
const CQUART_FGROUPS = (s::Symbol->Symbol("G_",s)).(CQUART_FNAMES)

const CQUART_FCONTROLS_ALT = :Lmkat
const CQUART_FNAMES_ALT = (s::Symbol->Symbol(s,"_ac")).(CQUART_FNAMES)
const CQUART_FGROUPS_ALT = (s::Symbol->Symbol("G_",s)).(CQUART_FNAMES_ALT)

const USED_COMPA_VALUE_FIELDS = [:che, :ch, :at, :dltt, :dlc, :ceq, :txditc,
  :txdb, :pstk, :lt]
const USED_COMPQ_VALUE_FIELDS = [:rdq]
#const USED_COMP_VALUE_FIELDS = [USED_COMPA_VALUE_FIELDS; USED_COMPQ_VALUE_FIELDS;]
const USED_CRSP_FIELDS = [:ret, :price, :shrout]

#final columns for the decompressed file
const RETAINED_COLS_CCM = [:gvkey, :lpermno, :linkeffdate, :linkenddate]

const RETAINED_COLS_COMP = [:gvkey; :fyear; :fyrmonth;
  :linkeffdate; :linkenddate; :Nadate;
  :ddate; :adate; :lpermno; :begindate; :enddate; USED_COMPA_VALUE_FIELDS;
  :bknegcash; :cash; :netdebt; :bkequity; :bkflev; :bkliab; :mkequity; :op;]

#not sure if we really need much of these
const RETAINED_COLS_UNIV = [RETAINED_COLS_COMP;
  :compid; :permno; :date; :ret; :mktcap; :net]

#Only keep what we need, also serves as a check to make sure we ahve what we think
const RETAINED_COLS_DATA_SERIES = unique([:gvkey; :fyear; :fyrmonth;
  :linkeffdate; :linkenddate; :Nadate;
  :ddate; :adate; :lpermno; :begindate; :enddate; :mkequity; :Lmkequity; :mkat;
  CQUART_FCONTROLS; CQUART_FCONTROLS_ALT;
  CQUART_FVALS; (s->Symbol(:G_,s)).(CQUART_FNAMES); (s->Symbol(:G_,s)).(CQUART_FNAMES_ALT);
  :vol252d; :Lvol252d; :LLvol252d;
  :Dat; :LDat; :ret12m; :Lret12m; :LLret12m; :bm; :Lbm; :op; :Lop])

#Table controls
const CONTROL0 = Vector{Symbol}()
const T1_CONTROL4 = [:Lbm; :Lmkequity; :LDat; :Lop]
const T1_CONTROL10 = [T1_CONTROL4; :L12ret12m; :L24ret12m;
  :fyret12m; :L12vol252d; :L24vol252d; :fyvol252d]

#@inline Base.isfinite(::Missing)=missing

#NOTE: The below adjusts for logged vs unlogged returns
#const RETAINED_FACT_COLS = (( RETURN_FIELD == :lret) ?
#  [RETURN_FIELD; FACTOR_FIELDS; :lrfr, :lmkt] :
#  [RETURN_FIELD; FACTOR_FIELDS; :rfr, :mkt]

#consider promiting to the below to finometrics
@inline obliterate!(df::DataFrame)::DataFrame = select!(df, Not(names(df)))

#checks if a column or multiple columns are a key
@inline function checkdfkey(df::AbstractDataFrame, keys::Union{Symbol, Vector{Symbol}})
  N::Int = size(df,1)
  Nunique::Int = length(groupby(df, keys))

  (N==Nunique) || error("$keys is not a key!!")

  return nothing
end

#converts gzip to lz4
function gz2lz4(file::String; filetype::String = "csv", location::String = "",
  gzextension::String = "gz", lz4extension::String = "lz4",
  deleteold::Bool = false)

  infile::String = "$location\\$file.$filetype.$gzextension"
  outfile::String = "$location\\$file.$filetype.$lz4extension"

  @assert isfile(infile)

  #Below doesnt work
  indata = open(infile)
  outdata = open(outfile, "w")
  s = GzipDecompressorStream(LZ4CompressorStream(outdata))
  write(s, indata)
  close(s)
end

function archivefile(file::String; inpath::String=DATA_PATH, inext::String = "",
  archivepath::String = ARCHIVE_PATH)

  (length(inext) > 0) && (file = "$file.$inext")
  inp::String = "$inpath\\$file"
  outp::String = "$archivepath\\$file.gz"

  open(GzipCompressorStream, outp, "w") do sout
    open(inp) do sin
      write(sout, sin)
    end
  end

  return nothing
end

function unarchivefile(file::String; outpath::String=DATA_PATH, outext::String = "",
  archivepath::String = ARCHIVE_PATH)

  (length(outext) > 0) && (file = "$file.$outext")
  outp::String = "$outpath\\$file"
  inp::String = "$archivepath\\$file.gz"

  open(GzipDecompressorStream, inp) do sin
    open(outp, "w") do sout
      write(sout, sin)
    end
  end

end

function recompressfiles(files = FILES_TO_RECOMPRESS; datapath=DATA_PATH)::Nothing

  for f ∈ files
    gz2lz4(f, location=datapath)
  end

  return nothing
end

recompressfile(file; datapath=DATA_PATH)::Nothing = recompressfiles(
  [file]; datapath=datapath)

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
### Working with cooler files, uses cool https://github.com/open2c/cooler

using PyCall
using DataFrames
using OffsetArrays, CoordinateTransformations, Rotations
using ImageTransformations

loadcooler() = pyimport("cooler")
loadcoolerfile(coolerfile; cooler=loadcooler()) = cooler.Cooler(coolerfile)

### load bin dirs
coolercoords(chrom, loc) = (replace(chrom, "chr" => ""), first(loc), last(loc))
function bindf(chrom, loc, cool; ext=0.5, saveload=true, projdir=getprojectdir())

    dffile =  joinpath(projdir, "results", "cooler", string("isletcooler_bindf_", chrom, "_", first(loc), "_", last(loc), "", ".tsv.gz"))
    if isnothing(cool) || !saveload
        if !isfile(dffile)
            error("Error: cannot load precalculated contact matrix bins $dffile, please load cooler file")
        end
        return CSV.read(dffile, DataFrame)
    end

    dx = Int(round((last(loc) - first(loc))*ext))
    xloc = (first(loc) - dx):(last(loc) + dx)

    bdf = cool.bins().fetch(coolercoords(chrom, xloc))
    cols = Symbol.(bdf.columns)
    data = [[v for v in bdf[c]] for c in cols]
    df = DataFrame(data, cols)
    df.binmid = div.(df.start .+ df.end, 2)
    df

    if saveload
        mkpath(dirname(dffile))
        CSV.write(dffile, df, compress=true)
    end
    df
end

## cooler matrix
function coolmatrix(chrom, loc, cool; balance="KR", ext=0.5, saveload=true, projdir=getprojectdir())

    matfile =  joinpath(projdir, "results", "cooler", string("isletcooler_", chrom, "_", first(loc), "_", last(loc), "", ".tsv.gz"))
    if isnothing(cool) || !saveload
        if !isfile(matfile)
            error("Error: cannot load precalculated contact matrix file $matfile, please load cooler file")
        end
        return Matrix{Float64}(CSV.read(matfile, DataFrame)[!, 2:end])
    end

    dx = Int(round((last(loc) - first(loc))*ext))
    xloc = (first(loc) - dx):(last(loc) + dx)
    CM = cool.matrix(balance=balance).fetch(coolercoords(chrom, xloc))

    if saveload
        mkpath(dirname(matfile))
        CSV.write(matfile, DataFrame(CM, :auto), compress=true)
    end
    CM
end

### rotate matrix to be 2d
function rotcooler(xp, CM ; kwargs...)
    R = warp(CM, recenter(RotMatrix(-π/4), [0, 0]) ; kwargs...)
    n = cld(size(R, 1), 2)
    R = R.parent[(n+1):end, :]
    nxp = range(first(xp), last(xp), length=size(R, 2))
    dx = xp[2] - xp[1]
    nyp = dx*(1:size(R, 1))
    nxp, nyp, R
end

## quantify islet loops
function quantisletloops(hicdata, cool, binsize=5000, ext=2, saveload=true, projdir=getprojectdir())
    datafile = joinpath(projdir, "results", "cooler", string("islet_hicquant.tsv.gz"))
    if isnothing(cool) || !saveload
        if !isfile(datafile)
            error("Error: cannot load precalculated quantification of islet loops $datafile, please load cooler file")
        end
        return CSV.read(datafile, DataFrame)
    end

    dx = binsize*ext
    μ = Float64[]
    @showprogress for row in eachrow(hicdata)
        locA = (row.startA - dx):(row.stopA + dx)
        locB = (row.startB - dx):(row.stopB + dx)
        m = cool.matrix(balance="KR").fetch(coolercoords(row.chrA, locA), coolercoords(row.chrB, locB))
        push!(μ, mean(m))
    end

    hicdata.LoopQuant = μ
    if saveload
        mkpath(dirname(datafile))
        CSV.write(datafile, hicdata, delim='\t', compress=true)
    end
    hicdata
end

### Plot HiC matrix over a region
function coolerpanel(chrom, loc, isletcooler, hicdata, regioncoords;  ext=0.75, plotroi=true)

    bincoords = bindf(chrom, loc, isletcooler, ext=ext)
    CM = coolmatrix(chrom, loc, isletcooler, ext=ext)
    nxp, nyp, R = rotcooler(bincoords.binmid, CM)
    cl=(0, 2.3139)
    ph = heatmap(nxp, nyp, log10.(R .+ 1), ylims=(0, 4.0e+5), colorbar=false, c=:inferno, xticks=false, clims=cl) ## clims hardcoded to compare to HK1 scale
    # cl = zlims()
    # @show cl
    yticks!((0:200:400)*1e+3, string.(0:.2:.4))

    hicloops = @subset(hicdata, :Study .== "IsletHiC", :chrA .== chrom, (:startA .∈ Ref(loc)) .| (:stopB .∈ Ref(loc)) .| (:startB .∈ Ref(loc)) .| (:stopA .∈ Ref(loc))) ## islet loops
    plottriloops!(hicloops, c=:white, lab="Islet Loops", plotlines=true)

    plot!(ylabel="Distance (Mb)")

    plotroi && vline!([mean(last(regioncoords["roi"]))], c=:white, lab="ROI", ls=:dash, title="Islets HiC")
    pg = plotgenemodels(loc, genecoords, fs=6)
    plot!(ylims=(-2, 2), yticks=false, yaxis=false, xticks=true, xlabel=string(chrom, ":", loc[1], "-", loc[2]), bottom_margin=5mm)
    plotroi && vline!([mean(last(hk1coords["roi"]))], c=:black, lab="", ls=:dash)
    pblank = plot([0], [0], zcolor=[NaN], colorbar=true, clims=cl, lab="", framestyle=:none, c=:inferno)
    lt = @layout [[a{0.65h}  ; b] c{0.05w}]

    p = plot(ph, pg, layout=lt, pblank, size=(900, 290), link=:x, xlims=(first(loc), last(loc)), fmt=:png, background_color_legend=:black, foreground_color_legend=:white, legendfontcolor=:white,
    left_margin=5mm, fontfamily="helvetica")

end

### plot triangluar loops over HiC Matrix
function plottriloops!(df; w=15000, np=100, c=:white, ss=1, plotlines=true, lab="")

    ps = Plots.partialcircle(0, 2π, np, w)
    ofs = ss*w/sqrt(2)

    for (k, row) in enumerate(eachrow(df))
        l = (row.startA + row.stopA)/2
        r = (row.startB + row.stopB)/2

        m, s = coordconv(l, r)
        plotlines && plot!([l, m - ofs], [0, s-ofs*sqrt(2)], c=c, lab="")
        plotlines && plot!([m+ofs, r], [s-ofs*sqrt(2), 0], c=c, lab="")
        plot!(first.(ps) .+ m, last.(ps).*sqrt(2) .+ s, c=c, lab=ifelse(k == 1, lab, ""))
    end
    plot!()
end

function coordconv(start, stop)
    m = (start+stop)/2
    s = (stop - start)/sqrt(2)
    m, s
end

iv(a, b) = a:b
function hk1_loop_summary(hicdata, hk1coords ; trf = x -> log10(x + 1))

    chrom, loc = hk1coords["roi"]
    ### hk1 loop

    hk1loops = @subset(hicdata, :chrA .== chrom, :chrB .== chrom, .!isempty.(intersect.(iv.(:startA, :stopB), Ref(loc))),  ((:stopB .- :startA) .< 1e+6))
    hk1loops.roidist = abs.(first(loc) .- hk1loops.startA) .+ abs.(last(loc) .- hk1loops.startB)
    sort!(hk1loops, :roidist) # take quantification of loop that tightest around the ROI
    lq = trf(first(hk1loops.LoopQuant))

    bins = range(0, 2.25, length=80)
    q = trf.(hicdata.LoopQuant[.!isnan.(hicdata.LoopQuant)])

    ph = stephist(q, bins = bins, c=:steelblue, fill=0, fillalpha=0.2, lab="", ylabel="Frequency", title="Histogram of log Contact Frequencies\nAll Loops Human Islets")
    vline!([lq], c=:black, lab="HK1 Loop", leg=:topleft, ls=:dash)
    ## ecdf
    ec = ecdf(q)

    pe = plot(bins, ec(bins), c=:steelblue, lab="", ylabel="Cumlative Proprtion", title="Empirical cdf log Contact Frequencies\nAll Loops Human Islets")
    yl = ylims()
    xl = xlims()

    plot!([lq, lq], [first(yl), ec(lq)], c=:black, lab="", ls=:dash)
    plot!([first(xl), lq], [ec(lq), ec(lq)], c=:black, lab=string("HK1 Loop: ", round(ec(lq), digits=2)), ls=:dash, leg=:left)
    plot!(xlims=xl, ylims=yl)
    plot(ph, pe, size=(790, 250), xlabel="log10 Loop Contact frequencies", bottom_margin=5mm, fontfamily="helvetica", left_margin=5mm, titlefont=font(12, "helvetica"), top_margin=5mm)

end


text:
using RadialBasisFunctionModels
using MLJBase

r = RadialBasisFunctionModels.RBFInterpolator(; kernel_name = :multiquadric, kernel_args = [1,])
X,y = @load_boston

R = machine(r, X, y)
fit!(R)

predict(R, X)


text:
# ITER parameters
δ = 0.33
ϵ = 0.32
κ = 1.7
B0 = 5.3
R0 = 6.2
qstar = 1.57
alpha = -0.155

S1 = solovev(B0, R0, ϵ, δ, κ, alpha, qstar,B0_dir=1,Ip_dir=1)
S2 = solovev(B0, R0, ϵ, δ, κ, alpha, qstar,B0_dir=-1,Ip_dir=1)
S3 = solovev(B0, R0, ϵ, δ, κ, alpha, qstar,B0_dir=1,Ip_dir=-1)
S4 = solovev(B0, R0, ϵ, δ, κ, alpha, qstar,B0_dir=-1,Ip_dir=-1)

btip1 = (1,1)
btip2 = (-1,1)
btip3 = (1,-1)
btip4 = (-1,-1)

cc1 = 3
cc2 = 3
cc3 = 3
cc4 = 3

r1 = magnetic_axis(S1) .+ (0.1,0.0)
r2 = magnetic_axis(S2) .+ (0.1,0.0)
r3 = magnetic_axis(S3) .+ (0.1,0.0)
r4 = magnetic_axis(S4) .+ (0.1,0.0)

test_data = ((cc1,S1,r1,btip1), (cc2,S2,r2,btip2), (cc3,S3,r3,btip3), (cc4,S4,r4,btip4))

@testset "Solov'ev Tests" begin

    @testset "Bt-Ip Signs" begin
        for (cc, S, r, btip) in test_data
            Bt, Ip = (sign(Bfield(S,r...)[2]), sign(Jfield(S,r...)[2]))
            @test (Bt, Ip) == btip
        end
    end

    @testset verbose = true "CurlB == μ₀J" begin
        for (cc, S, r, btip) in test_data
            @test curlB(S,r...) ≈ mu0*Jfield(S,r...) rtol=0.01
            @test Equilibrium.curlB_autodiff(S,r...) ≈ mu0*Jfield(S,r...) rtol=0.01
        end
    end

    @testset verbose = true "3D tests" begin
        for (cc, S, r, btip) in test_data
            @test Bfield(S,r...) ≈ Bfield(S,r[1],zero(r[1]),r[2]) rtol=0.01
            @test Jfield(S,r...) ≈ Jfield(S,r[1],zero(r[1]),r[2]) rtol=0.01

            @test gradB(S,r...) ≈ gradB(S,r[1],zero(r[1]),r[2]) rtol=0.01
            @test gradB(S,r[1],zero(r[1]),r[2]) ≈ Equilibrium.gradB_autodiff(S,r[1],zero(r[1]),r[2]) rtol=0.01

            @test curlB(S,r...) ≈ curlB(S,r[1],zero(r[1]),r[2]) rtol=0.01
            @test curlB(S,r[1],zero(r[1]),r[2]) ≈ Equilibrium.curlB_autodiff(S,r[1],zero(r[1]),r[2]) rtol=0.01
        end
    end

    @testset verbose = true "COCOS Consistancy" begin
        for (cc, S, r, btip) in test_data
            psi = S(r...)
            sigma_F = sign(poloidal_current(S,psi))
            sigma_pprime = sign(pressure_gradient(S,psi))
            sigma_q = sign(safety_factor(S, psi))
            sigma_dpsi = -sign(psi) # Solov'ev psi at boundary zero
            @test check_cocos(btip[1],btip[2], sigma_F, sigma_pprime, sigma_q, sigma_dpsi, cc; verbose=true)
        end
    end

    @testset "ForwardDiff" begin
        function opti(alpha)
            S = solovev(B0, R0, ϵ, δ, κ, alpha, qstar, B0_dir=1, Ip_dir=1)
            return (S.beta_p - 2)^2
        end

        @test (ForwardDiff.derivative(opti, alpha) !== nothing)
    end

end


text:
dfg = CloudGraphsDFG{NoSolverParams}("localhost", 7474, "neo4j", "test",
                            "testUser", "testRobot", "testSession",
                            nothing,
                            nothing,
                            IncrementalInference.decodePackedType,
                            IncrementalInference.rebuildFactorMetadata!,
                            solverParams=NoSolverParams())

# Nuke the user
clearUser!!(dfg)
@test lsSessions(dfg) == []
@test lsRobots(dfg) == []
@test !(Symbol(dfg.userId) in map(u -> u.id, lsUsers(dfg)))
# Create sentinel nodes using shortcut
createDfgSessionIfNotExist(dfg)
@test map(s -> s.id, lsSessions(dfg)) == [Symbol(dfg.sessionId)]
@test map(s -> s.id, lsRobots(dfg)) == [Symbol(dfg.robotId)]
@test Symbol(dfg.userId) in map(u -> u.id, lsUsers(dfg))
# Test that we can call it again.
createDfgSessionIfNotExist(dfg)
# And nuke it so we can try the longer functions.
clearUser!!(dfg)

# User, robot, and session
# TODO: Make easier ways to initialize these.
user = User(Symbol(dfg.userId), "Bob Zack", "Description", Dict{Symbol, String}())
robot = Robot(Symbol(dfg.robotId), user.id, "Test robot", "Description", Dict{Symbol, String}())
session = Session(Symbol(dfg.sessionId), robot.id, user.id, "Test Session", "Description", Dict{Symbol, String}())

@test createUser(dfg, user) == user
@test createRobot(dfg, robot) == robot
@test createSession(dfg, session) == session
@test map(s -> s.id, lsSessions(dfg)) == [Symbol(dfg.sessionId)]
@test map(s -> s.id, lsRobots(dfg)) == [Symbol(dfg.robotId)]
@test Symbol(dfg.userId) in map(u -> u.id, lsUsers(dfg))

# Test errors
dfgError = deepcopy(dfg)
# User/robot/session ID's can't start with numbers and can't have spaces.
dfgError.userId = "1testNope"
user = User(Symbol(dfgError.userId), "Bob Zack", "Description", Dict{Symbol, String}())
@test_throws Exception createUser(dfgError, user)

@test getUserData(dfg) == Dict{Symbol, String}()
@test getRobotData(dfg) == Dict{Symbol, String}()
@test getSessionData(dfg) == Dict{Symbol, String}()

# User/robot/session data
user.data = Dict{Symbol, String}(:a => "Hello", :b => "Goodbye")
robot.data = Dict{Symbol, String}(:c => "Hello", :d => "Goodbye")
session.data = Dict{Symbol, String}(:e => "Hello", :f => "Goodbye")

setUserData(dfg, user.data)
setRobotData(dfg, robot.data)
setSessionData(dfg, session.data)
@test getUserData(dfg) == user.data
@test getRobotData(dfg) == robot.data
@test getSessionData(dfg) == session.data

# Add some nodes.
v1 = addVariable!(dfg, :a, ContinuousScalar, labels = [:POSE])
v2 = addVariable!(dfg, :b, ContinuousScalar, labels = [:POSE])
v3 = addVariable!(dfg, :c, ContinuousScalar, labels = [:LANDMARK])
f1 = addFactor!(dfg, [:a; :b], LinearConditional(Normal(50.0,2.0)) )
f2 = addFactor!(dfg, [:b; :c], LinearConditional(Normal(50.0,2.0)) )

sessions = lsSessions(dfg)
@test map(s -> s.id, sessions) == [session.id]

# Pull and solve this graph
dfgLocal = GraphsDFG{SolverParams}(params=SolverParams())
DistributedFactorGraphs.getSubgraph(dfg, union(ls(dfg), lsf(dfg)), true, dfgLocal)

# Confirm that with sentinels we still have the same graph (doesn't pull in the sentinels)
@test symdiff(ls(dfgLocal), ls(dfg)) == []
@test symdiff(lsf(dfgLocal), lsf(dfg)) == []

# Solve it
tree, smtasks = solveTree!(dfgLocal)

# Make sure we can copy and solve normal orphaned sessions.
dfgOrphaned = deepcopy(dfg)
dfgOrphaned.sessionId = "doesntexist"
# Don't create a session
# Add some nodes.
v1 = addVariable!(dfgOrphaned, :a, ContinuousScalar, labels = [:POSE])
v2 = addVariable!(dfgOrphaned, :b, ContinuousScalar, labels = [:POSE])
v3 = addVariable!(dfgOrphaned, :c, ContinuousScalar, labels = [:LANDMARK])
f1 = addFactor!(dfgOrphaned, [:a; :b], LinearConditional(Normal(50.0,2.0)) )
f2 = addFactor!(dfgOrphaned, [:b; :c], LinearConditional(Normal(50.0,2.0)) )
# Solve it
dfgLocal = GraphsDFG{SolverParams}(params=SolverParams())
DistributedFactorGraphs.getSubgraph(dfgOrphaned, union(ls(dfgOrphaned), lsf(dfgOrphaned)), true, dfgLocal)
tree, smtasks = solveTree!(dfgLocal)

# If this passes without errors and we solve the graph, then all good.


text:
using Test
import Pluto
import Pluto: update_run!, WorkspaceManager, ClientSession, ServerSession, Notebook, Cell

@testset "Bonds" begin

    🍭 = ServerSession()
    🍭.options.evaluation.workspace_use_distributed = false
    fakeclient = ClientSession(:fake, nothing)
    🍭.connected_clients[fakeclient.id] = fakeclient

    @testset "AbstractPlutoDingetjes.jl" begin
        🍭.options.evaluation.workspace_use_distributed = true
        notebook = Notebook([
                Cell("""
                begin
                    import Pkg
                    try
                        Pkg.UPDATED_REGISTRY_THIS_SESSION[] = true
                    catch; end
                    Pkg.activate(mktempdir())
                    Pkg.add(Pkg.PackageSpec(name="AbstractPlutoDingetjes", rev="30f4f76"))
                    import AbstractPlutoDingetjes
                    const APD = AbstractPlutoDingetjes
                    import AbstractPlutoDingetjes.Bonds
                end
                """),
                Cell("""
                begin
                    struct HTMLShower
                        f::Function
                    end
                    Base.show(io::IO, m::MIME"text/html", hs::HTMLShower) = hs.f(io)
                end
                """),
                Cell("""
                APD.is_inside_pluto()
                """),
                Cell("""
                HTMLShower() do io
                    write(io, string(APD.is_inside_pluto(io)))
                end
                """),
                Cell("""
                HTMLShower() do io
                    write(io, string(APD.is_supported_by_display(io, Bonds.initial_value)))
                end
                """),
                Cell("""
                HTMLShower() do io
                    write(io, string(APD.is_supported_by_display(io, Bonds.transform_value)))
                end
                """),
                Cell("""
                HTMLShower() do io
                    write(io, string(APD.is_supported_by_display(io, Bonds.validate_value)))
                end
                """),
                Cell("""
                HTMLShower() do io
                    write(io, string(APD.is_supported_by_display(io, sqrt)))
                end
                """),
                Cell("""
                @bind x_simple html"<input type=range>"
                """),
                Cell("""
                x_simple
                """),
                Cell("""
                begin
                    struct OldSlider end

                    Base.show(io::IO, m::MIME"text/html", os::OldSlider) = write(io, "<input type=range value=1>")

                    Base.get(os::OldSlider) = 1
                end
                """),
                Cell("""
                @bind x_old OldSlider()
                """),
                Cell("""
                x_old
                """),
                Cell("""
                begin
                    struct NewSlider end

                    Base.show(io::IO, m::MIME"text/html", os::NewSlider) = write(io, "<input type=range value=1>")

                    Bonds.initial_value(os::NewSlider) = 1
                end
                """),
                Cell("""
                @bind x_new NewSlider()
                """),
                Cell("""
                x_new
                """),
                Cell("""
                begin
                    struct TransformSlider end

                    Base.show(io::IO, m::MIME"text/html", os::TransformSlider) = write(io, "<input type=range value=1>")

                    Bonds.initial_value(os::TransformSlider) = "x"
                    Bonds.transform_value(os::TransformSlider, from_js) = repeat("x", from_js)
                end
                """),
                Cell("""
                @bind x_transform TransformSlider()
                """),
                Cell("""
                x_transform
                """),
                Cell("""
                begin
                    struct BadTransformSlider end

                    Base.show(io::IO, m::MIME"text/html", os::BadTransformSlider) = write(io, "<input type=range value=1>")

                    Bonds.initial_value(os::BadTransformSlider) = "x"
                    Bonds.transform_value(os::BadTransformSlider, from_js) = repeat("x", from_js)
                end
                """),
                Cell("""
                @bind x_badtransform BadTransformSlider()
                """),
                Cell("""
                x_badtransform
                """),
                Cell("""
                count = Ref(0)
                """),
                Cell("""
                @bind x_counter NewSlider() #or OldSlider(), same idea
                """),
                Cell("""
                let
                    x_counter
                    count[] += 1
                end
                """),
                Cell("""
                @assert x_old == 1
                """),
                Cell("""
                @assert x_new == 1
                """),
                Cell("""
                @assert x_transform == "x"
                """),
            ])
        fakeclient.connected_notebook = notebook

        update_run!(🍭, notebook, notebook.cells)

        @test noerror(notebook.cells[1])
        @test noerror(notebook.cells[2])
        @test noerror(notebook.cells[3])
        @test noerror(notebook.cells[4])
        @test noerror(notebook.cells[5])
        @test noerror(notebook.cells[6])
        @test noerror(notebook.cells[7])
        @test noerror(notebook.cells[8])

        @test notebook.cells[3].output.body == "true"
        @test notebook.cells[4].output.body == "true"
        @test notebook.cells[5].output.body == "true"
        @test notebook.cells[6].output.body == "true"
        @test notebook.cells[7].output.body == "false"
        @test notebook.cells[8].output.body == "false"

        @test noerror(notebook.cells[9])
        @test noerror(notebook.cells[10])
        @test noerror(notebook.cells[11])
        @test noerror(notebook.cells[12])
        @test noerror(notebook.cells[13])
        @test noerror(notebook.cells[14])
        @test noerror(notebook.cells[15])
        @test noerror(notebook.cells[16])
        @test noerror(notebook.cells[17])
        @test noerror(notebook.cells[18])
        @test noerror(notebook.cells[19])
        @test noerror(notebook.cells[20])
        @test noerror(notebook.cells[21])
        @test noerror(notebook.cells[22])
        @test noerror(notebook.cells[23])
        @test noerror(notebook.cells[24])
        @test noerror(notebook.cells[25])
        @test noerror(notebook.cells[26])
        @test noerror(notebook.cells[27])
        @test noerror(notebook.cells[28])

        function set_bond_value(name, value, is_first_value=false)
            notebook.bonds[name] = Dict("value" => value, "is_first_value" => is_first_value)
            Pluto.set_bond_values_reactive(;
                session=🍭,
                notebook=notebook,
                bound_sym_names=[name],
                run_async=false,
            )
        end

        @test notebook.cells[10].output.body == "missing"
        set_bond_value(:x_simple, 1, true)
        @test notebook.cells[10].output.body == "1"

        @test notebook.cells[13].output.body == "1"
        set_bond_value(:x_old, 1, true)
        @test notebook.cells[13].output.body == "1"
        set_bond_value(:x_old, 99, false)
        @test notebook.cells[13].output.body == "99"

        @test notebook.cells[16].output.body == "1"
        set_bond_value(:x_new, 1, true)
        @test notebook.cells[16].output.body == "1"
        set_bond_value(:x_new, 99, false)
        @test notebook.cells[16].output.body == "99"

        @test notebook.cells[19].output.body == "\"x\""
        set_bond_value(:x_transform, 1, true)
        @test notebook.cells[19].output.body == "\"x\""
        set_bond_value(:x_transform, 3, false)
        @test notebook.cells[19].output.body == "\"xxx\""

        @test notebook.cells[22].output.body != "missing"
        set_bond_value(:x_badtransform, 1, true)
        @test notebook.cells[22].output.body != "missing"
        set_bond_value(:x_badtransform, 7, false)
        @test notebook.cells[22].output.body != "missing"

        @test notebook.cells[25].output.body == "1"
        set_bond_value(:x_counter, 1, true)
        @test notebook.cells[25].output.body == "1"
        set_bond_value(:x_counter, 7, false)
        @test notebook.cells[25].output.body == "2"

        WorkspaceManager.unmake_workspace((🍭, notebook))
        🍭.options.evaluation.workspace_use_distributed = false

        @test jl_is_runnable(notebook.path)

        # should_not_run = """
        # @assert 1 == 2
        # """
        # p = tempname()
        # write(p, should_not_run)
        # @test !jl_is_runnable(p)

    end
end


text:
@testset "celegansneural" begin

s = open(s -> read(s, String), "gml/celegansneural.gml")
parse_dict(s; debug=false)
@time x = parse_dict(s; debug=false)
@test length(x) == 2
@test x[:Creator] == "Mark Newman on Thu Aug 31 12:59:09 2006"
@test length(x[:graph][1]) == 3
@test length(x[:graph][1][:node]) == 297
@test length(x[:graph][1][:edge]) == 2359

end


text:
#=
 HOW TO SHOW AN S-EXPRESSION
 ≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡

  • cap lines to 80 characters (if possible)
  • overflow all s-expressions if not
  • if head is short [≤ 3], then keep first argument on the same line

 performance isn’t a concern but if it’s too shitty we can consider ropes
=#

struct ShowListContext
    indent::Int
    limit::Int
end

# arbitrarily decide that there’s always room for 5 more
space(ctx::ShowListContext) = max(5, ctx.limit - ctx.indent)

# performance is really not our concern
sindent(ctx::ShowListContext) = " " ^ ctx.indent

# “if we append or prepend, we should indent, not indented”
# no. indent is also a noun. indented makes it clear what this does
indented(ctx::ShowListContext, i=2) = ShowListContext(ctx.indent + i, ctx.limit)

# putting an ‘s’ in front of all functions that return strings is terrible style
# but sometimes terrible style is the best we’ve got
sprintwidth(α) = sum(charwidth(c) for c in unparse(α))

spprintall(ctx::ShowListContext, α) = join((β -> spprint(ctx, β)) ∘ α, '\n')

spprint(ctx::ShowListContext, α) = sindent(ctx) * unparse(α)

function spprint(ctx::ShowListContext, α::Cons)
    if sprintwidth(α) ≤ space(ctx) || length(α) ≤ 1
        # verbatim
        sindent(ctx) * unparse(α)
    elseif (spw = sprintwidth(car(α))) ≤ 3 && length(α) ≥ 3
        # (car cadr ¶ ...) format
        ctxi = indented(ctx, spw + 2)
        sindent(ctx) * "(" * unparse(car(α)) * " " * unparse(cadr(α)) * "\n" *
                spprintall(ctxi, cddr(α)) * ")"
    else
        # (car ¶ ...) format
        ctxi = indented(ctx, 2)
        sindent(ctx) * "(" * unparse(car(α)) * "\n" *
                spprintall(ctxi, cdr(α)) * ")"
    end
end

spprint(α::List) = spprint(ShowListContext(0, 80), α)

function Base.show(io::IO, α::List)
    if get(io, :multiline, false)
        print(io, spprint(α))
    else
        print(io, unparse(α))
    end
end


text:
export solve_modelKTR
function solve_modelKTR(nlp_stop, X :: PDataK, α:: T) where T
    # target value should be close to satisfy α=||d||
    target = [ ( abs( α - X.norm_dirs[i]) )   for i = 1 : X.nshifts ];
    # pick the closest shift to the target within positive definite H+λI
    # before, check the shift = 0 for direction within the trust region

    if (X.indmin == 0)  #  try Newton's direction
        X.indmin = 1
        if (X.positives[1]==1) & (X.norm_dirs[1] <= α)
            X.d = X.xShift[:,1]
            X.λ = 0.0
            return X.d, NaN * rand(length(X.d)), X.λ
        end
    end
    start = X.indmin
    bidon,indmin = findmin(target[X.positives[start:end]])
    X.indmin = start + indmin - 1
    p_imin = X.positives[X.indmin]
    X.d = X.xShift[:,p_imin]
    # Refine this rough constrained direction

    X.λ = X.shifts[p_imin]
    #println("In SolveModel   p_imin = $(p_imin),  α = $α")
    #println("In SolveModel   $(X.Ndirs[p_imin]-α) , $(X.Ndirs[p_imin+1]-α)  ,  α = $α")

    #end
    return X.d, NaN * rand(length(X.d)), X.λ
end

# To replace.
#    1. Form the signed ( α - X.norm_dirs[i]) ) within positives of the form +++...----
#    2. If starts with a +, interpolate the interval where sign changes
#       else (hard case) obtain from Krylov shift solver negative curvature directions
#            and extrapolate to the boundary the short (longest) direction
#            in the descent negative curvature direction
#
# Keep on treating λ=0 special for Newton's direction (no extrapolation to boundary)


text:
# AIML Core

# Thoughts: data structure for currently loaded aiml graph

max_stars_repo_name: seantfish/chatbot

max_issues_repo_name: seantfish/chatbot

max_forks_repo_name: seantfish/chatbot


text:
function ML_get_f_g(featuresTP::Matrix{Float64}, labels::Vector{Float64}, W::Matrix{Float64}, index::Int64)
	numFeatures = size(featuresTP)[1]
	numClasses = div(length(W),numFeatures)
	W=reshape(W,(numFeatures,numClasses))
	class=labels[index]
	x = featuresTP[:,index]
	a=exp(x'*W)
	aDb=a/sum(a)
	g=x*aDb
	g[:,round(Int,class)]-=x
	(-log(aDb[round(Int,class)]),reshape(g, (numFeatures*numClasses,1)),aDb)
end

function ML_restore_gradient(featuresTP::Matrix{Float64}, labels::Vector{Float64}, aDb, index::Int64)
	numFeatures = size(featuresTP)[1]
	numClasses = 129
	class=labels[index]
	x = featuresTP[:,index]
	g=x*aDb
	g[:,class]-=x
	reshape(g, (numFeatures*numClasses,1))
end

function ML_get_f(features::Matrix{Float64}, labels::Vector{Float64}, W::Matrix{Float64})
	numFeatures = size(features)[2]
	numClasses = div(length(W),numFeatures)
	W=reshape(W,(numFeatures,numClasses))
	f=0
	@inbounds @simd for k=1:length(labels)
		class=labels[k]
		x=features[k,:]
		a=exp(x*W)
		aDb=a/sum(a)
		f+=log(aDb[round(Int,class)])
	end
	(-f/length(labels))
end

# Returns f, percent correctly classified
function ML_for_output(features::Matrix{Float64}, labels::Vector{Float64}, W::Matrix{Float64})
	numFeatures = size(features)[2]
	numClasses = div(length(W),numFeatures)
	W=reshape(W,(numFeatures,numClasses))
	g=zeros(size(W))
	f=0
	pcc=0
	@inbounds @simd for k=1:length(labels)
		class=labels[k]
		x=features[k,:]
		a=exp(x*W)
		(v,chosenclass) = findmax(a);
		if chosenclass == class
			pcc+=1
		end
		aDb=a/sum(a)
		f+=log(aDb[round(Int,class)])
	end
	(-f/length(labels),pcc/length(labels))
end


text:
using  JITrench

# This macro automatically define the derivatives function.
macro diff!(ex)
    f_name = ex.args[1]
    df_name = Symbol(String(f_name) * "′")
    quote
        function $(esc(df_name))(x)
            x_var = Variable(x)
            y = $(esc(ex.args[1]))(x_var)
            backward!(y)
            return x_var.grad.values
        end
    end
end

f(x) = 2x^2 + 4x + 3

@diff! f(x)
# f′(x) = 4x + 4

f′(2)
# 12

f′(3)
# 16


text:
function pow(b, n, m)
    ret = 1
    while n != 0
        if n & 1 != 0
            ret = (ret*b) % m
        end
        b = b^2 % m
        n = n >> 1
    end
    ret
end

function solve()
    m = 1000000007
    n = parse(Int, readline())
    ret = pow(10, n, m) - 2*pow(9, n, m) + pow(8, n, m)
    ret %= m
    ret = (ret+m) % m
end

println(solve())

max_stars_repo_name: aishikawa/atcoder-julia

max_issues_repo_name: aishikawa/atcoder-julia

max_forks_repo_name: aishikawa/atcoder-julia