text:
module MonotoneCubicInterpolation
    export interpolateCubicHermite, interpolate_derivative_CubicHermite
    """
    Interpolate using cubic Hermite splines. The breakpoints in arrays xbp and ybp are assumed to be sorted.
    Evaluate the function in all points of the array xeval.
    Methods:
        :linear                yuck
        :finitedifference      classic cubic interpolation, no tension parameter
        :cardinal              cubic cardinal splines, uses tension parameter which must be between [0,1]
        :fritschcarlson"        monotonic - tangents are first initialized, then adjusted if they are not monotonic
        :fritschbutland        monotonic - faster algorithm (only requires one pass) but somewhat higher apparent "tension"
    Sources:
        Fritsch & Carlson (1980), "Monotone Piecewise Cubic Interpolation", doi:10.1137/0717021.
        Fritsch & Butland (1984), "A Method for Constructing Local Monotone Piecewise Cubic Interpolants", doi:10.1137/0905021.
    """

    function interpolateCubicHermite(xeval, xbp, ybp, method, tension)
        # first we need to determine tangents (m)
        n = length(xbp)
        m, delta = calcTangents(xbp, ybp, method, tension)

        c = zeros(n-1)
        d = zeros(n-1)
        for k = 1:n-1
            if method == :linear
                m[k] = delta[k]
                c[k] = 0
                d[k] = 0
                continue
            end
            xdiff = xbp[k+1] - xbp[k]
            c[k] = (3*delta[k] - 2*m[k] - m[k+1]) / xdiff
            d[k] = (m[k] + m[k+1] - 2*delta[k]) / xdiff / xdiff
        end
        len = length(xeval)
        f = zeros(len)
        k = 1
        for i = 1:len
            x = xeval[i]
            if (x < xbp[1]) || (x > xbp[n])
                throw(string("interpolateCubicHermite: x value ", x ," outside breakpoint range [", xbp[1] ,", " ,xbp[n], "]"))
            end
            while (k < n) && (x > xbp[k+1])
                k+=1
            end
            xdiff = x - xbp[k]
            f[i] = ybp[k] + m[k]*xdiff + c[k]*xdiff*xdiff + d[k]*xdiff*xdiff*xdiff
        end
        return f
    end

    function calcTangents(x, y, method, tension)
        n = length(x)
        delta = zeros(n-1)
        m = zeros(n)
        for k = 1:n-1 #   (var k=0; k < n-1; k++) {
            deltak = (y[k+1] - y[k]) / (x[k+1] - x[k])
            delta[k] = deltak
            if k == 1    # left endpoint, same for all methods
                m[k] = deltak
            else
                if method == :cardinal
                    m[k] = (1 - tension) * (y[k+1] - y[k-1]) / (x[k+1] - x[k-1])
                elseif method == :fritschbutland
                    alpha = (1 + (x[k+1] - x[k]) / (x[k+1] - x[k-1])) / 3  # Not the same alpha as below.
                    if delta[k-1] * deltak <= 0
                        m[k] = 0
                    else
                        m[k] = delta[k-1] * deltak / (alpha*deltak + (1-alpha)*delta[k-1])
                    end
                elseif method == :fritschcarlson
                    # If any consecutive secant lines change sign (i.e. curve changes direction), initialize the tangent to zero.
                    # This is needed to make the interpolation monotonic. Otherwise set tangent to the average of the secants.
                    if delta[k-1] * deltak < 0
                        m[k] = 0
                    else
                        m[k] = (delta[k-1] + deltak) / 2
                    end
                elseif method == :finitedifference    # FiniteDifference
                    m[k] = (delta[k-1] + deltak) / 2
                elseif method == :linear
                else
                    throw("Unknown method")
                end
            end
        end
        m[n] = delta[n-1]

        if method == :fritschcarlson
            """
            Fritsch & Carlson derived necessary and sufficient conditions for monotonicity in their 1980 paper. Splines will be
            monotonic if all tangents are in a certain region of the alpha-beta plane, with alpha and beta as defined below.
            A robust choice is to put alpha & beta within a circle around origo with radius 3. The FritschCarlson algorithm
            makes simple initial estimates of tangents and then does another pass over data points to move any outlier tangents
            into the monotonic region. FritschButland & Steffen algorithms make more elaborate first estimates of tangents that
            are guaranteed to lie in the monotonic region, so no second pass is necessary.
            """
            # Second pass of FritschCarlson: adjust any non-monotonic tangents.
            for k = 1:n-1 #(var k=0; k < n-1; k++) {
                deltak = delta[k]
                if (deltak == 0)
                    m[k] = 0
                    m[k+1] = 0
                    continue
                end
                alpha = m[k] / deltak
                beta = m[k+1] / deltak
                tau = 3 / sqrt(alpha.^2 + beta.^2);
                if (tau < 1)       # if we're outside the circle with radius 3 then move onto the circle
                    m[k] = tau * alpha * deltak
                    m[k+1] = tau * beta * deltak
                end
            end
        end
        return m, delta
    end

end # module


text:
using SocialSim

sweep_temp(20,5,50)

ext: jl

max_stars_repo_name: luca-mm/Bachelor-Thesis

max_issues_repo_name: luca-mm/Bachelor-Thesis

max_forks_repo_name: luca-mm/Bachelor-Thesis


text:
module ESDLZarr
using ESDL
import ZarrNative: ZGroup, zopen
import ESDL.CubeAPI: getNv, gettoffsnt, getvartype, Cube
import ESDL.Cubes: cubechunks, iscompressed
using Dates

struct ZarrConfig
  start_time::Date
  end_time::Date
  grid_x0::Int
  grid_y0::Int
  grid_height::Int
  grid_width::Int
  spatial_res::Float64
  temporal_res::Int
  chunk_sizes::NTuple{3,Int}
  compressed::Bool
end
function ZarrConfig(g::ZGroup,firstvar)
  timeax = g["time"]
  timevals = timeax[:]
  refdate = Date(split(timeax.attrs["units"]," ")[3])
  start_time = refdate+Dates.Day(first(timevals))-Day(4)
  end_time = refdate+Dates.Day(last(timevals))-Day(4)
  lonax = g["lon"]
  latax = g["lat"]
  sres = abs(latax[2][1]-latax[1][1])
  tres = timevals[2][1]-timevals[1][1]
  cs = g[firstvar].metadata.chunks
  compressed = !isa(g[firstvar].metadata.compressor,ZarrNative.NoCompressor)
  ZarrConfig(start_time,end_time,0,0,length(latax),2*length(latax),sres,tres,cs,compressed)
end

struct ZarrCube <: ESDL.CubeAPI.UCube
  path::String
  group::ZGroup
  dataset_files::Vector{String}
  var_name_to_var_index::Dict{String,Int}
  config::ZarrConfig
end
function ZarrCube(p::String)
  Cube(zopen(p))
end
function Cube(zg::ZGroup)
  spattempars = filter(v->ndims(v[2])==3,zg.arrays)
  varlist = collect(keys(spattempars))
  vni = Dict(i[2]=>i[1] for i in enumerate(varlist))
  ZarrCube(zg.storage.folder,zg,varlist,vni,ZarrConfig(zg,varlist[1]))
end
getvartype(c::ZarrCube, n::String)=eltype(c.group[n])
cubechunks(c::ZarrCube)=c.config.chunk_sizes
iscompressed(c::ZarrCube)=c.config.compressed
function ESDL.Cubes._read(s::ESDL.CubeAPI.SubCube{<:Any,ZarrCube},t::AbstractArray,r::CartesianIndices)
  grid_y1,grid_y2,grid_x1,grid_x2 = s.sub_grid
  y1,i1,y2,i2,ntime,NpY           = s.sub_times
  toffs,nt= gettoffsnt(s,r)
  toffs = (y1 - year(s.cube.config.start_time))*NpY + i1 + toffs
  rcor = CartesianIndices((r.indices[1].+(grid_x1-1), r.indices[2].+(grid_y1-1), (toffs:toffs.+nt-1)))
  #voffs,nv = getNv(r)
  singvar_zarr(t,s.cube,s.variable,rcor)
end
function ESDL.Cubes._read(s::ESDL.CubeAPI.SubCubeV{<:Any,ZarrCube},t::AbstractArray,r::CartesianIndices)
  grid_y1,grid_y2,grid_x1,grid_x2 = s.sub_grid
  y1,i1,y2,i2,ntime,NpY           = s.sub_times

  toffs,nt= gettoffsnt(s,r)
  toffs = (y1 - year(s.cube.config.start_time))*NpY + i1 + toffs
  voffs,nv = getNv(r)
  #@show (r.indices[1].+(grid_x1-1), r.indices[2].+(grid_y1-1), toffs:toffs+nt-1)
  rcor = CartesianIndices((r.indices[1].+(grid_x1-1), r.indices[2].+(grid_y1-1), toffs:toffs+nt-1))
  #@show rcor.indices
  for (iiv,iv)=enumerate((voffs+1):(voffs+nv))
    singvar_zarr(view(t,:,:,:,iiv),s.cube,s.variable[iv],rcor)
  end
end
function ESDL.Cubes._read(s::ESDL.CubeAPI.SubCubeStatic{<:Any,ZarrCube},t::AbstractArray,r::CartesianIndices)
  grid_y1,grid_y2,grid_x1,grid_x2 = s.sub_grid
  y1,i1,y2,i2,ntime,NpY           = s.sub_times
  toffs = (y1 - year(s.cube.config.start_time))*NpY + i1
  rcor = CartesianIndices((r.indices[1].+(grid_x1-1), r.indices[2].+(grid_y1-1),(toffs):(toffs)))
  singvar_zarr(t,s.cube,s.variable,rcor)
end

import ZarrNative
function singvar_zarr(outar,cube,variable,r)
  ZarrNative.readblock!(outar, cube.group[variable],r)
  mv  = cube.group[variable].metadata.fill_value
  if !isa(mv,Nothing)
    replace!(outar,mv=>missing)
  end
end

end # module


text:
"""
    transitiveclosure!(g, selflooped=false)

Compute the transitive closure of a directed graph, using DFS.
If `selflooped` is true, add self loops to the graph.

### Performance
Time complexity is ``\\mathcal{O}(|E||V|)``.

### Implementation Notes
This version of the function modifies the original graph.
"""
function transitiveclosure! end
@traitfn function transitiveclosure!(g::::IsDirected, selflooped=false)
    scc = strongly_connected_components(g)
    cg = condensation(g, scc)
    tp = reverse(topological_sort_by_dfs(cg))
    sr = [Vector{eltype(cg)}() for _ in vertices(cg)]

    x = selflooped ? 0 : 1
    for comp in scc
        for j in 1:(length(comp)-x)
            for k in (j+x):length(comp)
                add_edge!(g, comp[j], comp[k])
                add_edge!(g, comp[k], comp[j])
            end
        end
    end
    for u in tp
        for v in outneighbors(cg, u)
            union!(sr[u], sr[v], [v])
        end
    end
    for i in vertices(cg)
        for u in scc[i]
            for j in sr[i]
                for v in scc[j]
                    add_edge!(g, u, v)
                end
            end
        end
    end
    return g
end

"""
    transitiveclosure(g, selflooped=false)

Compute the transitive closure of a directed graph, using DFS.
Return a graph representing the transitive closure. If `selflooped`
is `true`, add self loops to the graph.

### Performance
Time complexity is ``\\mathcal{O}(|E||V|)``.

# Examples
```jldoctest
julia> using LightGraphs

julia> barbell = blockdiag(complete_digraph(3), complete_digraph(3));

julia> add_edge!(barbell, 1, 4);

julia> collect(edges(barbell))
13-element Array{LightGraphs.SimpleGraphs.SimpleEdge{Int64},1}:
 Edge 1 => 2
 Edge 1 => 3
 Edge 1 => 4
 Edge 2 => 1
 Edge 2 => 3
 Edge 3 => 1
 Edge 3 => 2
 Edge 4 => 5
 Edge 4 => 6
 Edge 5 => 4
 Edge 5 => 6
 Edge 6 => 4
 Edge 6 => 5

julia> collect(edges(transitiveclosure(barbell)))
21-element Array{LightGraphs.SimpleGraphs.SimpleEdge{Int64},1}:
 Edge 1 => 2
 Edge 1 => 3
 Edge 1 => 4
 Edge 1 => 5
 Edge 1 => 6
 Edge 2 => 1
 Edge 2 => 3
 Edge 2 => 4
 Edge 2 => 5
 Edge 2 => 6
 Edge 3 => 1
 Edge 3 => 2
 Edge 3 => 4
 Edge 3 => 5
 Edge 3 => 6
 Edge 4 => 5
 Edge 4 => 6
 Edge 5 => 4
 Edge 5 => 6
 Edge 6 => 4
 Edge 6 => 5
```
"""
function transitiveclosure(g::DiGraph, selflooped = false)
    copyg = copy(g)
    return transitiveclosure!(copyg, selflooped)
end

"""
    transitivereduction(g; selflooped=false)

Compute the transitive reduction of  a directed graph. If the graph contains
cycles, each strongly connected component is replaced by a directed cycle and
the transitive reduction is calculated on the condensation graph connecting the
components. If `selflooped` is true, self loops on strongly connected components
of size one will be preserved.

### Performance
Time complexity is ``\\mathcal{O}(|V||E|)``.

# Examples
```jldoctest
julia> using LightGraphs

julia> barbell = blockdiag(complete_digraph(3), complete_digraph(3));

julia> add_edge!(barbell, 1, 4);

julia> collect(edges(barbell))
13-element Array{LightGraphs.SimpleGraphs.SimpleEdge{Int64},1}:
 Edge 1 => 2
 Edge 1 => 3
 Edge 1 => 4
 Edge 2 => 1
 Edge 2 => 3
 Edge 3 => 1
 Edge 3 => 2
 Edge 4 => 5
 Edge 4 => 6
 Edge 5 => 4
 Edge 5 => 6
 Edge 6 => 4
 Edge 6 => 5

julia> collect(edges(transitivereduction(barbell)))
7-element Array{LightGraphs.SimpleGraphs.SimpleEdge{Int64},1}:
 Edge 1 => 2
 Edge 1 => 4
 Edge 2 => 3
 Edge 3 => 1
 Edge 4 => 5
 Edge 5 => 6
 Edge 6 => 4
```
"""
function transitivereducion end
@traitfn function transitivereduction(g::::IsDirected; selflooped::Bool=false)
    scc = strongly_connected_components(g)
    cg = condensation(g, scc)

    reachable = Vector{Bool}(undef, nv(cg))
    visited = Vector{Bool}(undef, nv(cg))
    stack = Vector{eltype(cg)}(undef, nv(cg))
    resultg = SimpleDiGraph{eltype(g)}(nv(g))

# Calculate the transitive reduction of the acyclic condensation graph.
    @inbounds(
    for u in vertices(cg)
        fill!(reachable, false) # vertices reachable from u on a path of length >= 2
        fill!(visited, false)
        stacksize = 0
        for v in outneighbors(cg,u)
      @simd for w in outneighbors(cg, v)
                if !visited[w]
                    visited[w] = true
                    stacksize += 1
                    stack[stacksize] = w
                end
            end
        end
        while stacksize > 0
            v = stack[stacksize]
            stacksize -= 1
            reachable[v] = true
      @simd for w in outneighbors(cg, v)
                if !visited[w]
                    visited[w] = true
                    stacksize += 1
                    stack[stacksize] = w
                end
            end
        end
# Add the edges from the condensation graph to the resulting graph.
  @simd for v in outneighbors(cg,u)
            if !reachable[v]
                add_edge!(resultg, scc[u][1], scc[v][1])
            end
        end
    end)

# Replace each strongly connected component with a directed cycle.
    @inbounds(
    for component in scc
        nvc = length(component)
        if nvc == 1
            if selflooped && has_edge(g, component[1], component[1])
                add_edge!(resultg, component[1], component[1])
            end
            continue
        end
        for i in 1:(nvc-1)
            add_edge!(resultg, component[i], component[i+1])
        end
        add_edge!(resultg, component[nvc], component[1])
    end)

    return resultg
end


text:
using IntervalArithmetic
using IntervalRootFinding

χ_λ(χ::AbstractArray, λ::Union{Float64,Interval{Float64}}) = map(χi -> 1.0 / ((1.0 / χi) + λ), χ)
χ_λ2(χ, λ) = map(χi -> 1.0 / ((1.0 / χi) + λ), χ)
function χ_λ!(χ_λ, χ, λ::Union{Float64,Interval{Float64}})
    for i in eachindex(χ_λ)
        χ_λ[i] = 1.0 / ((1.0 / χ[i]) + λ)
    end
end

dχ_λ(χ, λ::Union{Float64,Interval{Float64}}) = map(χi -> - ((1.0 / χi) + λ)^(-2), χ)

dΣch_λ_amp(G_plus_νq, γch, dχch_λ, qNorm) = -sum(G_plus_νq .* γch .* dχch_λ)*qNorm
dΣsp_λ_amp(G_plus_νq, γsp, dχsp_λ, qNorm) = -1.5*sum(G_plus_νq .* γsp .* dχsp_λ)*qNorm

function new_χλ(χ_in::SharedArray{Complex{Float64},2}, λ::Float64, sP::SimulationParameters)
    res = SharedArray{eltype(χ_in), ndims(χ_in)}(size(χ_in)...)
    if sP.χFillType == zero_χ_fill
        res[usable_ω,:] =  χ_λ(χ_in[uable_ω,:], λ)
    elseif sP.χFillType == lambda_χ_fill
        res =  χ_λ(χ_in, λ)
    else
        res[:] = deepcopy(χsp)
        res[usable_ω,:] =  χ_λ(χ_in[usable_ω,:], λ)
    end
    return res
end

function calc_λsp_rhs_usable(impQ_sp::ImpurityQuantities, impQ_ch::ImpurityQuantities, nlQ_sp::NonLocalQuantities, nlQ_ch::NonLocalQuantities, kGrid::T1, mP::ModelParameters, sP::SimulationParameters) where T1 <: ReducedKGrid
    usable_ω = intersect(nlQ_sp.usable_ω, nlQ_ch.usable_ω)
    χch_ω = kintegrate(kGrid, nlQ_ch.χ[usable_ω,:], dim=2)[:,1]
    @warn "currently using min(usable_sp, usable_ch) = min($(nlQ_sp.usable_ω),$(nlQ_ch.usable_ω)) = $(usable_ω) for all calculations. relax this?"

    #sh = get_sum_helper(usable_ω, sP, :b)
    iωn = 1im .* 2 .* (-sP.n_iω:sP.n_iω)[usable_ω] .* π ./ mP.β
    χch_ω_sub = subtract_tail(χch_ω, impQ_ch.tailCoeffs[3], iωn)
    χch_sum = real(sum_freq(χch_ω_sub, [1], Naive(), mP.β, corr=-impQ_ch.tailCoeffs[3]*mP.β^2/12)[1])

    #χupup_ω = 0.5*(χLocsp_ω .+ χLocch_ω)
    #χupup_DMFT_ω_sub = subtract_tail(χupup_ω[loc_range], E_kin_ED, iωn)
    #imp_density = real(sum_freq(χupup_DMFT_ω_sub, [1], Naive(), mP.β, corr=-E_kin_ED*mP.β^2/12))

    rhs = ((sP.tc_type_b != :nothing && sP.λ_rhs == :native) || sP.λ_rhs == :fixed) ? mP.n * (1 - mP.n/2) - χch_sum : real(impQ_ch.χ_loc + impQ_sp.χ_loc - χch_sum)
    #@info "tc:  $(sP.tc_type), rhs =  $rhs , $χch_sum , $(real(impQ_ch.χ_loc)) , $(real(impQ_sp.χ_loc)), $(real(χch_sum))"

    @info """Found usable intervals for non-local susceptibility of length
          sp: $(nlQ_sp.usable_ω), length: $(length(nlQ_sp.usable_ω))
          ch: $(nlQ_ch.usable_ω), length: $(length(nlQ_ch.usable_ω))
          usable: $(usable_ω), length: $(length(usable_ω))
          χch sum = $(χch_sum), rhs = $(rhs)"""

    return rhs, usable_ω
end

function λsp(χr::Array{Float64,2}, iωn::Array{Complex{Float64},1}, EKin::Float64,
                            rhs::Float64, kGrid::T1, mP::ModelParameters) where T1 <: ReducedKGrid
    #tmp = Array{Float64,2}(undef, size(χr)...)
    f(λint) = sum_freq(subtract_tail(kintegrate(kGrid, χ_λ(χr, λint), dim=2)[:,1],EKin, iωn), [1], Naive(), mP.β, corr=-EKin*mP.β^2/12)[1] - rhs
    df(λint) = sum_freq(kintegrate(kGrid, -χ_λ(χr, λint) .^ 2, dim=2)[:,1], [1], Naive(), mP.β)[1]

    nh    = ceil(Int64, size(χr,1)/2)
    χ_min    = -minimum(1 ./ χr[nh,:])
    λsp_old = newton_right(χr, f, df, χ_min)
    return λsp_old
end

function calc_λsp_correction(χ_in::SharedArray{Complex{Float64},2}, usable_ω::AbstractArray{Int64},
                            searchInterval::AbstractArray{Float64,1}, EKin::Float64,
                            rhs::Float64, kGrid::T1, mP::ModelParameters, sP::SimulationParameters) where T1 <: ReducedKGrid
    χr    = real.(χ_in[usable_ω,:])
    iωn = 1im .* 2 .* (-sP.n_iω:sP.n_iω)[usable_ω] .* π ./ mP.β
    f(λint) = sum_freq(subtract_tail(kintegrate(kGrid, χ_λ(χr, λint), dim=2)[:,1],EKin, iωn), [1], Naive(), mP.β, corr=-EKin*mP.β^2/12)[1] - rhs
    df(λint) = sum_freq(kintegrate(kGrid, -χ_λ(χr, λint) .^ 2, dim=2)[:,1], [1], Naive(), mP.β)[1]

    #TODO: new method needs testing
     X = @interval(0.02,20.0)
     tol = maximum([1e-7, (1e-4)*abs(searchInterval[2] - searchInterval[1])])
     #r = roots(f, df, X, Newton, tol)
     r2 = find_zeros(f, 0.00, 20.0, verbose=true)
     @info "Method 2 root:" r2
     #ln = empty(r2) ? 0.0 : mid(maximum(filter(x->!isempty(x),interval.(r))))
     #@info "Interval" ln

    nh    = ceil(Int64, size(χr,1)/2)
    χ_min    = -minimum(1 ./ χr[nh,:])
    λsp_old = newton_right(χr, f, df, χ_min)
    # if isempty(r)
    #   @warn "   ---> WARNING: no lambda roots with new method found!!!"
    # end
    # λsp = mid(maximum(filter(x->!isempty(x),interval.(r))))
    # if !isempty(r) && abs2(λsp_old - λsp) > 1e-5
    #    @warn "   ---> WARNING: old and new λ not matching!!! $(λsp_old) != $(λsp)"
    # end
    @info "Found λsp " λsp_old
    return λsp_old, new_χλ(χ_in, λsp_old, sP)
end

function extended_λ(nlQ_sp::NonLocalQuantities, nlQ_ch::NonLocalQuantities, bubble::BubbleT,
        Gνω::AbstractArray{Complex{Float64},2}, FUpDo::AbstractArray{Complex{Float64},3},
        Σ_loc_pos::AbstractArray{Complex{Float64},1}, Σ_ladderLoc::AbstractArray{Complex{Float64},1},kGrid::ReducedKGrid, tail_coeffs_upup, tail_coeffs_updo, mP::ModelParameters, sP::SimulationParameters; λsp_guess=0.1)
    # --- prepare auxiliary vars ---
    νZero = sP.n_iν
    ωZero = sP.n_iω
    ωindices = intersect(nlQ_sp.usable_ω, nlQ_ch.usable_ω)
    iωn = 1im .* 2 .* (-sP.n_iω:sP.n_iω)[ωindices] .* π ./ mP.β
    sh_f = get_sum_helper(2*sP.n_iν, sP, :f)
    sh_b = Naive() #get_sum_helper(length(ωindices), sP, :b)

    tmp = SharedArray{Complex{Float64},3}(length(ωindices), size(bubble,2), size(bubble,3))
    Σ_ladder_ω = SharedArray{Complex{Float64},3}(length(ωindices), size(bubble,2), sP.n_iν-sP.shift*(trunc(Int,sP.n_iω/2) + 1))
    χupdo_ω = SharedArray{eltype(nlQ_sp.χ),1}(length(ωindices))
    χupup_ω = SharedArray{eltype(nlQ_sp.χ),1}(length(ωindices))
    χsp_λ = SharedArray{eltype(nlQ_sp.χ),2}(length(ωindices),size(nlQ_sp.χ,2))
    χch_λ = SharedArray{eltype(nlQ_sp.χ),2}(length(ωindices),size(nlQ_ch.χ,2))

    Σ_hartree = mP.n * mP.U/2
    E_pot_tail_c = (mP.U^2 * 0.5 * mP.n * (1-0.5*mP.n) .+ Σ_hartree .* (kGrid.ϵkGrid .+ Σ_hartree .- mP.μ))
    E_pot_tail = E_pot_tail_c' ./ (iν_array(mP.β, 0:sP.n_iν-1) .^ 2)
    E_pot_tail_inv = sum((mP.β/2)  .* [Σ_hartree .* ones(size(kGrid.ϵkGrid)), (-mP.β/2) .* E_pot_tail_c])

    Σ_internal!(tmp, ωindices, bubble, FUpDo, sh_f)

    function cond_both!(F, λ)
        χsp_λ = SharedArray(χ_λ(nlQ_sp.χ, λ[1]))
        χch_λ = SharedArray(χ_λ(nlQ_ch.χ, λ[2]))
        calc_Σ_ω!(Σ_ladder_ω, ωindices, ωZero, νZero, sP.shift, χsp_λ, χch_λ, nlQ_sp.γ, nlQ_ch.γ, Gνω, tmp, mP.U, kGrid, sP)
        Σ_new = permutedims(mP.U .* sum_freq(Σ_ladder_ω, [1], Naive(), mP.β)[1,:,:], [2, 1])
        Σ_λ = Σ_new .- Σ_ladderLoc[1:size(Σ_new,1)] .+ Σ_loc_pos[1:size(Σ_new,1)]

        G_λ = G_from_Σ(Σ_λ .+ Σ_hartree, kGrid.ϵkGrid, 0:size(Σ_new, 1)-1, mP)
        E_pot_DGA = calc_E_pot(kGrid, flatten_2D(G_λ), Σ_λ, view(E_pot_tail,1:size(Σ_new, 1),:), E_pot_tail_inv)

        for (wi,w) in enumerate(ωindices)
            χupup_ω[wi] = kintegrate(kGrid, χch_λ[w,:] .+ χsp_λ[w,:])[1] ./ 2
            χupdo_ω[wi] = kintegrate(kGrid, χch_λ[w,:] .- χsp_λ[w,:])[1] ./ 2
        end

        χch_ω_sub = subtract_tail(χupup_ω, tail_coeffs_upup[3], iωn)
        lhs_c1 = real(sum_freq(χupup_ω, [1], sh_b, mP.β, corr=-tail_coeffs_upup[3]*mP.β^2/12)[1])
        lhs_c2 = real(sum_freq(χupdo_ω, [1], sh_b, mP.β)[1])

        rhs_c1 = mP.n/2 * (1 - mP.n/2)
        rhs_c2 = E_pot_DGA/mP.U + (mP.n/2) * (mP.n/2)
        F[1] = lhs_c1 - rhs_c1
        F[2] = lhs_c2 - rhs_c2
        println("-> $(F[1]), $(F[2])")
    end
    return nlsolve(cond_both!, [ λsp_guess; 0.0]).zero
end

function λsp_correction_search_int(χr::AbstractArray{Float64,2}, kGrid::ReducedKGrid, mP::ModelParameters; init=nothing, init_prct::Float64 = 0.1)
    nh    = ceil(Int64, size(χr,1)/2)
    χ_min    = -minimum(1 ./ χr[nh,:])
    int = if init === nothing
        rval = χ_min > 0 ? χ_min + 10/kGrid.Ns + 20/mP.β : minimum([20*abs(χ_min), χ_min + 10/kGrid.Ns + 20/mP.β])
        [χ_min, rval]
    else
        [init - init_prct*abs(init), init + init_prct*abs(init)]
    end
    @info "found " χ_min ". Looking for roots in intervall $(int)"
    return int
end

function λ_correction!(impQ_sp, impQ_ch, FUpDo, Σ_loc_pos, Σ_ladderLoc, nlQ_sp::NonLocalQuantities, nlQ_ch::NonLocalQuantities,
                      bubble::BubbleT, Gνω::SharedArray{Complex{Float64},2},
                      kGrid::ReducedKGrid,
                      mP::ModelParameters, sP::SimulationParameters; init_sp=nothing, init_spch=nothing)
    @info "Computing λsp corrected χsp, using " sP.χFillType " as fill value outside usable ω range."
    rhs,usable_ω_λc = calc_λsp_rhs_usable(impQ_sp, impQ_ch, nlQ_sp, nlQ_ch, kGrid, mP, sP)
    searchInterval_sp = λsp_correction_search_int(real.(nlQ_sp.χ[usable_ω_λc,:]), kGrid, mP, init=init_sp)
    searchInterval_spch = [-Inf, Inf]
    λsp,χsp_λ = calc_λsp_correction(nlQ_sp.χ, usable_ω_λc, searchInterval_sp, impQ_sp.tailCoeffs[3] , rhs, kGrid, mP, sP)
    #@info "Computing λsp corrected χsp, using " sP.χFillType " as fill value outside usable ω range."
    λ_new = extended_λ(nlQ_sp, nlQ_ch, bubble, Gνω, FUpDo, Σ_loc_pos, Σ_ladderLoc, kGrid, impQ_sp.tailCoeffs, [0.0,0.0,0.0,0.0,0.0], mP, sP; λsp_guess=λsp - 0.001)
    #λ_new = [0.0, 0.0]
    if sP.λc_type == :sp
        nlQ_sp.χ = χsp_λ
        nlQ_sp.λ = λsp
    elseif sP.λc_type == :sp_ch
        nlQ_sp.χ = SharedArray(χ_λ(nlQ_sp.χ, λ_new[1]))
        nlQ_sp.λ = λ_new[1]
        nlQ_ch.χ = SharedArray(χ_λ(nlQ_ch.χ, λ_new[2]))
        nlQ_ch.λ = λ_new[2]
    end
    @info "new lambda correction: λsp=$(λ_new[1]) and λch=$(λ_new[2])"
    return λsp, λ_new, searchInterval_sp, searchInterval_spch
end

function newton_right(χr::Array{Float64,2}, f::Function, df::Function,
                            start::Float64; nsteps=5000, atol=1e-11)
    done = false
    δ = 0.1
    x0 = start + δ
    xi = x0
    i = 1
    while !done
        fi = f(xi)
        dfi = df(xi)
        xlast = xi
        xi = x0 - fi / dfi
        (abs2(xi-x0) < atol) && break
        if xi < x0               # only ever search to the right!
            δ  = δ/2.0
            x0  = start + δ      # reset with smaller delta
            xi = x0
        else
            x0 = xi
        end
        (i >= nsteps ) && (done = true)
        i += 1
    end
    return xi
end


text:
using Yao
using LuxurySparse
using LinearAlgebra
using TropicalYao
using Viznet, Compose

export solve, SquareLattice, ChimeraLattice
export sgbonds, sgvertices

export Gh, Gvb, Ghb, G16, Gcut, Gcp, Greset
export vertextensor, bondtensor

function _init_reg(::Type{T}, L::Int, usecuda::Val{:false}) where T
    ArrayReg(ones(T, 1<<L))
end

function solve(sg::AbstractSpinglass{LT,T}; usecuda=false) where {LT,T}
    solve(Tropical{T}, sg; usecuda=usecuda)
end

function solve_and_count(sg::Spinglass{LT,T}; usecuda=false) where {LT,T}
    solve(CountingTropical{T}, sg; usecuda=usecuda)
end

struct SpinglassOptConfig{LT,T}
    sg::Spinglass{LT,T}
    eng::T
    grad_Js::Vector{T}
    grad_hs::Vector{T}
end

export tofile, fromfile
function fromfile(::Type{T}, filename::String, lattice) where T
    nv = length(lattice)
    ne = length(sgbonds(lattice))
    Js = zeros(T, ne)
    hs = zeros(T, nv)
    Js_grad = zeros(T, ne)
    hs_grad = zeros(T, nv)
    open(filename, "r") do f
        for i=1:length(Js)
            l = strip(readline(f))
            j, g = split(l, ' ')
            Js[i] = parse(T, j)
            Js_grad[i] = parse(T, g)
        end
        for i=1:length(hs)
            l = strip(readline(f))
            h, g = split(l, ' ')
            hs[i] = parse(T, h)
            hs_grad[i] = parse(T, g)
        end
    end
    SpinglassOptConfig(Spinglass(lattice, Js, hs), T(0), Js_grad, hs_grad)
end

function tofile(filename::String, res::SpinglassOptConfig)
    sg = res.sg
    open(filename, "w") do f
        for (J, g) in zip(sg.Js, res.grad_Js)
            write(f, "$J $g\n")
        end
        for (h, g) in zip(sg.hs, res.grad_hs)
            write(f, "$h $g\n")
        end
    end
    return res
end

tofile(filename::String) = res->tofile(filename, res)

include("square.jl")
include("chimera.jl")
include("second_neighbor.jl")
include("cylinder.jl")
include("cubic.jl")

function assign_Js(lt::Viznet.AbstractLattice, g::AbstractVector{T}) where T
    grid = zeros(length(lt))
    grid[1] = 1
    configs = length(lt)
    bonds = sgbonds(lt)
    remain_bonds = collect(1:length(bonds))
    for i=1:10
        nrem = length(remain_bonds)
        for (i_, k) in enumerate(Base.Iterators.reverse(copy(remain_bonds)))
            i, j = bonds[k]
            res = assign_one!(grid, i, j, g[k])
            if res
                deleteat!(remain_bonds, nrem-i_+1)
            end
        end
        isempty(remain_bonds) && break
    end
    if !isempty(remain_bonds)
        error("bonds $remain_bonds can not be set!")
    end
    return grid
end

function assign_Js_hs(lt::Viznet.AbstractLattice, grad_Js::AbstractVector{T}, grad_hs) where T
    grid = zeros(length(lt))
    for (v, hg) in zip(sgvertices(lt), grad_hs)
        grid[v] = hg
    end
    for ((i, j), Jg) in zip(sgbonds(lt), grad_Js)
        assign_one!(grid, i, j, Jg)
    end
    return grid
end

function assign_one!(grid, x, y, g)
    if grid[x] == 0 && grid[y] == 0
        return false
    elseif grid[x] == 0
        grid[x] = sign(g)*grid[y]
    elseif grid[y] == 0
        grid[y] = sign(g)*grid[x]
    else
        if grid[y] != sign(g)*grid[x]
            @warn "Grid assign σ($y) = $(grid[y]) and σ($x) = $(grid[x]) inconsistent with bond gradient $(g)!"
        end
    end
    return true
end


max_stars_repo_path: test/1 Algebraic functions/1.1 Binomial products/1.1.1 Linear/1.1.1.7 P(x) (a+b x)^m (c+d x)^n (e+f x)^p (g+h x)^q.jl


text:
using Test
using SafeTestsets

@time @safetestset "ESN constructors" begin include("constructors/test_esn_constructors.jl") end
@time @safetestset "HESN constructors" begin include("constructors/test_hesn_constructors.jl") end
@time @safetestset "DAFESN constructors" begin include("constructors/test_dafesn_constructors.jl") end
@time @safetestset "ESN input layers" begin include("fixed_layers/test_esn_input_layers.jl") end
@time @safetestset "ESN reservoirs" begin include("fixed_layers/test_esn_reservoirs.jl") end
@time @safetestset "Non Linear Algorithms" begin include("extras/test_nla.jl") end
@time @safetestset "Extended states" begin include("extras/test_extended_states.jl") end
@time @safetestset "MLJ Linear Models for ESN" begin include("training/test_mlj_lm.jl") end
@time @safetestset "ESGP" begin include("training/test_esgp.jl") end
@time @safetestset "RECA_discrete constructors" begin include("constructors/test_recadiscrete_constructors.jl") end
@time @safetestset "RECA gol constructors" begin include("constructors/test_recagol_constructors.jl") end
@time @safetestset "reca discrete predict" begin include("training/test_recadiscrete.jl") end
@time @safetestset "reca gol predict" begin include("training/test_recagol.jl") end
@time @safetestset "RMM constructors" begin include("constructors/test_rmm_constructors.jl") end
@time @safetestset "GRUESN constructors" begin include("constructors/test_gruesn_constructors.jl") end


text:
import ..InterpolationAndGrids: interpolate_and_bin

"""
    resample(udata::AbstractUncertainIndexValueDataset, regularization_scheme::InterpolateAndBin{Linear})

Draw a single realisation of `udata` and interpolate-and-bin the data according to the
provided regularization scheme. Assumes points in `udata` are independent and sorts the draw
according to the index values before interpolating. See also [`InterpolateAndBin`](@ref).

## Example

```julia
npts = 50
y = rand(npts)

N = Normal(0, 1)

for t in 3:npts
    y[t,1] = 0.7*y[t-1,1] - 0.35*y[t-2,1] + rand(N)
end

# Assume data are unevenly spaced
time = sample(1.0:npts*5, npts, ordered = true, replace = false)

# Assign some uncertainties to both time indices and values and gather
# in an UncertainIndexValueDataset
utime = UncertainValue.(Normal.(time, 2))
uy = UncertainValue.(Normal.(y, 0.1))
udata = UncertainIndexValueDataset(utime, uy)

# Interpolation-and-binning scheme. First interpolate to a very fine grid,
# then gather the points falling in each of the coarser bins and summarise
# each bin using the mean of the points in each bin.
left_bin_edges = 0:10:npts*5
r = InterpolateAndBin(mean, left_bin_edges, Linear(), 0:0.1:1000, Flat(OnGrid()))

# The binned time axis:
time_binned = left_bin_edges[1:end-1] .+ step(left_bin_edges)/2

# Get a set corresponding resampled (interpolated+binned) values
y_binned = resample(udata, r)

# Plot some interpolated+binned draws
time_binned = left_bin_edges[1:end-1] .+ step(left_bin_edges)/2

p = plot(xlabel = "time", ylabel = "value")
for i = 1:100
    plot!(time_binned, resample(udata, r), lw = 0.3, α = 0.2, ms = 0.1, c = :red,
        marker = stroke(0.1), label = "")
end
plot!(time, y, c = :black, lw = 1, ms = 2, marker = stroke(2.0, :black), label = "")
plot!(udata, c = :black, lw = 1, ms = 2, marker = stroke(0.1, :black), [0.05, 0.95], [0.05, 0.95])
vline!(left_bin_edges, c = :black, α = 0.3, lw = 0.3, label = "")
```
"""
function resample(udata::AbstractUncertainIndexValueDataset, regularization_scheme::InterpolateAndBin{Linear})
    inds, vals = resample(udata)

    sortidxs = sortperm(inds)

    r = regularization_scheme
    vals_binned = interpolate_and_bin(r.f, r.left_bin_edges, inds[sortidxs], vals[sortidxs],
        r.intp, r.intp_grid, r.extrapolation_bc)

    return vals_binned
end

export resample


text:
# -----------------------------------------------------------------------
# Functions that provide the control of alls the stepper motors
# -----------------------------------------------------------------------
# Required packages
using StepperControl, Optim
using BSON: @load
include("cylinder_structure.jl")

"""
    step2coord_x(;sprx=2048, R=80.0, xorig=-41.5)

## Description
Generates a  function which converts steps into displacement in x direction
## Arguments
- sprx: steps per revolution in x system
- R: bell crank radius
- xorig: horizontal distance from the center of rotation to the origin
## Example
```jldoctest
julia> f = step2coord_x();

julia> f(20)
4.323700717312371
```
"""
function step2coord_x(;sprx=2048, R=80.0, xorig=-41.5)
    θ0 = asin(xorig/R)
    x(v) = R*(sin(θ0 + 2*π*v/sprx) - sin(θ0))
    return x
end

"""
    coord2step_x(;sprx=2048, R=80.0, xorig=-41.5)

## Description
Generates a  function which converts displacement in x direction into steps
## Arguments
- sprx: steps per revolution in x system
- R: bell crank radius
- xorig: horizontal distance from the center of rotation to the origin
## Example
```jldoctest
julia> f = coord2step_x();

julia> f(4.32)
20
```
"""
function coord2step_x(;sprx=2048, R=80.0, xorig=-41.5)
    θ0 = asin(xorig/R)
    x(v) = Int(round( (sprx/(2π)) * (asin(v/R + sin(θ0)) - θ0) ))
    return x
end

"""
    step2coord_y(;sprx=2048, spry=2048, R=80.0, r=12.5, xorig=-41.5)

## Description
Generates a  function which converts steps into displacement in y direction
## Arguments
- sprx: steps per revolution in x system
- spry: steps per revolution in y system
- R: bell crank radius
- r: gear radius
- xorig: horizontal distance from the center of rotation to the origin
## Example
```jldoctest
julia> f = step2coord_y();

julia> f([20,100])
1.5861345582862714
```
"""
function step2coord_y(;sprx=2048, spry=2048, R=80.0, r=12.5, xorig=-41.5)
    θ0 = asin(xorig/R)
    y(v) = 2π*v[2]*r/spry - R*(cos(θ0 + 2*π*v[1]/sprx) - cos(θ0))
    return y
end

"""
    coord2step_y(;sprx=2048, spry=2048, R=80.0, r=12.5, xorig=-41.5)

## Description
Generates a  function which converts displacement in y direction into steps
## Arguments
- sprx: steps per revolution in x system
- spry: steps per revolution in y system
- R: bell crank radius
- r: gear radius
- xorig: horizontal distance from the center of rotation to the origin
## Example
```jldoctest
julia> f = coord2step_y();

julia> f(4.32, 1.58)
100
```
"""
function coord2step_y(;sprx=2048, spry=2048, R=80.0, r=12.5, xorig=-41.5)
    θ0 = asin(xorig/R)
    y(v) = Int(round( (spry/(2π*r)) * (v[2] + R*(cos(asin(v[1]/R + sin(θ0))) - cos(θ0))) ))
    return y
end

"""
    beststrategy(dev::StepperSystem, posnew::OneCylType)

## Description
Find best moving strategy to avoid collisions
## Arguments
- dev: stepper system from StepperControl library
- posnew: new CylType position
"""
beststrategy(dev::StepperSystem, posnew::OneCylType) = 1:2

function beststrategy(dev::StepperSystem, posnew::TwoCylType)

    posold = getpos(dev)

    if posnew[1] >= posold[1]
        if posnew[1] >= 55.0
            strategy = [3,4,2,1]
        else
            strategy = [3,4,1,2]
        end
    else
        if posnew[1] <= 41.0
            strategy = [1,2,4,3]
        else
            strategy = 1:4
        end
    end

    return strategy
end

function beststrategy(dev::StepperSystem, posnew::ThreeCylType)

    posold = getpos(dev)

    if posnew[1] >= posold[1]
        if posnew[1]+posnew[5] <= 75.0
            strategy = [5,6,3,4,2,1]
        else
            strategy = [5,6,4,3,2,1]
        end
    else
        if posnew[1] >= 42.0
            strategy = [1,2,3,4,6,5]
        else
            strategy = [1,2,4,3,6,5]
        end
    end

    return strategy
end

"""
    open_onecylinder(sprx=2121, spry=2018; R=80.0, r=12.5, xorig=-41.5, testnocon=false)

## Description
Function to get connection and configure stepper system for one cylinder
## Arguments
- sprx: steps per revolution in x direction
- spry: steps per revolution in y direction
- R: bell crank length
- r: gear radius
- xorig: horizontal distance from de center of the rotor to the origin
- testnocon: if true then no connection is given (for testing purposes)
"""
function open_onecylinder(sprx=2099, spry=2007; R=80.0, r=12.5, xorig=-41.5, testnocon=false)

    # getting connection
    if testnocon
        dev = stepper_open(2, testnocon=true)
    else
        dev = stepper_open(2)
    end

    # transformation functions
    s2cx = step2coord_x(sprx=sprx, R=R, xorig=xorig)
    s2cy = step2coord_y(sprx=sprx, spry=spry, R=R, r=r, xorig=xorig)
    c2sx = coord2step_x(sprx=sprx, R=R, xorig=xorig)
    c2sy = coord2step_y(sprx=sprx, spry=spry, R=R, r=r, xorig=xorig)

    stepper_config!(dev, motorID=["x1","y1"], step2coord=[s2cx,s2cy], coord2step=[c2sx,c2sy], depend=[1,[1,2]])

    return dev
end

"""
    open_cylinderNN(NNmodelfile::AbstractString; xorig=-59.0, motorID=["x2","y2"], testnocon=false)

## Description
Function to get connection and configure stepper system for one cylinder considering neural net models
## Arguments
- NN: neural network
- xorig: horizontal distance from de center of the rotor to the origin
- motorID: stepper motor ids
- testnocon: if true then no connection is given (for testing purposes)
"""
function open_cylinderNN(NN; xorig=-59.0, motorID=["x2","y2"], testnocon=false)

    # getting connection
    if testnocon
        dev = stepper_open(2, testnocon=true)
    else
        dev = stepper_open(2)
    end

    # ideal models
    θ0 = asin(xorig/80)
    s2cx(sx) = 80*(sin(θ0+ 2π*sx/2094) - sin(θ0))
    s2cy(s) = 2π*12.5*s[2]/1994 - 80*(cos(θ0 + 2*π*s[1]/2094) - cos(θ0)) + NN(s[1])
    c2sx(x) = Int(round( (2094/(2π))*(asin(x/80+sin(θ0)) - θ0) ))

    approx_c2sy(c) = (1994/(2π*12.5))*(c[2] + 80*(cos(asin(c[1]/80+sin(θ0)))-cos(θ0)))
    function c2sy(c)
        sx = c2sx(c[1])
        minval(sy) = sum(abs2, c[2] - s2cy([sx,sy]))
        r = optimize(v -> minval(first(v)), [approx_c2sy(c)], LBFGS())
        return Int(round(r.minimizer[1]))
    end

    stepper_config!(dev, motorID=motorID, step2coord=[s2cx,s2cy], coord2step=[c2sx,c2sy], depend=[1,[1,2]])

    return dev
end

"""
    open_twocylinder(; testnocon=false)

"""
function open_twocylinder(; testnocon=false)

    # getting connection
    if testnocon
        dev = stepper_open(4, testnocon=true)
    else
        dev = stepper_open(4)
    end

    # transformation functions cylinder 1
    θ0₁ = asin(-41.5/80.0)
    s2cx1(sx) = 80*(sin(θ0₁+ 2π*sx/2099) - sin(θ0₁))
    s2cy1(s) = 2π*12.5*s[2]/2007 - 80*(cos(θ0₁ + 2*π*s[1]/2099) - cos(θ0₁))
    c2sx1(x) = Int(round( (2099/(2π))*(asin(x/80+sin(θ0₁)) - θ0₁) ))
    c2sy1(c) = Int(round( (2007/(2π*12.5))*(c[2] + 80*(cos(asin(c[1]/80+sin(θ0₁)))-cos(θ0₁))) ))

    # transformation functions cylinder 2
    θ0₂ = asin(-59.0/80.0)
    s2cx2(sx) = 80*(sin(θ0₂+ 2π*sx/2094) - sin(θ0₂))
    s2cy2(s) = 2π*12.5*s[2]/1994 - 80*(cos(θ0₂ + 2*π*s[1]/2094) - cos(θ0₂)) + NNyx2(s[1])
    c2sx2(x) = Int(round( (2094/(2π))*(asin(x/80+sin(θ0₂)) - θ0₂) ))
    approx_c2sy2(c) = (1994/(2π*12.5))*(c[2] + 80*(cos(asin(c[1]/80+sin(θ0₂)))-cos(θ0₂)))
    function c2sy2(c)
        sx = c2sx2(c[1])
        minval(sy) = sum(abs2, c[2] - s2cy2([sx,sy]))
        r = optimize(v -> minval(first(v)), [approx_c2sy2(c)], LBFGS())
        return Int(round(r.minimizer[1]))
    end

    stepper_config!(dev, motorID=["x1","y1","x2","y2"], step2coord=[s2cx1,s2cy1,s2cx2,s2cy2],
                    coord2step=[c2sx1,c2sy1,c2sx2,c2sy2], depend=[1,[1,2],3,[3,4]])

    return dev
end

"""
    open_threecylinder(; testnocon=false)

"""
function open_threecylinder(; testnocon=false)

    # getting connection
    if testnocon
        dev = stepper_open(6, testnocon=true)
    else
        dev = stepper_open(6)
    end

    # transformation functions cylinder 1
    θ0₁ = asin(-41.5/80.0)
    s2cx1(sx) = 80*(sin(θ0₁+ 2π*sx/2099) - sin(θ0₁))
    s2cy1(s) = 2π*12.5*s[2]/2009 - 80*(cos(θ0₁ + 2*π*s[1]/2099) - cos(θ0₁)) + NNyx1(s[1])
    c2sx1(x) = Int(round( (2099/(2π))*(asin(x/80+sin(θ0₁)) - θ0₁) ))
    #c2sy1(c) = Int(round( (2009/(2π*12.5))*(c[2] + 80*(cos(asin(c[1]/80+sin(θ0₁)))-cos(θ0₁))) ))
    approx_c2sy1(c) = (2009/(2π*12.5))*(c[2] + 80*(cos(asin(c[1]/80+sin(θ0₂)))-cos(θ0₂)))
    function c2sy1(c)
        sx = c2sx1(c[1])
        minval(sy) = sum(abs2, c[2] - s2cy1([sx,sy]))
        r = optimize(v -> minval(first(v)), [approx_c2sy1(c)], LBFGS())
        return Int(round(r.minimizer[1]))
    end

    # transformation functions cylinder 2
    θ0₂ = asin(-59.0/80.0)
    s2cx2(sx) = 80*(sin(θ0₂+ 2π*sx/2176) - sin(θ0₂))
    s2cy2(s) = 2π*12.5*s[2]/1985 - 80*(cos(θ0₂ + 2*π*s[1]/2176) - cos(θ0₂)) + NNyx2(s[1])
    c2sx2(x) = Int(round( (2176/(2π))*(asin(x/80+sin(θ0₂)) - θ0₂) ))
    approx_c2sy2(c) = (1985/(2π*12.5))*(c[2] + 80*(cos(asin(c[1]/80+sin(θ0₂)))-cos(θ0₂)))
    function c2sy2(c)
        sx = c2sx2(c[1])
        minval(sy) = sum(abs2, c[2] - s2cy2([sx,sy]))
        r = optimize(v -> minval(first(v)), [approx_c2sy2(c)], LBFGS())
        return Int(round(r.minimizer[1]))
    end

    # transformation functions cylinder 3
    θ0₃ = asin(-75.0/80.0)
    s2cx3(sx) = 80*(sin(θ0₃+ 2π*sx/2097) - sin(θ0₃))
    s2cy3(s) = 2π*12.5*s[2]/2041 - 80*(cos(θ0₃ + 2*π*s[1]/2097) - cos(θ0₃)) + NNyx3(s[1])
    c2sx3(x) = Int(round( (2097/(2π))*(asin(x/80+sin(θ0₃)) - θ0₃) ))
    approx_c2sy3(c) = (2041/(2π*12.5))*(c[2] + 80*(cos(asin(c[1]/80+sin(θ0₃)))-cos(θ0₃)))
    function c2sy3(c)
        sx = c2sx3(c[1])
        minval(sy) = sum(abs2, c[2] - s2cy3([sx,sy]))
        r = optimize(v -> minval(first(v)), [approx_c2sy3(c)], LBFGS())
        return Int(round(r.minimizer[1]))
    end

    stepper_config!(dev, motorID=["x1","y1","x2","y2","x3","y3"], step2coord=[s2cx1,s2cy1,s2cx2,s2cy2,s2cx3,s2cy3],
                    coord2step=[c2sx1,c2sy1,c2sx2,c2sy2,c2sx3,c2sy3], depend=[1,[1,2],3,[3,4],5,[5,6]])

    return dev
end

"""
    movecyl!(dev::StepperSystem, posnew::OneCylType; relat=false, strategy=1:2)

## Description
Function to move the cylinders system
## Arguments
- dev: stepper system from StepperControl library
- posnew: new CylType position
- strategy: stepper motor trigger order. If missing, then beststrategy() is called.
- relat: if true then the relative movement is performed. Default=false.
- method: method provided by the StepperControl library
"""
function movecyl!(dev::StepperSystem, posnew::CylType; strategy=missing,
                    relat=false, method=StepperControl.manhattan_msg)

    coords = convert(Array, posnew)

    if typeof(posnew) == TwoCylType
        coords[3] += coords[1]
        coords[4] += coords[2]
    elseif typeof(posnew) == ThreeCylType
        coords[3] += coords[1]
        coords[4] += coords[2]
        coords[5] += coords[3]
        coords[6] += coords[4]
    end

    if ismissing(strategy)
        strategy = beststrategy(dev,posnew)
    end

    msg = stepper_move!(dev, coords, relat=relat, order=strategy, method=method)

    println(msg)
end


text:
using OrdinaryDiffEq, Test
f = (u,p,t) -> u*p[1]
prob = ODEProblem(f,1.01,(0.0,1.0))
@test_throws ErrorException solve(prob,Tsit5())


text:
# get bare propagator
function get_G_bare(
    Λ :: Float64,
    w :: Float64
    ) :: Float64

    val = 0.0

    if abs(w) > 1e-8
        val = -expm1(-w^2 / Λ^2) / w
    end

    return val
end

# get dressed propagator
function get_G(
    Λ :: Float64,
    w :: Float64,
    m :: Mesh,
    a :: Action
    ) :: Float64

    val = 0.0

    if abs(w) > 1e-8
        Σ     = get_Σ(w, m, a)
        G0    = get_G_bare(Λ, w)
        denom = 0.0

        # ensure causality
        if w * Σ >= 0.0
            denom = 1.0 / G0 + Σ
        else
            denom = 1.0 / G0
        end

        val = 1.0 / denom
    end

    return val
end

# get single scale propagator
function get_S(
    Λ :: Float64,
    w :: Float64,
    m :: Mesh,
    a :: Action
    ) :: Float64

    val = 0.0

    if abs(w) > 1e-8
        G   = get_G(Λ, w, m, a)
        G0  = get_G_bare(Λ, w)
        val = (G / G0)^2 * exp(-w^2 / Λ^2) * 2.0 * w / Λ^3
    end

    return val
end

# get differentiated propagator bubble
function get_propagator_kat(
    Λ  :: Float64,
    w1 :: Float64,
    w2 :: Float64,
    m  :: Mesh,
    a  :: Action,
    da :: Action
    )  :: Float64

    dΣ  = get_Σ(w1, m, da)
    S1  = get_S(Λ, w1, m, a)
    G1  = get_G(Λ, w1, m, a)
    G2  = get_G(Λ, w2, m, a)
    val = (S1 + G1^2 * dΣ) * G2 / (2.0 * pi)

    return val
end

# get undifferentiated propagator bubble
function get_propagator(
    Λ  :: Float64,
    w1 :: Float64,
    w2 :: Float64,
    m  :: Mesh,
    a  :: Action
    )  :: Float64

    G1  = get_G(Λ, w1, m, a)
    G2  = get_G(Λ, w2, m, a)
    val = G1 * G2 / (2.0 * pi)

    return val
end

# compute boundary corrections for vertex integrals
function compute_corrs_kat!(
    Λ     :: Float64,
    m     :: Mesh,
    a     :: Action,
    da    :: Action,
    corrs :: Array{Float64, 3},
    Γ_tol :: NTuple{2, Float64}
    )     :: Nothing

    # reset buffer
    @turbo corrs .= 0.0

    # compute boundary corrections
    @sync for i in size(corrs, 3)
        Threads.@spawn begin
            # compute boundary corrections for s channel
            s_propagator  = v -> get_propagator_kat(Λ, v + 0.5 * m.Ωs[i], 0.5 * m.Ωs[i] - v, m, a, da) + get_propagator_kat(Λ, 0.5 * m.Ωs[i] - v, v + 0.5 * m.Ωs[i], m, a, da)
            s_val         = m.Ωs[end] + m.νs[end]
            s_bound_minus = s_propagator(-s_val)
            s_bound_plus  = s_propagator( s_val)

            if abs(s_bound_minus) > 1e-8
                corrs[1, 1, i] = quadgk(s_propagator, Inf, -s_val, atol = Γ_tol[1], rtol = Γ_tol[2])[1] / s_bound_minus
            end

            if abs(s_bound_plus) > 1e-8
                corrs[2, 1, i] = quadgk(s_propagator, s_val, Inf, atol = Γ_tol[1], rtol = Γ_tol[2])[1] / s_bound_plus
            end

            # compute boundary corrections for t channel
            t_propagator  = v -> get_propagator_kat(Λ, v + 0.5 * m.Ωt[i], v - 0.5 * m.Ωt[i], m, a, da) + get_propagator_kat(Λ, v - 0.5 * m.Ωt[i], v + 0.5 * m.Ωt[i], m, a, da)
            t_val         = m.Ωt[end] + m.νt[end]
            t_bound_minus = t_propagator(-t_val)
            t_bound_plus  = t_propagator( t_val)

            if abs(t_bound_minus) > 1e-8
                corrs[1, 2, i] = quadgk(t_propagator, -Inf, -t_val, atol = Γ_tol[1], rtol = Γ_tol[2])[1] / t_bound_minus
            end

            if abs(t_bound_plus) > 1e-8
                corrs[2, 2, i] = quadgk(t_propagator, t_val, Inf, atol = Γ_tol[1], rtol = Γ_tol[2])[1] / t_bound_plus
            end

            # compute boundary corrections for u channel
            u_propagator  = v -> get_propagator_kat(Λ, v - 0.5 * m.Ωs[i], v + 0.5 * m.Ωs[i], m, a, da) + get_propagator_kat(Λ, v + 0.5 * m.Ωs[i], v - 0.5 * m.Ωs[i], m, a, da)
            u_val         = m.Ωs[end] + m.νs[end]
            u_bound_minus = u_propagator(-u_val)
            u_bound_plus  = u_propagator( u_val)

            if abs(u_bound_minus) > 1e-8
                corrs[1, 3, i] = quadgk(u_propagator, -Inf, -u_val, atol = Γ_tol[1], rtol = Γ_tol[2])[1] / u_bound_minus
            end

            if abs(u_bound_plus) > 1e-8
                corrs[2, 3, i] = quadgk(u_propagator, u_val, Inf, atol = Γ_tol[1], rtol = Γ_tol[2])[1] / u_bound_plus
            end
        end
    end

    return nothing
end

# compute boundary corrections for vertex integrals
function compute_corrs!(
    Λ     :: Float64,
    m     :: Mesh,
    a     :: Action,
    corrs :: Array{Float64, 3},
    Γ_tol :: NTuple{2, Float64}
    )     :: Nothing

    # reset buffer
    @turbo corrs .= 0.0

    # compute boundary corrections
    @sync for i in size(corrs, 3)
        Threads.@spawn begin
            # compute boundary corrections for s channel
            s_propagator  = v -> -get_propagator(Λ, v + 0.5 * m.Ωs[i], 0.5 * m.Ωs[i] - v, m, a)
            s_val         = m.Ωs[end] + m.νs[end]
            s_bound_minus = s_propagator(-s_val)
            s_bound_plus  = s_propagator( s_val)

            if abs(s_bound_minus) > 1e-8
                corrs[1, 1, i] = quadgk(s_propagator, -Inf, -s_val, atol = Γ_tol[1], rtol = Γ_tol[2])[1] / s_bound_minus
            end

            if abs(s_bound_plus) > 1e-8
                corrs[2, 1, i] = quadgk(s_propagator, s_val, Inf, atol = Γ_tol[1], rtol = Γ_tol[2])[1] / s_bound_plus
            end

            # compute boundary corrections for t channel
            t_propagator  = v -> -get_propagator(Λ, v + 0.5 * m.Ωt[i], v - 0.5 * m.Ωt[i], m, a)
            t_val         = m.Ωt[end] + m.νt[end]
            t_bound_minus = t_propagator(-t_val)
            t_bound_plus  = t_propagator( t_val)

            if abs(t_bound_minus) > 1e-8
                corrs[1, 2, i] = quadgk(t_propagator, -Inf, -t_val, atol = Γ_tol[1], rtol = Γ_tol[2])[1] / t_bound_minus
            end

            if abs(t_bound_plus) > 1e-8
                corrs[2, 2, i] = quadgk(t_propagator, t_val, Inf, atol = Γ_tol[1], rtol = Γ_tol[2])[1] / t_bound_plus
            end

            # compute boundary corrections for u channel
            u_propagator  = v -> -get_propagator(Λ, v - 0.5 * m.Ωs[i], v + 0.5 * m.Ωs[i], m, a)
            u_val         = m.Ωs[end] + m.νs[end]
            u_bound_minus = u_propagator(-u_val)
            u_bound_plus  = u_propagator( u_val)

            if abs(u_bound_minus) > 1e-8
                corrs[1, 3, i] = quadgk(u_propagator, -Inf, -u_val, atol = Γ_tol[1], rtol = Γ_tol[2])[1] / u_bound_minus
            end

            if abs(u_bound_plus) > 1e-8
                corrs[2, 3, i] = quadgk(u_propagator, u_val, Inf, atol = Γ_tol[1], rtol = Γ_tol[2])[1] / u_bound_plus
            end
        end
    end

    return nothing
end


text:
requiredTypes = [ "CareerStatus" ]

for reqType in requiredTypes
    if !isdefined( Symbol( uppercase( string( reqType[ 1 ] ) ) * reqType[ 2:end ] ) )
        include( joinpath( dirname( Base.source_path() ), "..", "types", reqType * ".jl" ) )
    end  # if !isdefined( Symbol( ...
end

function get_AcademicLevel(cs::CareerStatus)
    return cs.AcademicLevel
end

function set_AcademicLevel(cs::CareerStatus, Val::AbstractLevel)
    cs.AcademicLevel = Val
end

function get_Rank(cs::CareerStatus)
    return cs.Rank
end

function set_Rank(cs::CareerStatus, Val::AbstractLevel)
    cs.Rank = Val
end

function get_Assignment(cs::CareerStatus)
    return cs.Assignment
end

function set_Assignment(cs::CareerStatus, Val::AbstractLevel)
    cs.Assignment = Val
end

function get_Affiliation(cs::CareerStatus)
    return cs.Affiliation
end

function set_Affiliation(cs::CareerStatus, Val::AbstractLevel)
    cs.Affiliation = Val
end

function get_Job(cs::CareerStatus)
    return cs.Job
end

function set_Job(cs::CareerStatus, Val::AbstractLevel)
    cs.Job = Val
end

function get_SubJob(cs::CareerStatus)
    return cs.SubJob
end

function set_SubJob(cs::CareerStatus, Val::AbstractLevel)
    cs.SubJob = Val
end

function set_MinStay(st::CareerStatus, min::Int)
    st.MinStay = min
end

function set_MaxStay(st::CareerStatus, max::Int)
    st.MaxStay = max
end

function get_MinStay(st::CareerStatus)
    st.MinStay
end

function get_MaxStay(st::CareerStatus)
    st.MaxStay
end

function get_Name_Level(cs::CareerStatus, typ::DataType)
    if typ == AcademicLevel
        return cs.AcademicLevel.Name
    elseif typ == Rank
        return cs.Rank.Name
    elseif typ == Assignment
        return cs.Assignment.Name
    elseif typ == Affiliation
        return cs.Affiliation.Name
    elseif typ == Job
        return cs.Job.Name
    elseif typ == SubJob
        return cs.SubJob.Name
    end
end

function get_Level(cs::CareerStatus, typ::DataType)
    if typ == AcademicLevel
        return cs.AcademicLevel
    elseif typ == Rank
        return cs.Rank
    elseif typ == Assignment
        return cs.Assignment
    elseif typ == Affiliation
        return cs.Affiliation
    elseif typ == Job
        return cs.Job
    elseif typ == SubJob
        return cs.SubJob
    end
end

max_stars_repo_name: Omazaria/CareerPathBeDef

max_issues_repo_name: Omazaria/CareerPathBeDef

max_forks_repo_name: Omazaria/CareerPathBeDef


text:
using DataFrames, FixedEffectModels, StatsModels
N = 10000000
K = 100
id1 = Int.(rand(1:(N/K), N))
id2 = Int.(rand(1:K, N))
w = cos.(id1)

x1 = 5 * cos.(id1) + 5 * sin.(id2) + randn(N)
x2 =  cos.(id1) +  sin.(id2) + randn(N)
y= 3 .* x1 .+ 5 .* x2 .+ cos.(id1) .+ cos.(id2).^2 .+ randn(N)

df = DataFrame(id1 = categorical(id1), id2 = categorical(id2), x1 = x1, x2 = x2, w = w, y = y)
@time completecases(df, [:y, :x1, :x2, :id2])

@time ModelFrame(@formula(y~x1 + x2), df)

@time reg(df, @model(y ~ x1 + x2))
# 0.582029 seconds (852 allocations: 535.311 MiB, 18.28% gc time)
@time reg(df, @model(y ~ x1 + x2, vcov = cluster(id2)))
# 0.809652 seconds (1.55 k allocations: 649.821 MiB, 14.40% gc time)
@time reg(df, @model(y ~ x1 + x2, fe = id1))
# 1.655732 seconds (1.21 k allocations: 734.353 MiB, 16.88% gc time)
@time reg(df, @model(y ~ x1 + x2, fe = id1, vcov = cluster(id1)))
#  2.113248 seconds (499.66 k allocations: 811.751 MiB, 15.08% gc time)
@time reg(df, @model(y ~ x1 + x2, fe = id1 + id2))
# 3.553678 seconds (1.36 k allocations: 1005.101 MiB, 10.55% gc time))

# more regressors
df.x3 =  cos.(id1) + sin.(id2) + randn(N)
df.x4 =  cos.(id1) + sin.(id2) + randn(N)
df.x5 =  cos.(id1) + sin.(id2) + randn(N)
df.x6 =  cos.(id1) + sin.(id2) + randn(N)
df.x7 =  cos.(id1) + sin.(id2) + randn(N)
@time reg(df, @model(y ~ x1 + x2 + x3 + x4 + x5 + x6 + x7, subset = x3 .>= 0.5))

@time reg(df, @model(y ~ x1 + x2 + x3 + x4 + x5 + x6 + x7, fe = id1, subset = x3 .>= 0.5))
#  4.064132 seconds (2.24 k allocations: 965.500 MiB, 12.50% gc time)

df.id1 = categorical(mod.(1:size(df, 1), Ref(5)))
df.id4 = categorical(mod.(1:size(df, 1), Ref(4)))
@time reg(df, @model(y ~ id4, fe = id1), method = :lsmr)

# Benchmark Parallel
df.id3 = categorical(Int.(rand(1:15, N)))
df.x3 =  cos.(id1) + sin.(id2) + randn(N)
sort!(df, [:id1])
@time reg(df, @model(y ~ x1 + id3, fe = id1 + id2 + id2&x3, weights = w), method = :lsmr)
@time reg(df, @model(y ~ x1 + id3, fe = id1 + id2 + id2&x3, weights = w), method = :lsmr_parallel)
@time reg(df, @model(y ~ x1 + id3, fe = id1 + id2 + id2&x3, weights = w), method = :lsmr_threads)

#check that as fast as lm with no fixed effects
df.id3 = categorical(Int.(rand(1:15, N)))
@time reg(df, @model(y ~ x1 + id3))
using GLM
@time lm(@formula(y ~ x1 + id3), df)

# compare with R data.table
using DataFrames
function f(y, id)
    cache = zeros(eltype(y), length(id.pool))
    for i in 1:length(y)
        cache[id.refs[i]] += y[i]
    end
    out = similar(y)
    for i in 1:length(y)
        out[i] = cache[id.refs[i]]
    end
    return out
end
N = 100000000
K = 100
y = randn(N)
id1 = categorical(Int.(rand(1:(N/K), N)))
id2 = categorical(Int.(rand(1:K, N)))
@time f(y, id1)
#1.9
@time f(y, id2)
#0.5

# in R
N = 100000000
K = 100
y = rnorm(N)
id1 = sample(floor(N/K), N, TRUE)
id2 = sample(K, N, TRUE)
dt = data.table(y = y, id1 = id1, id2 = id2)
setkey(dt, id1)
system.time(dt[, list(y1 = mean(y)), by = id1])
#1.3
setkey(dt, id2)
system.time(dt[, list(y1 = mean(y)), by = id2])
#1.3


text:
include("utilities.jl")

using SEAL
using Printf

"""
    example_rotation_bfv()

Perform some rotation on data encryped with the BFV scheme. This routine is based on the file
`native/examples/5_rotation.cpp` of the original SEAL library and should yield the exact same
output.

* [SEAL](https://github.com/microsoft/SEAL)
* [native/examples/5_rotation.cpp](https://github.com/microsoft/SEAL/blob/master/native/examples/5_rotation.cpp)
"""
function example_rotation_bfv()
  print_example_banner("Example: Rotation / Rotation in BFV")

  parms = EncryptionParameters(SchemeType.bfv)

  poly_modulus_degree = 8192
  set_poly_modulus_degree!(parms, poly_modulus_degree)

  set_coeff_modulus!(parms, coeff_modulus_bfv_default(poly_modulus_degree))
  set_plain_modulus!(parms, plain_modulus_batching(poly_modulus_degree, 20))

  context = SEALContext(parms)
  print_parameters(context)
  println()

  keygen = KeyGenerator(context)
  secret_key_ = secret_key(keygen)
  public_key_ = PublicKey()
  create_public_key!(public_key_, keygen)
  relin_keys_ = RelinKeys()
  create_relin_keys!(relin_keys_, keygen)
  encryptor = Encryptor(context, public_key_)
  evaluator = Evaluator(context)
  decryptor = Decryptor(context, secret_key_)

  batch_encoder = BatchEncoder(context)
  slot_count_ = slot_count(batch_encoder)
  row_size = div(slot_count_, 2)
  println("Plaintext matrix row size: ", row_size)

  pod_matrix = zeros(UInt64, slot_count_)
  pod_matrix[1] = 0
  pod_matrix[2] = 1
  pod_matrix[3] = 2
  pod_matrix[4] = 3
  pod_matrix[row_size + 1] = 4
  pod_matrix[row_size + 2] = 5
  pod_matrix[row_size + 3] = 6
  pod_matrix[row_size + 4] = 7

  println("Input plaintext matrix:")
  print_matrix(pod_matrix, row_size)

  plain_matrix = Plaintext()
  print_line(@__LINE__)
  println("Encode and encrypt.")
  encode!(plain_matrix, pod_matrix, batch_encoder)
  encrypted_matrix = Ciphertext()
  encrypt!(encrypted_matrix, plain_matrix, encryptor)
  println("    + Noise budget in fresh encryption: ",
          invariant_noise_budget(encrypted_matrix, decryptor),
          " bits")
  println()

  galois_keys = GaloisKeys()
  create_galois_keys!(galois_keys, keygen)

  print_line(@__LINE__)
  println("Rotate rows 3 steps left.")
  rotate_rows_inplace!(encrypted_matrix, 3, galois_keys, evaluator)
  plain_result = Plaintext()
  println("    + Noise budget after rotation: ",
          invariant_noise_budget(encrypted_matrix, decryptor),
          " bits")
  println("    + Decrypt and decode ...... Correct.")
  decrypt!(plain_result, encrypted_matrix, decryptor)
  decode!(pod_matrix, plain_result, batch_encoder)
  print_matrix(pod_matrix, row_size)

  print_line(@__LINE__)
  println("Rotate columns.")
  rotate_columns_inplace!(encrypted_matrix, galois_keys, evaluator)
  println("    + Noise budget after rotation: ",
          invariant_noise_budget(encrypted_matrix, decryptor),
          " bits")
  println("    + Decrypt and decode ...... Correct.")
  decrypt!(plain_result, encrypted_matrix, decryptor)
  decode!(pod_matrix, plain_result, batch_encoder)
  print_matrix(pod_matrix, row_size)

  print_line(@__LINE__)
  println("Rotate rows 4 steps right.")
  rotate_rows_inplace!(encrypted_matrix, -4, galois_keys, evaluator)
  plain_result = Plaintext()
  println("    + Noise budget after rotation: ",
          invariant_noise_budget(encrypted_matrix, decryptor),
          " bits")
  println("    + Decrypt and decode ...... Correct.")
  decrypt!(plain_result, encrypted_matrix, decryptor)
  decode!(pod_matrix, plain_result, batch_encoder)
  print_matrix(pod_matrix, row_size)
end

"""
    example_rotation_ckks()

Perform some rotation on data encryped with the CKKS scheme. This routine is based on the file
`native/examples/5_rotation.cpp` of the original SEAL library and should yield the exact same
output.

* [SEAL](https://github.com/microsoft/SEAL)
* [native/examples/5_rotation.cpp](https://github.com/microsoft/SEAL/blob/master/native/examples/5_rotation.cpp)
"""
function example_rotation_ckks()
  print_example_banner("Example: Rotation / Rotation in CKKS")

  parms = EncryptionParameters(SchemeType.ckks)

  poly_modulus_degree = 8192
  set_poly_modulus_degree!(parms, poly_modulus_degree)
  set_coeff_modulus!(parms, coeff_modulus_create(poly_modulus_degree, [40, 40, 40, 40, 40]))

  context = SEALContext(parms)
  print_parameters(context)
  println()

  keygen = KeyGenerator(context)
  public_key_ = PublicKey()
  create_public_key!(public_key_, keygen)
  secret_key_ = secret_key(keygen)
  relin_keys_ = RelinKeys()
  create_relin_keys!(relin_keys_, keygen)
  galois_keys_ = GaloisKeys()
  create_galois_keys!(galois_keys_, keygen)
  encryptor = Encryptor(context, public_key_)
  evaluator = Evaluator(context)
  decryptor = Decryptor(context, secret_key_)

  encoder = CKKSEncoder(context)

  slot_count_ = slot_count(encoder)
  println("Number of slots: ", slot_count_)
  input = collect(range(0.0, 1.0, length=slot_count_))
  println("Input vector:")
  print_vector(input, 3, 7)

  initial_scale = 2.0^50

  print_line(@__LINE__)
  println("Encode and encrypt.")
  plain = Plaintext()
  encode!(plain, input, initial_scale, encoder)
  encrypted = Ciphertext()
  encrypt!(encrypted, plain, encryptor)

  rotated = Ciphertext()
  print_line(@__LINE__)
  println("Rotate 2 steps left.")
  rotate_vector!(rotated, encrypted, 2, galois_keys_, evaluator)
  println("    + Decrypt and decode ...... Correct.")
  decrypt!(plain, rotated, decryptor)
  result = similar(input)
  decode!(result, plain, encoder)
  print_vector(result, 3, 7)

  return
end

"""
    example_rotation()

Perform some rotation on data encryped with the BFV and CKKS schemes. This routine is based on the
file `native/examples/5_rotation.cpp` of the original SEAL library and should yield the exact same
output.

* [SEAL](https://github.com/microsoft/SEAL)
* [native/examples/5_rotation.cpp](https://github.com/microsoft/SEAL/blob/master/native/examples/5_rotation.cpp)
"""
function example_rotation()
  print_example_banner("Example: Rotation")

  example_rotation_bfv()
  example_rotation_ckks()
end


text:
module AnimatedPlots

importall SFML

include("types.jl")
include("graph.jl")
include("plotwindow.jl")
include("axis.jl")
include("gifmaker.jl")

windows = PlotWindow[]

function plot(graph::Graph, name, width, height)
	window = create_window(name=name, width=width, height=height)
	add_graph(window, graph)
	redraw(window)
	t = open_window(window)
	window.task = t
	nothing
end
function plot(fun::Function, name, width, height)
	plot(StaticGraph(fun), name, width, height)
end

function plot(graph::Graph; new_window=false)
	if length(windows) > 0 && !new_window
		add_graph(windows[end], graph)
		redraw(windows[end])
		nothing
	else
		plot(graph, "Plot", 800, 600)
	end
end
function plot(fun::Function; new_window=false)
	if length(windows) > 0 && !new_window
		plot(StaticGraph(fun))
	else
		plot(fun, "Plot", 800, 600)
	end
end

function remove(graph::Graph)
	splice!(windows[end].graphs, find(windows[end].graphs .== graph)[1])
	nothing
end

function create_window(window::RenderWindow)
	plotwindow = PlotWindow(window)
	redraw(plotwindow)
	push!(windows, plotwindow)
	return plotwindow
end

function create_window(; name="Plot", width=800, height=600)
	window = PlotWindow(name, width, height)
	redraw(window)
	push!(windows, window)
	return window
end

function open_window(window::PlotWindow)
	@async begin
		while isopen(window)
			sleep(0)
			check_input(window)
			redraw(window)
			clear(window.renderwindow, SFML.white)
			draw(window)
			display(window.renderwindow)
		end
		splice!(windows, find(windows .== window)[1])
	end
end

function current_window()
	return windows[end]
end

function follow(graph::AnimatedGraph)
	follow(current_window(), graph)
end

function unfollow()
	unfollow(current_window())
end

function screenshot(filename::AbstractString)
	image = capture(current_window().renderwindow)
	save_to_file(image, filename)
	println("Saved plot to '$filename'")
end

function make_gif(width, height, duration, filename="plot.gif", delay="0.06")
	make_gif(current_window(), width, height, duration, filename, delay)
end

function set_zoom(zoom::Real)
	current_window().ppu = zoom
	redraw(current_window(), true)
end

export open_window, create_window, plot, remove, current_window, get_window,
follow, unfollow, screenshot, make_gif, set_zoom

end


text:
using LinearAlgebra
using Plots
using Revise

includet("src/IterativeMethods.jl")
includet("src/InternalArray.jl")

x0 = [0, 300]
v0 = [100, 0]
m = 1
c = 0.1
g = [0, -9.8]

# Solving for speed
v = eulermethod(0, v0) do t, v
	g
end
p1 = plot(InternalArray(v, 1), InternalArray(v, 2), label = "v (no)")

# Solving por position
x = eulermethod(0, x0) do t, x
	v0 + t*g
end
p2 = plot(InternalArray(x, 1), InternalArray(x, 2), label = "x (no)")

# Solving for speed with friction
v = eulermethod(0, v0) do t, v
	g-c/m*norm(v)*v
end
p3 = plot(InternalArray(v, 1), InternalArray(v, 2), label = "v (yes)")

# Solving for position with friction
x = eulermethod(0, [x0; v0]) do t, x
	v = [x[3], x[4]]
	return [v; g-c/m*norm(v)*v]
end
p4 = plot(InternalArray(x, 1), InternalArray(x, 2), label = "x (yes)")

plot(p1, p2, p3, p4, plot_title = "Position and speed with and without friction")

max_stars_repo_name: useredsa/numeric-differential-equations

max_issues_repo_name: useredsa/numeric-differential-equations

max_forks_repo_name: useredsa/numeric-differential-equations


text:
mutable struct CircularDeque{T}
    buffer::Vector{T}
    capacity::Int
    n::Int
    first::Int
    last::Int
end

"""
    CircularDeque{T}(n)

Create a double-ended queue of maximum capacity `n`, implemented as a circular buffer. The element type is `T`.
"""
CircularDeque{T}(n::Int) where {T} = CircularDeque(Vector{T}(undef, n), n, 0, 1, n)

Base.length(D::CircularDeque) = D.n
Base.eltype(::Type{CircularDeque{T}}) where {T} = T

"""
    capacity(D::CircularDeque)

Return the capacity of the circular deque
"""
capacity(D::CircularDeque) = D.capacity

function Base.empty!(D::CircularDeque)
    D.n = 0
    D.first = 1
    D.last = D.capacity
    return D
end

Base.isempty(D::CircularDeque) = D.n == 0

"""
    first(D::CircularDeque)

Get the item at the front of the queue.
"""
@inline function first(D::CircularDeque)
    @boundscheck D.n > 0 || throw(BoundsError())
    return @inbounds D.buffer[D.first]
end

"""
    last(D::CircularDeque)

Get the item from the back of the queue.
"""
@inline function last(D::CircularDeque)
    @boundscheck D.n > 0 || throw(BoundsError())
    return @inbounds D.buffer[D.last]
end

@inline function Base.push!(D::CircularDeque, v)
    @boundscheck D.n < D.capacity || throw(BoundsError()) # prevent overflow
    D.n += 1
    tmp = D.last + 1
    D.last = ifelse(tmp > D.capacity, 1, tmp)  # wraparound
    @inbounds D.buffer[D.last] = v
    return D
end

@inline Base.@propagate_inbounds function Base.pop!(D::CircularDeque)
    v = last(D)
    D.n -= 1
    tmp = D.last - 1
    D.last = ifelse(tmp < 1, D.capacity, tmp)
    v
end

"""
    pushfirst!(D::CircularDeque, v)

Add an element to the front.
"""
@inline function pushfirst!(D::CircularDeque, v)
    @boundscheck D.n < D.capacity || throw(BoundsError())
    D.n += 1
    tmp = D.first - 1
    D.first = ifelse(tmp < 1, D.capacity, tmp)
    @inbounds D.buffer[D.first] = v
    D
end

"""
    popfirst!(D::CircularDeque)

Remove the element at the front.
"""
@inline Base.@propagate_inbounds function popfirst!(D::CircularDeque)
    v = first(D)
    D.n -= 1
    tmp = D.first + 1
    D.first = ifelse(tmp > D.capacity, 1, tmp)
    v
end

# getindex sans bounds checking
@inline function _unsafe_getindex(D::CircularDeque, i::Integer)
    j = D.first + i - 1
    if j > D.capacity
        j -= D.capacity
    end
    @inbounds ret = D.buffer[j]
    return ret
end

@inline function Base.getindex(D::CircularDeque, i::Integer)
    @boundscheck 1 <= i <= D.n || throw(BoundsError())
    return _unsafe_getindex(D, i)
end

# Iteration via getindex
@inline function iterate(d::CircularDeque, i = 1)
    i == d.n + 1 ? nothing : (_unsafe_getindex(d, i), i+1)
end

function Base.show(io::IO, D::CircularDeque{T}) where T
    print(io, "CircularDeque{$T}([")
    for i = 1:length(D)
        print(io, D[i])
        i < length(D) && print(io, ',')
    end
    print(io, "])")
end


text:
using CTechCommon
# using FactCheck
using Base.Test

# @testset "all" begin

# @testset "arrays" begin
#   r = 1:10
#   x = collect(r)
#   @test get(getith(r,3)) == 3
#   @test get(getith(x,3)) == 3
#   @test isnull(getith(x,11)) == true

#   z = ([1,"one"], [2,"two"])
#   z1, z2 = unzip(z)
#   @test z1 == [1,2]
#   @test z2 == ["one","two"]
#   @test sizes(z) == ((2,),(2,))
#   @test mapf((sin,cos), 1) == [sin(1), cos(1)]

#   m = rand(3,3)
#   v = rand(3)
#   @test row(m, 1)[2] == m[1,2]
#   @test col(m, 1)[2] == m[2,1]
#   vw = rows(m, 1:2)
#   vw[1,1] = 0
#   @test m[1,1] == 0
#   @test nrows(m) == 3
#   @test ncols(m) == 3
#   @test nrows(v) == 3
#   @test ncols(v) == 1

#   m = reshape(collect(1:4),2,2)
#   @test addOnes(collect(1:3)) == Int[1,2,3,1]
#   @test addOnes(collect(1:3.)) == Float64[1,2,3,1]
#   @test addOnes(m) == Int[1 3 1; 2 4 1]

#   @test getPctOfInt(0.5, 10) == 5
#   @test getPctOfInt(1.5, 10) == 10
#   @test splitRange(10, 0.5) == (1:5, 6:10)
#   @test splitMatrixRows(m, 0.5)[1] == rows(m, 1:1)

#   @test stringfloat(0.1111) == "0.111"
#   @test stringfloat(0.1111, 1) == "0.1"
#   @test stringfloats([0.1111, 0.9999], 2) == "[0.11, 1.00]"
# end

# facts("misc") do
#   @test donothing(5,5) == nothing
#   @test nop(5) == 5
#   @test returntrue(5) == true
# end

@testset "price" begin
  @test Price(50.55) == CTechCommon.makePrice(5055, 2)
  @test Price(50.55, 4) == Price(50.55, 2)

  fp = 55.55555
  p = Price(fp)
  @test p.priceLong == 555556
  @test p.precision == 4
  @test p.multiplier == 10_000
  @test CTechCommon.getLong(p, 2) == 5556
  @test float(p + Price(1.01,2)) ≈ 56.5656

  p = Price(fp, 2)
  rounded = round(fp, 2)
  @test p + 0.001 == p
  @test float(p + 0.005) ≈ rounded + 0.01
  @test float(p * 2) ≈ rounded * 2
end

@testset "time" begin
  # function jacob_calcSecondsEpochToMidnight(secondsSinceEpoch::Integer)
  #   utc = DateTime(Date(Dates.unix2datetime(secondsSinceEpoch)))
  #   adjustment = Dates.Second(div(Dates.value(now() - now(Dates.UTC)),1000))
  #   return Dates.value(utc+adjustment)
  # end

  year, month, day = 2013,1,1
  @test CTechCommon.getHoursAdjustmentFromUTC(year,month,day) == 5

  @test TimeOfDay("10") == 10 * nanosInOneHour
  @test TimeOfDay("10:30") == 10 * nanosInOneHour + 30 * nanosInOneMinute
  @test TimeOfDay("10:30:05") == 10 * nanosInOneHour + 30 * nanosInOneMinute + 5 * nanosInOneSecond
  @test TimeOfDay("10:30:05.00001") == 10 * nanosInOneHour + 30 * nanosInOneMinute + 5 * nanosInOneSecond + 10 * nanosInOneMicro

  t = TimeOfDay("10:30:05")
  buf = IOBuffer()
  print(buf, t)
  @test String(buf) == "10:30:05.000000"

  tdiff = TimeOfDay(nanosInOneMinute)
  trange = t : tdiff : t + 5 * nanosInOneMinute
  @test length(trange) == 6
  @test trange[2] == TimeOfDay("10:31:05")

  @test iszero(TimeOfDay(0)) == true
  @test iszero(TimeOfDay(1)) == false
  @test_throws CTechCommon.NegativeTimeOfDayError{Int64} TimeOfDay(-1)

end

@testset "fixedsym" begin
  str = "abc"

  n = 6
  buf = IOBuffer(n)
  for c in str
    write(buf, UInt8(c))
  end
  for i in buf.size+1:n
    write(buf, UInt8(0))
  end

  buf.ptr = 1
  sym = Symbol6(buf)
  @test sizeof(sym) == n
  @test String(sym) == str

  n = 8
  buf = IOBuffer(n)
  for c in str
    write(buf, UInt8(c))
  end
  for i in buf.size+1:n
    write(buf, UInt8(0))
  end

  buf.ptr = 1
  sym = Symbol8(buf)
  @test sizeof(sym) == n
  @test String(sym) == str
end

@testset "bufferedio" begin
end

@testset "trie" begin
end

@testset "Logger" begin
  @test log_severity() == InfoSeverity
  f = open(tempname(),"w")
  @test log_io!(f) == nothing
  @test log_io() == f
  close(f)
  log_severity!(DebugSeverity)
  @test log_severity() == DebugSeverity
  log_severity!(ErrorSeverity)
  @test log_severity() == ErrorSeverity
  log_severity!(InfoSeverity)
end

@testset "markets" begin
  @test typeof(getFee(EDGX, false)) == Float64
  # @test EDGA < EDGX
  @test EXCH_SORTED_BY_TAKE_FEE == [EDGA,EDGX]

  id = generateOID()
  @test generateOID() == id + 1
  @test generateOID() == id + 2

  s = "MSFT"
  @test String(Ticker(s)) == s
  @test Ticker(s) == s
  @test Ticker(s) < "NSFT"
end

type TmpVal; val::Int; end
updatetmp(tmpval::TmpVal, newval::Integer) = (tmpval.val = newval; nothing)

type Counter
  n::Int
  Counter() = new(0)
end

@testset "pubsub" begin
  # TODO test:
  #   various filters
  #   ordering of connections
  #   function types, argument lists
  #   unregistering

  # simple test function to add the sum of the args to the counter
  handlemsg(c::Counter, args...) = (c.n += sum(args))

  reset_hub()
  pub1 = Publisher(handlemsg, Filters([:key1, 5, 6], [:key2, "hi", "yo"]))

  c1 = Counter()
  sub1 = subscribe(handlemsg, c1)

  # sub1 should add 1 to c1
  publish(pub1, 1)
  @test c1.n == 1

  # now both sub1 and sub2 should be subscribed to pub1, so this should add 9 (2+3+4) twice
  sub2 = subscribe(handlemsg, c1, Filters([:key1, 6, 7, 8]))
  publish(pub1, 2,3,4)
  @test c1.n == 19

  # unregistering sub1, so the next pub should add 10 once
  unregister(sub1)
  publish(pub1, 5,5)
  @test c1.n == 29

  # sub3 should not be subscribed to pub1
  sub3 = subscribe(handlemsg, c1, Filters([:key2, "word"]))
  publish(pub1, 10)
  @test c1.n == 39

end

const tdiff = TimeOfDay("00:00:30")
const tstart = TimeOfDay("9:35:30")
const tstop = TimeOfDay("9:39:15")

@testset "scheduler" begin

  function periodiccallback(c, x...)
    c.n += sum(x)
    @LOG c x
    if c.n == 11
      stopScheduler(TimeOfDay("9:39:15"))
      LOG("stopScheduler")
    end
    schedule(NOW() + tdiff, pub, 2)
  end

  # set up a simple pub/sub, where the listener object is a Counter.
  # every time the event triggers, it reschedules tdiff (30 seconds) ahead.
  # the first pub at tstart has 1 as the arg, the rest have 2.
  # counter should reach 11 on the 6th pub at 9:38, at which time it schedules
  # the stop for tstop.
  # 2 more callbacks are scheduled and run before we hit the stop event at 9:39:15
  reset_hub()
  c = Counter()
  pub = Publisher(periodiccallback)
  sub = subscribe(periodiccallback, c)

  initScheduler()               # resets the scheduler (unnecessary for first pass)
  schedule(tstart, pub, 1)      # schedules first pub
  stopped = processAllEvents()  # main loop

  @test stopped == true
  @test NOW() == tstop
  @test c.n == 15
end

# end # all

# FactCheck.exitstatus()


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "GnuPG"
version = v"2.2.27"

# Collection of sources required to build libgcrypt
sources = [
    ArchiveSource("https://gnupg.org/ftp/gcrypt/gnupg/gnupg-$(version).tar.bz2",
                  "34e60009014ea16402069136e0a5f63d9b65f90096244975db5cea74b3d02399"),
    DirectorySource("./bundled"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/gnupg-*/
if [[ "${target}" != ${MACHTYPE} ]]; then
    # Delete `gpg-error-config` of the host to prevent it from being picked up
    # when configuring the package
    rm "${host_bindir}/gpg-error-config"
    # Use Windows LDAP
    FLAGS=(LDAPLIBS="-lwldap32")
fi
if [[ "${target}" == *86*-linux-gnu ]]; then
    # We have an old glibc which doesn't have `IN_EXCL_UNLINK`
    FLAGS=(ac_cv_func_inotify_init=no)
    # Add -lrt dependency to fix the error
    #     undefined reference to `clock_gettime'
    atomic_patch -p1 ../patches/intel-linux-gnu-add-rt-lib.patch
fi
./configure --prefix=${prefix} \
    --build=${MACHTYPE} \
    --host=${target} \
    "${FLAGS[@]}"
make -j${nproc}
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line.  We are manually disabling
# many platforms that do not seem to work.
platforms = supported_platforms()

# The products that we will ensure are always built
products = [
    ExecutableProduct("dirmngr", :dirmngr),
    ExecutableProduct("dirmngr-client", :dirmngr_client),
    ExecutableProduct("gpg", :gpg),
    ExecutableProduct("gpg-agent", :gpg_agent),
    ExecutableProduct("gpgconf", :gpgconf),
    ExecutableProduct("gpg-connect-agent", :gpg_connect_agent),
    ExecutableProduct("gpgscm", :gpgscm),
    ExecutableProduct("gpgsm", :gpgsm),
    # ExecutableProduct("gpgsplit", :gpgsplit),
    ExecutableProduct("gpgtar", :gpgtar),
    ExecutableProduct("gpgv", :gpgv),
    # ExecutableProduct("gpg-wks-server", :gpg_wks_server),
    ExecutableProduct("kbxutil", :kbxutil),
    # ExecutableProduct("watchgnupg", :watchgnupg),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    # We need this to run a host yat2m executable
    HostBuildDependency("Libgpg_error_jll"),
    # We need this to run a host msgfmt executable
    HostBuildDependency("Gettext_jll"),
    Dependency("GnuTLS_jll"),
    Dependency("Libksba_jll"),
    Dependency("Libgcrypt_jll"),
    Dependency("Libgpg_error_jll"),
    Dependency("nPth_jll"),
    Dependency("Zlib_jll"),
    Dependency("Libassuan_jll"),
    Dependency("OpenLDAPClient_jll"),
    Dependency("Bzip2_jll"),
    Dependency("SQLite_jll"),
    Dependency("libusb_jll"),
    Dependency("Nettle_jll", v"3.4.1", compat="~3.4.1"),
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: bzinberg/Yggdrasil

max_issues_repo_name: bzinberg/Yggdrasil

max_forks_repo_name: bzinberg/Yggdrasil


text:
##################################################
## Test of the translation of a solid expansion ##
##################################################
@testset "Translation of Rlm" begin

    ε = 1000*eps(Float64)
    @polyvar x y z;

    Wert = zeros(25,2)
    v = randn(3) # Verschiebung
    pkt = randn(3) # Testpunkt

    vp = pkt+v;

    for l=0:4
        for m=-l:l
            q = SphericalHarmonicExpansions.zlm(l,m,x,y,z);
            p = SphericalHarmonicExpansions.translateRlm(l,m,v[1],v[2],v[3]);

            polyRlm = @fastfunc q;
            polyTranslateRlm = @fastfunc p;

            Wert[l*(l+1)+m+1,1] = polyRlm(vp[1],vp[2],vp[3]);
            Wert[l*(l+1)+m+1,2] = polyTranslateRlm(pkt[1],pkt[2],pkt[3]);
        end
    end

    @test isapprox(Wert[:,1],Wert[:,2],atol = ε)
end

#################################################
## Test of the translation of the coefficients ##
#################################################
@testset "Translation of Clm" begin

    @polyvar x y z

    ###################################
    # Translation of the coefficients #
    # with v and after that with -v   #
    ###################################
    # Shift vector
    v = randn(3)

    # Coefficients
    L = 6
    C = SphericalHarmonicCoefficients(randn((L+1)^2))

    # Translation
    coeffsTrans = translation(C,v)
    coeffsBack = translation(coeffsTrans,-v)

    # Test
    @test isapprox(C,coeffsBack)

    ############################################
    # Konkretes Beispiel für die Funktion xyz  #
    # solid & spherical Expansion              #
    ############################################
    # Coefficients
    Cspher = SphericalHarmonicCoefficients(zeros(25));
    Cspher[3,-2] = 1/sqrt(15)*sqrt((4*pi)/(2*3+1));
    Csolid = SphericalHarmonicCoefficients(zeros(25),1.0,true);
    Csolid[3,-2] = 1/sqrt(15);

    # Translation
    v = [1.,-1.,2.]
    coeffsTransSpher = translation(Cspher,v)
    coeffsTransSolid = translation(Csolid,v)

    # Correct solution
    CtransSpher = zeros(25)
    CtransSpher[1:11] = sqrt(4pi/3).*[-2*sqrt(3), 2, -1, -2, 2/sqrt(5), 1/sqrt(5), 0, -1/sqrt(5), 0, 0, 1/sqrt(35)]
    CtransSpher = SphericalHarmonicCoefficients(CtransSpher)
    CtransSolid = zeros(25)
    CtransSolid[1:11] = [-2, 2, -1, -2, 2/sqrt(3), 1/sqrt(3), 0, -1/sqrt(3), 0, 0, 1/sqrt(15)]
    CtransSolid = SphericalHarmonicCoefficients(CtransSolid,1.0,true)

    # Test
    @test isapprox(CtransSpher,coeffsTransSpher)
    @test isapprox(CtransSolid,coeffsTransSolid)
end