text:
include("model.jl")
using Flux, Flux.Tracker
import Random

Random.seed!(0);

# Generate Data
N = 1000
x = randn(N)
b0 = 2.0
b1 = -3.0
sig = 0.5
y = b0 .+ b1 .* x .+ randn(N) * sig

minibatch_size = 100
vp = VP()
loss(y, x) = -elbo(y, x, vp, N) / N

params = Tracker.Params([getfield(vp, fn) for fn in fieldnames(typeof(vp))])
# grads = Tracker.gradient(() -> loss(y, x), params)

opt = ADAM(1e-1)
niters = 5000

@time for i in 1:niters
  idx = sample(1:N, minibatch_size)
  Flux.train!(loss, params, [(y[idx], x[idx])], opt)
  if i % (niters / 10) == 0
    println("Progress: $(i)/$(niters) | ELBO: $(-loss(y, x).data)")
  end
end

sig_post = exp.([rsample(vp.log_sig) for b in 1:1000])
println("b0_mean: $(vp.b0[1].data) | b1_mean: $(vp.b1[1].data) | sig_mean: $(mean(sig_post).data)")
println("Loss: $(loss(y, x).data)")


text:
# This file uses constants defined in constant.jl
#
# Follow Bryan and Cox (1972) : An Approximate Equation of State for Numerical Models of Ocean Circulation
# Here we use Table 3, and arbitrary pick Z=0m as the formula.
#

const T_ref =   13.5
const S_ref =   32.6
const ρ_ref = 1024.458

const _ρ1 = -.20134    / ρ_ref
const _ρ2 =  .77096    / ρ_ref
const _ρ3 = -.49261e-2 / ρ_ref
const _ρ4 =  .46092e-3 / ρ_ref
const _ρ5 = -.20105e-2 / ρ_ref
const _ρ6 =  .36597e-4 / ρ_ref
const _ρ7 =  .47371e-5 / ρ_ref
const _ρ8 =  .37735e-4 / ρ_ref
const _ρ9 =  .65493e-5 / ρ_ref

@inline function TS2b(T::Float64, S::Float64)
    ΔT = T - T_ref
    ΔS = S - S_ref
    return - g * ( _ρ1*ΔT + _ρ2*ΔS +  _ρ3*(ΔT^2) + _ρ4*(ΔS^2) + _ρ5*ΔT*ΔS + _ρ6*(ΔT^3) + _ρ7*(ΔS^2)*ΔT + _ρ8*(ΔT^2)*ΔS + _ρ9*(ΔS^3) )
end

@inline function TS2α(T, S)
    ΔT = T - T_ref
    ΔS = S - S_ref
    return  - ( _ρ1 + 2*_ρ3*ΔT + _ρ5*ΔS + 3*_ρ6*(ΔT^2) + _ρ7*(ΔS^2) + 2*_ρ8*ΔS*ΔT )
end

@inline function TS2β(T, S)
    ΔT = T - T_ref
    ΔS = S - S_ref
    return  _ρ2 + 2*_ρ4*ΔS + _ρ5*ΔT + 2*_ρ7*ΔS*ΔT + _ρ8*(ΔT^2) + 3*_ρ9*(ΔS^2)
end

ext: jl

max_stars_repo_name: meteorologytoday/IOM

max_issues_repo_name: meteorologytoday/IOM

max_forks_repo_name: meteorologytoday/IOM


text:
### Fetch required packages ###
using DiffEqBase, Catalyst, DiffEqJump, Random, Statistics, Test
using ModelingToolkit: get_states, get_ps
using StableRNGs
rng = StableRNG(12345)

### Declares a test network. ###
higher_order_network_1 = @reaction_network begin
    p,                                                          ∅ ⟼ X1
    r1,                                                         2X1 ⟼ 3X2
    mm(X1,r2,K),                                                3X2 ⟼ X3 + 2X4
    r3,                                                         X3 + 2X4 ⟼ 3X5 + 3X6
    r4*X2,                                                      3X5 + 3X6 ⟼ 3X5 + 2X7 + 4X8
    r5,                                                         3X5 + 2X7 + 4X8 ⟼ 10X9
    r6,                                                         10X9 ⟼ X10
    d,                                                          2X10 ⟼ ∅
end p r1 r2 K r3 r4 r5 r6 d

### Tests that deterministic and stochastic differential functions are identical. ###

higher_order_network_2 = @reaction_network begin
    p,                                                          ∅ ⟾ X1
    r1*X1^2/factorial(2),                                       2X1 ⟾ 3X2
    mm(X1,r2,K)*X2^3/factorial(3),                              3X2 ⟾ X3 + 2X4
    r3*X3*X4^2/factorial(2),                                    X3 + 2X4 ⟾ 3X5 + 3X6
    r4*X2*X5^3*X6^3/(factorial(3)*factorial(3)),                3X5 + 3X6 ⟾ 3X5 + 2X7 + 4X8
    r5*X5^3*X7^2*X8^4/(factorial(3)*factorial(2)*factorial(4)), 3X5 + 2X7 + 4X8 ⟾ 10X9
    r6*X9^10/factorial(10),                                     10X9 ⟾ X10
    d*X10^2/factorial(2),                                       2X10 ⟾ ∅
end p r1 r2 K r3 r4 r5 r6 d

f1 = ODEFunction(convert(ODESystem,higher_order_network_1),jac=true)
f2 = ODEFunction(convert(ODESystem,higher_order_network_2),jac=true)
g1 = SDEFunction(convert(SDESystem,higher_order_network_1))
g2 = SDEFunction(convert(SDESystem,higher_order_network_2))
for factor in [1e-2, 1e-1, 1e0, 1e1, 1e2, 1e3]
    u0 = factor*rand(rng,length(get_states(higher_order_network_1)))
    p = factor*rand(rng,length(get_ps(higher_order_network_2)))
    t = rand(rng)
    @test all(abs.(f1(u0,p,t) .- f2(u0,p,t)) .< 100*eps())
    @test all(abs.(f1.jac(u0,p,t) .- f2.jac(u0,p,t)) .< 100*eps())
    @test all(abs.(g1(u0,p,t) .- g2(u0,p,t)) .< 100*eps())
end

### Tests that the discrete jump systems are equal. ###
higher_order_network_3 = @reaction_network begin
    p,                                                  ∅ ⟼ X1
    r1*binomial(X1,2),                                  2X1 ⟾ 3X2
    mm(X1,r2,K)*binomial(X2,3),                         3X2 ⟾ X3 + 2X4
    r3*binomial(X3,1)*binomial(X4,2),                   X3 + 2X4 ⟾ 3X5 + 3X6
    r4*X2*binomial(X5,3)*binomial(X6,3),                3X5 + 3X6 ⟾ 3X5 + 2X7 + 4X8
    r5*binomial(X5,3)*binomial(X7,2)*binomial(X8,4),    3X5 + 2X7 + 4X8 ⟾ 10X9
    r6*binomial(X9,10),                                 10X9 ⟾ X10
    d*binomial(X10,2),                                  2X10 ⟾ ∅
end p r1 r2 K r3 r4 r5 r6 d

for factor in [1e-1, 1e0], repeat = 1:5
    u0 = rand(rng,1:Int64(factor*100),length(get_states(higher_order_network_1)))
    p = factor*rand(rng,length(get_ps(higher_order_network_3)))
    prob1 = JumpProblem(higher_order_network_1,DiscreteProblem(higher_order_network_1,u0,(0.,1000.),p),Direct())
    sol1 = solve(prob1,SSAStepper())
    prob2 = JumpProblem(higher_order_network_3,DiscreteProblem(higher_order_network_3,u0,(0.,1000.),p),Direct())
    sol2 = solve(prob2,SSAStepper())
    for i = 1:length(u0)
        vals1 = getindex.(sol1.u,i);
        vals2 = getindex.(sol1.u,i);
        (mean(vals2)>0.001) && @test 0.8 < mean(vals1)/mean(vals2) < 1.25
        (std(vals2)>0.001) && @test 0.8 < std(vals1)/std(vals2) < 1.25
    end
end


text:
module SnllsSEM

using StructuralEquationModels

include("imply/SNLLS.jl")
include("loss/SNLLS.jl")

export SNLLS, SemSNLLS

end # module


text:
#   ribbon modulation example
#

using ShapeFactory

app = dispatch("CATIA.Application")
sel = app.ActiveDocument.Selection
prt = app.ActiveDocument.Part

a0 = prt.Parameters.CreateDimension("from_angle", "ANGLE", 0)
a0.ValuateFromString("45deg")

a1 = prt.Parameters.CreateDimension("to_angle", "ANGLE", 0)
a1.ValuateFromString("1deg")

rw = prt.Parameters.CreateDimension("ribbon_width", "LENGTH", 0)
rw.ValuateFromString("300mm")

rh = prt.Parameters.CreateDimension("ribbon_height", "LENGTH", 0)
rh.ValuateFromString("1000mm")

ri = prt.Parameters.CreateDimension("ribbon_inset", "LENGTH", 0)
ri.ValuateFromString("50mm")

fi = prt.Parameters.CreateDimension("flange_inset", "LENGTH", 0)
fi.ValuateFromString("20mm")

rn = 10

gst = prt.HybridBodies.Add()
gst.Name = "Output"

fac = prt.HybridShapeFactory
xyp = prt.FindObjectByName("xy plane")

ydr = fac.AddNewDirectionByCoord(0, 1, 0)
zdr = fac.AddNewDirectionByCoord(0, 0, 1)

function Base.repr(idsp::Ptr{IDispatch})
    return join(split(idsp.Name, "\\")[2:end], "\\")
end

function relate(idsp::Ptr{IDispatch}, frm::String)
    name = replace(repr(idsp), "\\" => "_")
    prt.Relations.CreateFormula(name, "", idsp, frm)
end

function blend_surface(edges, flip)
    n = (1, 2)
    e = fac.AddNewExtrude.(edges, 0, (2 * -flip, 2 * flip), zdr)
    b = fac.AddNewBlend()
    b.Coupling = 4
    b.SetCurve.(n, edges)
    b.SetSupport.(n, e)
    b.SetContinuity.(n, 1)
    b.SetTransition.(n, (flip, -flip))
    b.RuledDevelopableSurface = true
    return b
end

function flange_surface(edges, plane, limit1, limit2, flip)
    pin = fac.AddNewPointOnCurveFromDistance(edges[2], 0, false)
    ppn = fac.AddNewLineAngle(edges[2], xyp, pin, false, 0, 1, 90, false)
    relate(pin.Ratio, "$(repr(fi))")
    ext = fac.AddNewExtrude(pin, 0, 1, zdr)
    relate(ext.EndOffset, "$(repr(rh)) * -$(flip)")
    fln = fac.AddNewExtrude(ext, 0, 10, fac.AddNewDirection(ppn))
    fln.SecondLimitType = 2
    fln.SecondUptoElement = plane
    s1 = fac.AddNewHybridSplit(fln, limit1, flip)
    s2 = fac.AddNewHybridSplit(s1, limit2, flip)
    gst.AppendHybridShape(s2)
    sel.Add.((fln, s1, pin, ppn, ext, edges...))
end

function main_edges(pts, angle)
    edg = fac.AddNewLineAngle.(ydr, xyp, pts[end-1:end], false, 0, 1, 1, (false, true))
    relate(edg[1].EndOffset, "$(repr(rw)) + $(repr(ri))")
    relate(edg[1].Angle, angle)
    relate(edg[2].BeginOffset, "-$(repr(ri))")
    relate(edg[2].EndOffset, "$(repr(rw))")
    relate(edg[2].Angle, angle)
    return edg
end

function ribbons()

    px = [fac.AddNewPointCoord(0, 0, 0)]
    dt = 1 / (rn-1)
    pl = nothing
    bl = nothing

    for n in 0:rn-1
        ndt = n * dt
        sgn = 1 - ((n % 2) * 2)
        ang = "$(repr(a1)) * $(ndt) + $(repr(a0)) * (1- $(ndt))"
        frm = "$(repr(rw))/cos($(ang))"
        frm =  n == 0 ? frm : "$(repr(px[end].Y)) + ($(frm))"

        push!(px, fac.AddNewPointCoord(0, 0, 0))
        relate(px[end].Z, "$((n+1) % 2) * -$(repr(rh))")
        relate(px[end].Y, frm)
        gst.AppendHybridShape(px[end])

        edges = main_edges(px[end-1:end], ang)
        gst.AppendHybridShape.(edges)

        blend = blend_surface(edges, sgn)
        gst.AppendHybridShape(blend)

        if pl != nothing
            flange_surface(edges, pl, bl, blend, sgn)
        end

        pl = fac.AddNewPlaneAngle(xyp, edges[1], 90, false)
        bl = blend
    end
end

ribbons()
sel.VisProperties.SetShow(1)
prt.Update()

max_stars_repo_name: zaika-generative/ShapeFactory

max_issues_repo_name: zaika-generative/ShapeFactory

max_forks_repo_name: zaika-generative/ShapeFactory


text:
import Base: @propagate_inbounds,getindex, setindex!,iterate,size,length,push!,
              collect,view

export BodyList, RigidMotionList, RigidTransformList, DirectlySpecifiedMotionList,
          getrange, numpts

abstract type SetOfBodies end

const LISTS = [:BodyList, :Body],
              [:RigidMotionList, :RigidBodyMotion],
              [:RigidTransformList, :RigidTransform],
              [:DirectlySpecifiedMotionList,:DirectlySpecifiedMotion]

for (listtype,listelement) in LISTS

  @eval struct $listtype <: SetOfBodies
      list :: Vector{$listelement}
  end

  @eval Base.eltype(f::$listtype) = Base.eltype(f.list)

  @eval $listtype() = $listtype($listelement[])

  @eval @propagate_inbounds getindex(A::$listtype, i::Int) = A.list[i]
  @eval @propagate_inbounds setindex!(A::$listtype, v::$listelement, i::Int) = A.list[i] = v
  @eval @propagate_inbounds getindex(A::$listtype, I...) = A.list[I...]
  @eval @propagate_inbounds setindex!(A::$listtype, v, I...) = A.list[I...] = v

  @eval iterate(A::$listtype) = iterate(A.list)
  @eval iterate(A::$listtype,I) = iterate(A.list,I)
  @eval size(A::$listtype) = size(A.list)
  @eval length(A::$listtype) = length(A.list)

  @eval push!(bl::$listtype,b::$listelement) = push!(bl.list,b)

end

numpts(bl::BodyList) = mapreduce(numpts,+,bl)

"""
    collect(bl::bodylist) -> Vector{Float64}, Vector{Float64}

Collect the inertial-space coordinates of all of the Lagrange points comprising
the bodies in body list `bl` and return each assembled set of coordinates as a vector.
"""
function collect(bl::BodyList)
    xtmp = Float64[]
    ytmp = Float64[]
    for b in bl
        append!(xtmp,b.x)
        append!(ytmp,b.y)
    end
    return xtmp,ytmp
end
collect(body::Body) = collect(BodyList([body]))

"""
    getrange(bl::BodyList,i::Int) -> Range

Return the subrange of indices in the global set of surface point data
corresponding to body `i` in a BodyList `bl`.
"""
function getrange(bl::BodyList,i::Int)
    i <= length(bl) || error("Unavailable body")
    first = 1
    j = 1
    while j < i
        first += length(bl[j])
        j += 1
    end
    last = first+length(bl[i])-1
    return first:last
end

"""
    view(f::AbstractVector,bl::BodyList,i::Int) -> SubArray

Provide a view of the range of values in vector `f` corresponding to the Lagrange
points of the body with index `i` in a BodyList `bl`.
"""
function Base.view(f::AbstractVector,bl::BodyList,i::Int)
    length(f) == numpts(bl) || error("Inconsistent size of data for viewing")
    return view(f,getrange(bl,i))
end

"""
    sum(f::AbstractVector,bl::BodyList,i::Int) -> Real

Compute a sum of the elements of vector `f` corresponding to body `i` in body
list `bl`.
"""
Base.sum(f::AbstractVector,bl::BodyList,i::Int) = sum(view(f,bl,i))

"""
    (tl::RigidTransformList)(bl::BodyList) -> BodyList

Carry out in-place transformations of each body in `bl` with the
corresponding transformation in `tl`.
"""
@inline function (tl::RigidTransformList)(bl::BodyList)
  length(tl) == length(bl) || error("Inconsistent lengths of lists")
  map((T,b) -> T(b),tl,bl)
end

"""
    vec(tl::RigidTransformList) -> Vector{Float64}

Returns a concatenation of length-3 vectors of the form [x,y,α] corresponding to the translation
and rotation specified by the given by the list of transforms `tl`.
"""
vec(tl::RigidTransformList) where {N} = mapreduce(vec,vcat,tl)

"""
    RigidTransformList(x::Vector)

Parses vector `x`, containing rigid-body configuration data ordered as
x, y, α for each body, into a `RigidTransformList`.
"""
function RigidTransformList(x::Vector{T}) where T <: Real
    nb, md = _length_and_mod(x)
    md == 0 || error("Invalid length of vector")
    first = 0
    tl = RigidTransformList()
    for i in 1:nb
        push!(tl,RigidTransform(view(x,first+1:first+CHUNK)))
        first += CHUNK
    end
    return tl
end

_length_and_mod(x::Vector{T}) where T <: Real = (n = length(x); return n ÷ CHUNK, n % CHUNK)

"""
    rigidbodyvelocity(ml::RigidMotionList,t::Real) -> Vector

Return the velocity components (as a vector) of a `RigidMotionList`
at the given time `t`.
"""
function rigidbodyvelocity(ml::Union{RigidMotionList,DirectlySpecifiedMotionList},t::Real)
    u = Float64[]
    for m in ml
      ui = rigidbodyvelocity(m,t)
      append!(u,ui)
    end
    return u
end


text:
export dubins, dubins_length, dubins_waypoints

dubins_length(q0::StaticVector{3}, qf::StaticVector{3}; r=1) = dubins(q0, qf, r=r).cost
function dubins_waypoints(q0::StaticVector{3}, qf::StaticVector{3}, dt_or_N; v=1, r=1)
    waypoints(SimpleCarDynamics{0,0}(), SE2State(q0), dubins(q0, qf, v=v, r=r).controls, dt_or_N)
end

function dubins((x0, y0, θ0)::StaticVector{3,T0}, (xf, yf, θf)::StaticVector{3,Tf}; v=1, r::R=1) where {T0,Tf,R}
    T = promote_type(T0, Tf, R)
    dx = (xf - x0)/r
    dy = (yf - y0)/r
    d = hypot(dx, dy)
    θ = atan(dy, dx)
    a = θ0 - θ
    b = θf - θ
    sa, ca = sincos(a)
    sb, cb = sincos(b)

    cmin = T(Inf)
    # LSL
    tmp = 2 + d*d - 2*(ca*cb + sa*sb - d*(sa - sb))
    if tmp >= 0
        θ = atan(cb - ca, d + sa - sb)
        t = mod2piF(-a + θ)
        p = sqrt(max(tmp, 0))
        q = mod2piF(b - θ)
        c = t + p + q

        cmin = c
        ctrl = SVector(
            carsegment2stepcontrol(1, t),
            carsegment2stepcontrol(0, p),
            carsegment2stepcontrol(1, q)
        )
    end

    # RSR
    tmp = 2 + d*d - 2*(ca*cb + sa*sb - d*(sb - sa))
    if tmp >= 0
        θ = atan(ca - cb, d - sa + sb)
        t = mod2piF(a - θ)
        p = sqrt(max(tmp, 0))
        q = mod2piF(-b + θ)
        c = t + p + q
        if c < cmin
            cmin = c
            ctrl = SVector(
                carsegment2stepcontrol(-1, t),
                carsegment2stepcontrol( 0, p),
                carsegment2stepcontrol(-1, q)
            )
        end
    end

    # RSL
    tmp = d*d - 2 + 2*(ca*cb + sa*sb - d*(sa + sb))
    if tmp >= 0
        p = sqrt(max(tmp, 0))
        θ = atan(ca + cb, d - sa - sb) - atan(T(2), p)
        t = mod2piF(a - θ)
        q = mod2piF(b - θ)
        c = t + p + q
        if c < cmin
            cmin = c
            ctrl = SVector(
                carsegment2stepcontrol(-1, t),
                carsegment2stepcontrol( 0, p),
                carsegment2stepcontrol( 1, q)
            )
        end
    end

    # LSR
    tmp = -2 + d*d + 2*(ca*cb + sa*sb + d*(sa + sb))
    if tmp >= 0
        p = sqrt(max(tmp, 0))
        θ = atan(-ca - cb, d + sa + sb) - atan(-T(2), p)
        t = mod2piF(-a + θ)
        q = mod2piF(-b + θ)
        c = t + p + q
        if c < cmin
            cmin = c
            ctrl = SVector(
                carsegment2stepcontrol( 1, t),
                carsegment2stepcontrol( 0, p),
                carsegment2stepcontrol(-1, q)
            )
        end
    end

    # RLR
    tmp = (6 - d*d  + 2*(ca*cb + sa*sb + d*(sa - sb)))/8
    if abs(tmp) < 1
        p = 2*T(pi) - acos(tmp)
        θ = atan(ca - cb, d - sa + sb)
        t = mod2piF(a - θ + p/2)
        q = mod2piF(a - b - t + p)
        c = t + p + q
        if c < cmin
            cmin = c
            ctrl = SVector(
                carsegment2stepcontrol(-1, t),
                carsegment2stepcontrol( 1, p),
                carsegment2stepcontrol(-1, q)
            )
        end
    end

    # LRL
    tmp = (6 - d*d  + 2*(ca*cb + sa*sb - d*(sa - sb)))/8
    if abs(tmp) < 1
        p = 2*T(pi) - acos(tmp)
        θ = atan(-ca + cb, d + sa - sb)
        t = mod2piF(-a + θ + p/2)
        q = mod2piF(b - a - t + p)
        c = t + p + q
        if c < cmin
            cmin = c
            ctrl = SVector(
                carsegment2stepcontrol( 1, t),
                carsegment2stepcontrol(-1, p),
                carsegment2stepcontrol( 1, q)
            )
        end
    end

    ctrl = scalespeed.(scaleradius.(ctrl, r), v)
    (cost=cmin*r/v, controls=ctrl)
end


text:
@testset "RDA IO, single voxel" begin
    header,data = load_rda("rda/single_voxel_example.rda")

    # Basic header reading
    @test header["PatientName"] == "Foo"
    @test header["PatientBirthDate"] == "19700101"
    @test header["TR"] == "2000.000000"
    @test header["PixelSpacing3D"] == "20.000000"
    @test length(header) == 83

    fake_fid = ComplexF64.(1:1024, 1024:-1:1)
    @test size(data) == size(fake_fid)
    @test data == fake_fid

    # Test exact roundtrip through save_rda
    io = IOBuffer()
    save_rda(io, header, data)
    @test read("rda/single_voxel_example.rda") == take!(io)
end

@testset "RDA IO, CSI" begin
    header,data = load_rda("rda/CSI_example.rda")

    # Check that data came through correctly shaped
    fake_csi_data = reshape(ComplexF64.(1:64*2*3*1, 64*2*3*1:-1:1), 64,2,3,1)
    @test size(data) == size(fake_csi_data)
    @test data == fake_csi_data

    # Test exact roundtrip through save_rda
    io = IOBuffer()
    save_rda(io, header, data)
    @test read("rda/CSI_example.rda") == take!(io)
end

max_stars_repo_name: TRIImaging/MagneticResonanceSignals

max_issues_repo_name: TRIImaging/MagneticResonanceSignals

max_forks_repo_name: TRIImaging/MagneticResonanceSignals


text:
"""
```
function old_to_new_cloud(cloud::DSGE.Cloud)::SMC.Cloud
```
Converter between types.
"""
function old_to_new_cloud(cloud::DSGE.Cloud)::SMC.Cloud
    return SMC.Cloud(cloud.particles, cloud.tempering_schedule,
                     cloud.ESS, cloud.stage_index, cloud.n_Φ, cloud.resamples,
                     cloud.c, cloud.accept, cloud.total_sampling_time)
end

"""
```
function Cloud(c::ParticleCloud)
function Cloud(c::Cloud)
```
Returns a Cloud type.
"""
function SMC.Cloud(cloud::ParticleCloud)
    return SMC.Cloud(hcat(Matrix{Float64}(get_vals(cloud)'), get_loglh(cloud), get_logprior(cloud),
                      get_old_loglh(cloud), get_accept(cloud), get_weights(cloud)),
                 cloud.tempering_schedule, cloud.ESS, cloud.stage_index, cloud.n_Φ,
                 cloud.resamples, cloud.c, cloud.accept, cloud.total_sampling_time)
end
function SMC.Cloud(cloud::SMC.Cloud)
    return cloud
end

function DSGE.Cloud(cloud::SMC.Cloud)
    return DSGE.Cloud(cloud.particles, cloud.tempering_schedule,
                     cloud.ESS, cloud.stage_index, cloud.n_Φ, cloud.resamples,
                     cloud.c, cloud.accept, cloud.total_sampling_time)
end
function SMC.Cloud(cloud::DSGE.Cloud)
    return SMC.Cloud(cloud.particles, cloud.tempering_schedule,
                     cloud.ESS, cloud.stage_index, cloud.n_Φ, cloud.resamples,
                     cloud.c, cloud.accept, cloud.total_sampling_time)
end


text:
# Sequential preprocessor for large sized text files
# Gather the sequences of the same length s.t. when it is called it brings batches of similar sequences with the help of iterables

import Base: start, next, done

const SOS = "<s>"
const EOS = "</s>"
const UNK = "<unk>"
const LLIMIT = 2

type Data
    word_to_index::Dict{AbstractString, Int}
    index_to_word::Vector{AbstractString}
    batchsize::Int
    serve_type::AbstractString
    sequences::Dict
end

function Data(datafile; word_to_index=nothing, vocabfile=nothing, serve_type="bitarray", batchsize=20)
    # TODO: right now Data type can only be a bit array because of Knet support, no need to following line
    # (serve_type == nothing) && error("Please specify the data serve type: onehot, bitarray or sequence")
    existing_vocab = (word_to_index != nothing)
    if !existing_vocab
        word_to_index = Dict{AbstractString, Int}(SOS=>1, EOS=>2, UNK=>3)
        vocabfile != nothing && info("Working with provided vocabfile : $vocabfile")
        vocabfile != nothing && (V = vocab_from_file(vocabfile))
    end

    stream = open(datafile)
    sequences = Dict{Int64, Array{Any, 1}}()
    for line in eachline(stream)
        words = Int32[]
        push!(words, word_to_index[SOS])
        for word in split(line)
            if !existing_vocab && vocabfile != nothing && !(word in V)
                word = UNK
            end
            if existing_vocab
                index = get(word_to_index, word, word_to_index[UNK])
            else
                index = get!(word_to_index, word, 1+length(word_to_index))
            end
            push!(words, index)
        end
        push!(words, word_to_index[EOS])
        skey = length(words)
        (skey == LLIMIT) && println("LLMIT needs to be checked") #'ceause we already put <s> and </s> tokens?
        (!haskey(sequences, skey)) && (sequences[skey] = Any[])
        push!(sequences[skey], words)
    end
    close(stream)
    vocabsize = length(word_to_index)
    index_to_word = Array(AbstractString, vocabsize)
    for (word, index) in word_to_index
        index_to_word[index] = word
    end
    Data(word_to_index, index_to_word, batchsize, serve_type, sequences)
end

function sentenbatch(nom::Array{Any,1}, from::Int, batchsize::Int, vocabsize::Int; serve_type="bitarr")
    total = length(nom)
    to = (from + batchsize - 1 < total) ? (from + batchsize - 1) : total

    # not to work with surplus sentences
    if (to-from + 1 < batchsize)
        warning("Surplus does not being cleaned correctly!")
        return (nothing, 1)
    end

    new_from = (to == total) ? 1 : (to + 1)
    seqlen = length(nom[1]) # TODO: get rid of length computation give it as an extra argument!
    sentences = nom[from:to]

    # If Knet supports lookup that part will work correctly, it is commented for speed purposes
    #if serve_type == "lookup"
    #    return (sentences, new_from)
    #end

    scount = batchsize # modified future code
    data = [ falses(scount, vocabsize) for i=1:seqlen ]
    for cursor=1:seqlen
        for row=1:scount
            index = sentences[row][cursor]
            data[cursor][row, index] = 1
        end
    end
    return (data, new_from)
end

"""Removes the surplus sentences randomly"""
function clean_seqdict!(seqdict::Dict{Int64,Array{Any,1}}, batchsize::Int)
    for seqlen in keys(seqdict)
        remain = rem(length(seqdict[seqlen]), batchsize)
        while remain != 0
            index = rand(1:length(seqdict[seqlen]))
            deleteat!(seqdict[seqlen], index)
            remain -= 1
        end
        if isempty(seqdict[seqlen])
            delete!(seqdict, seqlen)
        end
    end
end

function start(s::Data)
    sdict = deepcopy(s.sequences)
    clean_seqdict!(sdict, s.batchsize)
    @assert (!isempty(sdict)) "There is not enough data with that batchsize $(s.batchsize)"
    slens = collect(keys(sdict))
    seqlen = pop!(slens)
    from = nothing
    vocabsize = length(s.word_to_index)
    state = (sdict, seqlen, slens,from, vocabsize)
    return state
end

function next(s::Data, state)
    (sdict, seqlen, slens, from, vocabsize) = state

    if from == nothing
        (item, new_from) = sentenbatch(sdict[seqlen], 1, s.batchsize, vocabsize)
    elseif from == 1
        seqlen = pop!(slens)
        (item, new_from) = sentenbatch(sdict[seqlen], from, s.batchsize, vocabsize)
    else
        (item, new_from) = sentenbatch(sdict[seqlen], from, s.batchsize, vocabsize)
    end
    from = new_from
    state = (sdict, seqlen, slens, from, vocabsize)
    return (item, state)
end

function done(s::Data, state)
    (sdict, seqlen, slens, from, vocabsize) = state
    return isempty(slens) && (from == 1)
end

""" Creates a set that contains all the words in that file, vocab file given as each vocab in a single line
    sorted_counted represents pure create_vocab.sh output
"""
function vocab_from_file(vocabfile; sorted_counted=false)
    V = Set{AbstractString}()
    open(vocabfile) do file
        for line in eachline(file)
            line = split(line)
            if sorted_counted
                (length(line)>1) && push!(V, line[2])
            else
                !isempty(line) && push!(V, line[1])
            end
        end
    end
    return V
end

""" Builds the kth sentence from a given sequence """
function ibuild_sentence(tdata::Data, sequence::Array{BitArray{2},1}, kth::Int)
    sentence = Array{Any, 1}()
    for i=1:length(sequence)
        z = find(x->x==true, sequence[i][kth,:])
        @assert(length(z) == 1)
        append!(sentence, z)
    end
    ret = map(x->tdata.index_to_word[x], sentence)
    return ret
end

# FUTURE CODE: if one day knet8 allows us to change batchsize on the fly, following lines will implement surplus batch implementation, this code snippet would be put on sentenbatch
    # (length(sentences) != batchsize) && (println("I am using the surplus sentences:) $from : $to"))
    # scount = length(sentences) # it can be either batchsize or the surplus sentences
    # data = [ falses(scount, vocabsize) for i=1:seqlen ]
    # for cursor=1:seqlen
    #     for row=1:scount
    #         index = sentences[row][cursor]
    #         data[cursor][row, index] = 1
    #     end
    # end


text:
###########################################################################
#
#  Chinese remaindering
#
################################################################################

@doc Markdown.doc"""
    induce_crt(a::Generic.Poly{nf_elem}, p::fmpz, b::Generic.Poly{nf_elem}, q::fmpz) -> Generic.Poly{nf_elem}, fmpz

Given polynomials $a$ defined modulo $p$ and $b$ modulo $q$, apply the CRT
to all coefficients recursively.
Implicitly assumes that $a$ and $b$ have integral coefficients (ie. no
denominators).
"""
function induce_crt(a::Generic.Poly{nf_elem}, p::fmpz, b::Generic.Poly{nf_elem}, q::fmpz, signed::Bool = false)
  c = parent(a)()
  pi = invmod(p, q)
  mul!(pi, pi, p)
  pq = p*q
  if signed
    pq2 = div(pq, 2)
  else
    pq2 = fmpz(0)
  end
  for i=0:max(degree(a), degree(b))
    setcoeff!(c, i, induce_inner_crt(coeff(a, i), coeff(b, i), pi, pq, pq2))
  end
  return c, pq
end

@doc Markdown.doc"""
  induce_rational_reconstruction(a::Generic.Poly{nf_elem}, M::fmpz) -> bool, Generic.Poly{nf_elem}

Apply rational reconstruction to the coefficients of $a$. Implicitly assumes
the coefficients to be integral (no checks done)
returns true iff this is successful for all coefficients.
"""
function induce_rational_reconstruction(a::Generic.Poly{nf_elem}, M::fmpz)
  b = parent(a)()
  for i=0:degree(a)
    fl, x = rational_reconstruction(coeff(a, i), M)
    if fl
      setcoeff!(b, i, x)
    else
      return false, b
    end
  end
  return true, b
end

################################################################################
#
#  GCD
#
################################################################################

@doc Markdown.doc"""
  gcd(a::Generic.Poly{nf_elem}, b::Generic.Poly{nf_elem}) -> Generic.Poly{nf_elem}

Computes the greatest common divisor of $f$ and $g$ using a modular algorithm.
"""
function gcd(a::Generic.Poly{nf_elem}, b::Generic.Poly{nf_elem}, test_sqfr::Bool = false)
  # modular kronnecker assumes a, b !=n 0
  if iszero(a)
    if iszero(b)
      return b
    else
      return  inv(lead(b))*b
    end
  elseif iszero(b)
    return inv(lead(a))*a
  end
  if min(degree(a), degree(b)) >= 6 || degree(base_ring(a)) > 5 || test_sqfr
    g = gcd_modular_kronnecker(a, b, test_sqfr)
    test_sqfr && return g
    return inv(lead(g))*g  # we want it monic...
  else
    return gcd_euclid(a, b)
  end
end

################################################################################
#
#  Modular GCD
#
################################################################################

# There is some weird type instability
function gcd_modular(a::Generic.Poly{nf_elem}, b::Generic.Poly{nf_elem})
  # naive version, kind of
  # polys should be integral
  # rat recon maybe replace by known den if poly integral (Kronnecker)
  # if not monic, scale by gcd
  # remove content?
  a = a*(1//leading_coefficient(a))
  b = b*(1//leading_coefficient(b))
  global p_start
  p = p_start
  K = base_ring(parent(a))
  @assert parent(a) == parent(b)
  g = zero(a)
  d = fmpz(1)
  while true
    p = next_prime(p)
    me = modular_init(K, p)
    t = Hecke.modular_proj(a, me)
    fp = deepcopy(t)::Array{fq_nmod_poly, 1}  # bad!!!
    gp = Hecke.modular_proj(b, me)
    gp = [gcd(fp[i], gp[i]) for i=1:length(gp)]::Array{fq_nmod_poly, 1}
    gc = Hecke.modular_lift(gp, me)::Generic.Poly{nf_elem}
    if isone(gc)
      return parent(a)(1)
    end
    if d == 1
      g = gc
      d = fmpz(p)
    else
      if degree(gc) < degree(g)
        g = gc
        d = fmpz(p)
      elseif degree(gc) > degree(g)
        continue
      else
        g, d = induce_crt(g, d, gc, fmpz(p))
      end
    end
    fl, gg = induce_rational_reconstruction(g, d)
    if fl  # not optimal
      r = mod(a, gg)
      if iszero(r)
        r = mod(b, gg)
        if iszero(r)
          return gg
        end
      end
    end
  end
end

function _preproc_pol(a::Generic.Poly{nf_elem}, b::Generic.Poly{nf_elem})
  a1 = a*(1//leading_coefficient(a))
  da = Base.reduce(lcm, [denominator(coeff(a1, i)) for i=0:degree(a)])
  b1 = b*(1//leading_coefficient(b))
  db = Base.reduce(lcm, [denominator(coeff(b1, i)) for i=0:degree(b)])
  d = gcd(da, db)
  a2 = a1*da
  b2 = b1*db
  Kt = parent(a)
  K = base_ring(Kt)
  if isdefining_polynomial_nice(K)
    fsa = evaluate(derivative(K.pol), gen(K))*d
  else
    fsa = short_elem(different(any_order(K)))*d
  end
  return a2, b2, fsa
end

function gcd_euclid(a::AbstractAlgebra.PolyElem{nf_elem}, b::AbstractAlgebra.PolyElem{nf_elem})
   check_parent(a, b)
   if length(a) > length(b)
      (a, b) = (b, a)
   end
   while !iszero(a)
      (a, b) = (mod(b, a), a)
   end
   d = lead(b)
   return divexact(b, d)
end

#similar to gcd_modular, but avoids rational reconstruction by controlling
#a/the denominator
function gcd_modular_kronnecker(a::Generic.Poly{nf_elem}, b::Generic.Poly{nf_elem}, test_sqfr::Bool = false)
  # rat recon maybe replace by known den if poly integral (Kronnecker)
  # if not monic, scale by gcd
  # remove content?
  @assert parent(a) == parent(b)
  a, b, fsa = _preproc_pol(a, b)
  Kt = parent(a)
  K = base_ring(Kt)
  #now gcd(a, b)*fsa should be in the equation order...
  global p_start
  p = p_start

  g = zero(Kt)
  d = fmpz(1)
  last_g = zero(Kt)
  while true
    p = next_prime(p)
    F = GF(p, cached = false)
    Fx = PolynomialRing(F, "x", cached = false)[1]
    Fp = Fx(K.pol)
    if !issquarefree(Fp)
      continue
    end
    fs = factor_shape(Fp)
    if any(x -> x > 4, keys(fs))
      continue
    end
    local me, fp, gp, fsap
    try
      me = modular_init(K, p)
      fp = deepcopy(Hecke.modular_proj(a, me))  # bad!!!
      gp = Hecke.modular_proj(b, me)
      fsap = Hecke.modular_proj(fsa, me)
    catch ee
      if typeof(ee) != Hecke.BadPrime
        rethrow(ee)
      end
      continue
    end
    gp = fq_nmod_poly[fsap[i] * gcd(fp[i], gp[i]) for i=1:length(gp)]
    gc = Hecke.modular_lift(gp, me)
    if isconstant(gc)
      return one(Kt)
    end
    if test_sqfr
      return zero(Kt)
    end
    if isone(d)
      g = gc
      d = fmpz(p)
    else
      if degree(gc) < degree(g)
        g = gc
        d = fmpz(p)
      elseif degree(gc) > degree(g)
        continue
      else
        g, d = induce_crt(g, d, gc, fmpz(p), true)
      end
    end
    if g == last_g && iszero(mod(a, g)) && iszero(mod(b, g))
      return divexact(g, lead(g))
    else
      last_g = g
    end
  end
end

################################################################################
#
#  Modular extended GCD
#
################################################################################

#seems to be faster than gcdx - if problem large enough.
#rational reconstructio is expensive - enventually
#TODO: figure out the denominators in advance. Resultants?

function gcdx_modular(a::Generic.Poly{nf_elem}, b::Generic.Poly{nf_elem})
  a = a*(1//leading_coefficient(a))
  b = b*(1//leading_coefficient(b))
  global p_start
  p = p_start
  K = base_ring(parent(a))
  @assert parent(a) == parent(b)
  g = zero(a)
  d = fmpz(1)
  last_g = parent(a)(0)
  while true
    p = next_prime(p)
    me = modular_init(K, p)
    fp = deepcopy(Hecke.modular_proj(a, me))  # bad!!!
    gp = Hecke.modular_proj(b, me)
    ap = similar(gp)
    bp = similar(gp)
    for i=1:length(gp)
      gp[i], ap[i], bp[i] = gcdx(fp[i], gp[i])
    end
    gc = Hecke.modular_lift(gp, me)
    aa = Hecke.modular_lift(ap, me)
    bb = Hecke.modular_lift(bp, me)
    if d == 1
      g = gc
      ca = aa
      cb = bb
      d = fmpz(p)
    else
      if degree(gc) < degree(g)
        g = gc
        ca = aa
        cb = bb
        d = fmpz(p)
      elseif degree(gc) > degree(g)
        continue
      else
        g, dd = induce_crt(g, d, gc, fmpz(p))
        ca, dd = induce_crt(ca, d, aa, fmpz(p))
        cb, d = induce_crt(cb, d, bb, fmpz(p))
      end
    end
    fl, ccb = Hecke.induce_rational_reconstruction(cb, d)
    if fl
      fl, cca = Hecke.induce_rational_reconstruction(ca, d)
    end
    if fl
      fl, gg = Hecke.induce_rational_reconstruction(g, d)
    end
    if fl
      r = mod(a, g)
      if iszero(r)
        r = mod(b, g)
        if iszero(r) && ((cca*a + ccb*b) == gg)
          return gg, cca, ccb
        end
      end
    end
  end
end

#similar to gcd_modular, but avoids rational reconstruction by controlling
#a/the denominator using resultant. Faster than above, but still slow.
#mainly due to the generic resultant. Maybe use only deg-1-primes??
#fact: g= gcd(a, b) and 1= gcd(a/g, b/g) = u*(a/g) + v*(b/g)
#then u*res(a/g, b/g) is mathematically integeral, same for v
#scaling by f'(a) makes it i nthe equation order
#
# missing/ next attempt:
#  write invmod using lifting
#  write gcdx using lifting (lin/ quad)
#  try using deg-1-primes only (& complicated lifting)
#
function gcdx_mod_res(a::Generic.Poly{nf_elem}, b::Generic.Poly{nf_elem})
  @assert parent(a) == parent(b)
  a = a*(1//leading_coefficient(a))
  da = Base.reduce(lcm, [denominator(coeff(a, i)) for i=0:degree(a)])
  b = b*(1//leading_coefficient(b))
  db = Base.reduce(lcm, [denominator(coeff(a, i)) for i=0:degree(a)])
  d = gcd(da, db)
  a = a*da
  b = b*db
  Kt = parent(a)
  K = base_ring(Kt)
  fsa = change_base_ring(K, derivative(K.pol), parent = Kt)*d
  #now gcd(a, b)*fsa should be in the equation order...
  global p_start
  p = p_start
  g = zero(parent(a))
  d = fmpz(1)
  r = zero(K)
  fa = zero(parent(a))
  fb = zero(parent(b))
  last_g = (parent(a)(0), parent(a)(0), parent(a)(0), parent(a)(0))

  while true
    p = next_prime(p)
    me = modular_init(K, p)
    fp = deepcopy(Hecke.modular_proj(a, me))  # bad!!!
    gp = (Hecke.modular_proj(b, me))
    fsap = (Hecke.modular_proj(fsa, me))
    g_ = similar(fp)
    fa_ = similar(fp)
    fb_ = similar(fp)
    r_ = Array{fq_nmod}(length(fp))
    for i=1:length(gp)
      g_[i], fa_[i], fb_[i] = gcdx(fp[i], gp[i])
      r_[i] = resultant(div(fp[i], g_[i]), div(gp[i], g_[i]))
      g_[i] *= fsap[i]
      fa_[i] *= (fsap[i]*r_[i])
      fb_[i] *= (fsap[i]*r_[i])
    end
    rc = Hecke.modular_lift(r_, me)
    gc = Hecke.modular_lift(g_, me)
    fac = Hecke.modular_lift(fa_, me)
    fbc = Hecke.modular_lift(fb_, me)
    if d == 1
      g = gc
      r = rc
      fa = fac
      fb = fbc
      d = fmpz(p)
    else
      if degree(gc) < degree(g)
        g = gc
        r = rc
        fa = fac
        fb = fbc
        d = fmpz(p)
      elseif degree(gc) > degree(g)
        continue
      else
        g, d1 = induce_crt(g, d, gc, fmpz(p), true)
        fa, d1 = induce_crt(fa, d, fac, fmpz(p), true)

        r = Hecke.induce_inner_crt(r, rc, d*invmod(d, fmpz(p)), d1, div(d1, 2))
        fb, d = induce_crt(fb, d, fbc, fmpz(p), true)

      end
    end
    if (g, r, fa, fb) == last_g
      if g*r == fa*a + fb*b
        return g*r, fa, fb ## or normalise to make gcd monic??
      else
        last_g = (g, r, fa, fb)
      end
    else
      last_g = (g, r, fa, fb)
    end
  end
end

################################################################################
#
#  Equality modulo rational integer
#
################################################################################

function eq_mod(a::Generic.Poly{nf_elem}, b::Generic.Poly{nf_elem}, d::fmpz)
  e = degree(a) == degree(b)
  K= base_ring(parent(a))
  i=0
  while e && i<= degree(a)
    j = 0
    while e && j<degree(K)
      e = e && (numerator(coeff(coeff(a, i), j)) - numerator(coeff(coeff(b, i), j))) % d == 0
      j += 1
    end
    i += 1
  end
  return e
end

################################################################################
#
#  Predicates
#
################################################################################

function issquarefree(x::Generic.Poly{nf_elem})
  return isone(gcd(x, derivative(x), true))
end

################################################################################
#
#  Some weird conversion
#
################################################################################

function nf_poly_to_xy(f::PolyElem{Nemo.nf_elem}, Qxy::PolyRing, Qx::PolyRing)
  K = base_ring(f)
  Qy = parent(K.pol)
  y = gen(Qx)

  res = zero(Qxy)
  for i=degree(f):-1:0
    res *= y
    res += change_base_ring(Qx, Qy(coeff(f, i)), parent = Qxy)
  end
  return res
end

################################################################################
#
#  Modular resultant
#
################################################################################

function resultant_mod(f::Generic.Poly{nf_elem}, g::Generic.Poly{nf_elem})
  global p_start
  p = p_start
  K = base_ring(parent(f))
  @assert parent(f) == parent(g)
  r = K()
  d = fmpz(1)
  last_r = K()
  first = true
  while true
    p = next_prime(p)
    me = modular_init(K, p)
    fp = deepcopy(Hecke.modular_proj(f, me))  # bad!!!
    gp = Hecke.modular_proj(g, me)
    rp = Array{fq_nmod}(undef, length(gp))
    for i=1:length(gp)
      rp[i] = resultant(fp[i], gp[i])
    end
    rc = Hecke.modular_lift(rp, me)
    if d == 1
      r = rc
      d = fmpz(p)
    else
      r, d = induce_crt(r, d, rc, fmpz(p))
    end
    fl, ccb = Hecke.rational_reconstruction(r, d)
    if fl
      if first
        first = false
        last_r = ccb
#        println("first: $ccb")
      else
        if ccb == last_r
          return ccb
        else
#      println("fail2: $ccb")
          last_r = ccb
        end
      end
    else
#      println("fail")
    end
  end
end

function landau_mignotte_bound(f::PolyElem{nf_elem})
  Zx, x = PolynomialRing(FlintZZ, cached = false)
  g = Zx()
  for i=0:degree(f)
    setcoeff!(g, i, Hecke.upper_bound(sqrt(t2(coeff(f, i))), fmpz))
  end
  b = fmpz()
  ccall((:fmpz_poly_factor_mignotte, libflint), Nothing, (Ref{fmpz}, Ref{fmpz_poly}), b, g)
  return b
end

function cld_bound(f::PolyElem{nf_elem}, k::Array{Int, 1})
  @assert all(kk -> 0 <= kk < degree(f), k)
  Zx, x = PolynomialRing(FlintZZ, cached = false)
  g = Zx()
  for i=0:degree(f)
    setcoeff!(g, i, Hecke.upper_bound(sqrt(t2(coeff(f, i))), fmpz))
  end
  bb = fmpz[]
  for kk = k
    b = FlintZZ()
    ccall((:fmpz_poly_CLD_bound, libflint), Nothing, (Ref{fmpz}, Ref{fmpz_poly}, Int64), b, g, kk)
    push!(bb, b)
  end
  return bb
end
cld_bound(f::PolyElem{nf_elem}, k::Int) = cld_bound(f, [k])[1]

function cld_bound(f::fmpz_poly, k::Int)
  @assert 0 <= k < degree(f)
  b = FlintZZ()
  ccall((:fmpz_poly_CLD_bound, libflint), Nothing, (Ref{fmpz}, Ref{fmpz_poly}, Int64), b, f, k)
  return b
end
cld_bound(f::fmpz_poly, k::Array{Int, 1}) = map(x->cld_bound(f, x), k)


text:
module MaxHelpingHandNorthernLightsCustomColors

using ..Ahorn, Maple

@mapdef Effect "MaxHelpingHand/NorthernLightsCustomColors" NorthernLightsCustomColors(only::String="*", exclude::String="", gradientColor1::String="020825", gradientColor2::String="170c2f",
    colors::String="2de079,62f4f6,45bc2e,3856f0", displayBackground::Bool=true)

placements = NorthernLightsCustomColors

function Ahorn.canFgBg(effect::NorthernLightsCustomColors)
    return true, true
end

end

max_stars_repo_name: Kitty-Cats/MaxHelpingHand

max_issues_repo_name: Kitty-Cats/MaxHelpingHand

max_forks_repo_name: Kitty-Cats/MaxHelpingHand


text:
# This file is part of GenericSchur.jl, released under the MIT "Expat" license

struct Balancer{T}
    ilo :: Int
    ihi :: Int
    prow :: Vector{Int}
    pcol :: Vector{Int}
    D :: Vector{T}
    trivial :: Bool
end

# intended to be a translation of xGEBAL from LAPACK
# Copyright:
# Univ. of Tennessee
# Univ. of California Berkeley
# Univ. of Colorado Denver
# NAG Ltd.
# w/ some allowance for scheme from R.James et al. 2014 arxiv 1401.5766
"""
    balance!(A; scale=true, permute=true) => Abal, B::Balancer

Balance a matrix so that various operations are more stable.

If `permute`, then row and column permutations are found such
that `Abal` has the block form `[T₁ X Y; 0 C Z; 0 0 T₂]` where
`T₁` and `T₂` are upper-triangular.
If `scale`, then a diagonal similarity transform using powers of 2
is found such that the 1-norm of `Abal` is near unity. The transformations
are encoded into `B` so that they can be inverted for eigenvectors, etc.
Balancing typically improves the accuracy of eigen-analysis.
"""
function balance!(A::AbstractMatrix{T}; scale=true, permute=true,
                  p::Int=1, algo = :pr) where T
    n = checksquare(A)
    RT = real(T)
    β = RT(2)
    factor = RT(0.95)

    sfmin1 = floatmin(RT) / eps(RT)
    sfmax1 = one(RT) / sfmin1
    sfmin2 = sfmin1 * β
    sfmax2 = one(RT) / sfmin2

    D = ones(T,n)
    ilo = 1
    ihi = n
    trivial = true
    sp = fill(0,n)

    if permute
        # look for row permutations which make the lower portion UT
        ihi = n+1
        js = 0
        while ihi > 1
            ihi -= 1
            exch = false
            ms = 0
            for j=ihi:-1:1
                exch = true
                js = j
                for i=1:ihi
                    (i == j) && continue
                    if A[j,i] != 0
                        exch = false
                    end
                end
                if exch
                    ms = ihi
                    break
                end
            end
            if exch
                sp[ms] = js
                if js != ms
                    trivial = false
                    for i=1:ihi
                        A[i,js], A[i,ms] = A[i,ms], A[i,js]
                    end
                    for i=ilo:n
                        A[js,i], A[ms,i] = A[ms,i], A[js,i]
                    end
                end
            else
                break
            end
        end
        # look for column permutations which make the left portion UT
        ilo = 0
        while ilo < n
            ilo += 1
            exch = false
            ms = 0
            js = 0
            for j=ilo:ihi
                js = j
                exch = true
                for i=ilo:ihi
                    (i == j) && continue
                    if A[i,j] != 0
                        exch = false
                    end
                end
                if exch
                    ms = ilo
                    break
                end
            end
            if exch
                sp[ms] = js
                if ms != js
                    trivial = false
                    for i=1:ihi
                        A[i,js], A[i,ms] = A[i,ms], A[i,js]
                    end
                    for i=ilo:n
                        A[js,i], A[ms,i] = A[ms,i], A[js,i]
                    end
                end
            else
                break
            end
        end
    end

    sf(c,r) = (p == 1) ? (c + r) : c*c + r*r
    # abs(c) + abs(r), but args always nonneg

    if scale
        converged = false
        while !converged
            converged = true
            for i=ilo:ihi
                if algo == :pr
                    # LAPACK mixes norms
                    c = norm(view(A,ilo:ihi,i),2)
                    r = norm(view(A,i,ilo:ihi),2)
#                    c = cxnorm(A,ilo,ihi,i,p)
#                    r = rxnorm(A,ilo,ihi,i,p)
                else
                    c = norm(view(A,ilo:ihi,i),p)
                    r = norm(view(A,i,ilo:ihi),p)
                end
                ca,_ = _findamax(view(A,ilo:ihi,i))
                ra,_ = _findamax(view(A,i,ilo:ihi))
                if c == 0 || r == 0
                    continue
                end
                g = r / β
                s = sf(c,r)
                f = one(RT)
                while c < r/β
                    if c >= g || (max(f,c,ca) >= sfmax2) || (min(r,g,ra) <= sfmin2)
                        break
                    end
                    if isnan(c+f+ca+r+g+ra)
                        error("NaN encountered while balancing")
                    end
                    f *= β
                    c *= β
                    ca *= β
                    r /= β
                    g /= β
                    ra /= β
                end
                g = c / β
                while r <= c/β
                    if (g < r) || (max(r,ra) >= sfmax2) || (min(f,c,g,ca) <= sfmin2)
                        break
                    end
                    f /= β
                    c /= β
                    g /= β
                    ca /= β
                    r *= β
                    ra *= β
                end
                if f != one(RT)
                    trivial = false
                end
                (sf(c,r) >= factor * s) && continue
                # TODO: continue on [uo]flow
                converged = false
                D[i] *= f
                A[i,ilo:n] .*= (one(T) / f)
                A[1:ihi,i] .*= f
            end
        end # while !converged
    end # if scale
    prow = sp[1:ilo-1]
    pcol = sp[ihi+1:n]
    B = Balancer(ilo,ihi,prow,pcol,D,trivial)
    return A, B
end

# this is the appropriate transformation for right eigenvectors
# stored as columns of a matrix
function lmul!(B::Balancer, A::StridedMatrix{T}) where T
    n = checksquare(A)
    if B.trivial
        return A
    end
    ilo, ihi = B.ilo, B.ihi
    if ilo != ihi
        lmul!(Diagonal(B.D), A)
    end
    for j=ilo-1:-1:1
        m = B.prow[j]
        if m == j
            continue
        end
        for i=1:n
            A[j,i], A[m,i] = A[m,i], A[j,i]
        end
    end
    for j=ihi+1:n
        m = B.pcol[j-ihi]
        if m == j
            continue
        end
        for i=1:n
            A[j,i], A[m,i] = A[m,i], A[j,i]
        end
    end
    A
end

# this should be the appropriate transformation for left eigenvectors
# stored as columns of a matrix
function ldiv!(B::Balancer, A::StridedMatrix{T}) where T
    n = checksquare(A)
    if B.trivial
        return A
    end
    ilo, ihi = B.ilo, B.ihi
    if ilo != ihi
        rmul!(A, Diagonal(one(T) ./ B.D))
    end
    for j=ilo-1:-1:1
        m = B.prow[j]
        if m == j
            continue
        end
        for i=1:n
            A[j,i], A[m,i] = A[m,i], A[j,i]
        end
    end
    for j=ihi+1:n
        m = B.pcol[j-ihi]
        if m == j
            continue
        end
        for i=1:n
            A[j,i], A[m,i] = A[m,i], A[j,i]
        end
    end
    A
end

# James et al. claim Parlet & Reinsch advocate skipping diagonal elts
# but this is (no longer?) in LAPACK
function cxnorm(A::AbstractMatrix{T},i0,i1,i,p) where T
    RT = real(T)
        s = zero(RT)
        for j=i0:i1
            if j==i
                continue
            end
            if p == 1
                s += abs(A[j,i])
            else
                s += A[j,i]' * A[j,i]
            end
        end
        if p == 2
            s = sqrt(s)
        end
        return s
end
function rxnorm(A::AbstractMatrix{T},i0,i1,i,p) where T
    RT = real(T)
        s = zero(RT)
        for j=i0:i1
            if j==i
                continue
            end
            if p == 1
                s += abs(A[i,j])
            else
                s += A[i,j]' * A[i,j]
            end
        end
        if p == 2
            s = sqrt(s)
        end
        return s
end


text:
"""
Common exceptions
https://github.com/frictionlessdata/DataPackage.jl#exceptions
"""

struct PackageError <: Exception
    message::String
end

max_issues_repo_name: frictionlessdata/datapackage-jl

max_forks_repo_name: frictionlessdata/datapackage-jl


text:
module TestDay05

import AOC2021.Day05.part1, AOC2021.Day05.part2
using Test

test_input_raw = """
0,9 -> 5,9
8,0 -> 0,8
9,4 -> 3,4
2,2 -> 2,1
7,0 -> 7,4
6,4 -> 2,0
0,9 -> 2,9
3,4 -> 1,4
0,0 -> 8,8
5,5 -> 8,2"""

test_input = [string(x) for x in split(test_input_raw, "\n")]

@testset "Day 05" begin
    @testset "part 1" begin
        @test part1(test_input) == 5
    end
    @testset "part 2" begin
        @test part2(test_input) == 12
    end
end

end # module


text:
immutable Tagger <: Model; forw; back; pred; params;
    function Tagger(forw::Function, back::Function, pred::Function; o...)
        forw = Net(forw; o...)
        back = Net(back; o...)
        pred = Net(pred; o..., ninputs=2)
        par = vcat(params(forw), params(back), params(pred))
        new(forw, back, pred, par)
    end
end

params(m::Tagger)=m.params
reset!(m::Tagger)=(reset!(m.forw);reset!(m.back);reset!(m.pred))
test(m::Tagger, data, loss; o...)=(l=zeros(2); tagger_loop(m, data, loss; losscnt=l, o...); l[1]/l[2])
train(m::Tagger, data, loss; o...)=tagger_loop(m, data, loss; trn=true, o...)

function tagger_loop(m::Tagger, data, loss; gcheck=false, o...)                     # t3:7904
    x,ygold,mask = Any[],Any[],Any[]
    for item in data                                                            # each item contains one (x,y,mask) triple for a minibatch of tokens for time t
        if item != nothing
            push!(x, copy(item[1])); push!(ygold, copy(item[2])); push!(mask, copy(item[3]))
        else                                                                    # or an item can be nothing marking sentence end
            reset!(m)
            yforw = tagger_forw(m.forw, x; o...)                                # t3:1309
            yback = reverse(tagger_forw(m.back, reverse(x); o...))              # t3:1354
            ypred = tagger_forw(m.pred, yforw, yback; o...)                     # t3:323 pred network should accept two inputs
            tagger_loss(ypred, ygold, mask, loss; o...)
            tagger_bptt(m, ygold, mask, loss; o...)                             # t3:4858
            gcheck && break     # only do one sentence and no update when gradient checking
            tagger_update(m; o...)
            empty!(x); empty!(ygold); empty!(mask)
        end
    end
end

function tagger_forw(net::Net, inputs...; o...)
    N = length(inputs[1])
    ystack = cell(N)
    for n=1:N
        ypred = forw(net, map(x->x[n], inputs)...; seq=true, o...)
        ystack[n] = copy(ypred)
    end
    return ystack
end

function tagger_loss(ypred, ygold, mask, loss; losscnt=nothing, maxnorm=nothing, lossreport=0, o...)
    losscnt==nothing && return
    (yrows, ycols) = size2(ypred[1])
    for i=1:length(ypred)
        @assert (yrows, ycols) == size2(ypred[i])
        ntoks = (mask[i] == nothing ? ycols : sum(mask[i]))
        losscnt[1] += loss(ypred[i], ygold[i]; mask=mask[i], o...)
        losscnt[2] += ntoks/ycols
    end
    if lossreport > 0 && losscnt[2]*ycols > lossreport
        println((exp(losscnt[1]/losscnt[2]), losscnt[1]*ycols, losscnt[2]*ycols, (maxnorm==nothing ? [] : maxnorm)...))
        losscnt[1] = losscnt[2] = 0
        maxnorm == nothing || (maxnorm[1] = maxnorm[2] = 0)
    end
end

function tagger_bptt(m::Tagger, ygold, mask, loss; trn=false, o...)             # t3:4858
    trn || return
    N=length(ygold)
    gforw,gback = cell(N),cell(N)
    for n=N:-1:1
        (gf, gb) = back(m.pred, ygold[n], loss; seq=true, mask=mask[n], getdx=true, o...) # t3:620
        gforw[n],gback[n] = copy(gf),copy(gb)
    end
    for n=1:N
        back(m.back, gback[n]; seq=true, mask=mask[n], o...)                    # t3:2104
    end
    for n=N:-1:1
        back(m.forw, gforw[n]; seq=true, mask=mask[n], o...)                    # t3:2069
    end
end

function tagger_update(m::Tagger; gclip=0, trn=false, maxnorm=nothing, o...)
    trn || return
    g = nothing
    if gclip > 0
        g = gnorm(m)
        gclip=(g > gclip > 0 ? gclip/g : 0)
    end
    update!(m; gclip=gclip, o...)
    if maxnorm != nothing
        w = wnorm(m)
        g == nothing && (g = gnorm(m))
        w > maxnorm[1] && (maxnorm[1]=w)
        g > maxnorm[2] && (maxnorm[2]=g)
    end
end


text:
# CSCI5573Final/daggen.jl
# Licensed under the MIT License. See LICENSE.md file in the project root for
# full license information.

# This is all a bit ugly. Lets just not talk about it.
struct Edge
    ID1::Int64
    ID2::Int64
end

function ListToDictDAG(tasklist, dot_file="")
    dag = Dict{Int64,_Task}()
    comm_cost = Dict{Edge, Int64}()

    for task in tasklist[1:end-1]
        new_task = ParseTask(task)
        dag[new_task.ID] = new_task
    end

    new_task = ParseTask(tasklist[end], true)
    dag[new_task.ID] = new_task

    for (ID, _task) in dag
        for child in _task.Children
            add_dependency!(dag[child], ID)
        end
    end

    # writeDOT("pre", dag)

    # Extract edges
    for (ID, _task) in copy(dag)
        if _task.Type === :TRANSFER

            parent = _task.Dependencies[begin]
            child = _task.Children[begin]

            # Add edges. Just for giggles we do both directions
            comm_cost[Edge(parent, child)] = _task.Cost
            comm_cost[Edge(child, parent)] = _task.Cost

            # Remove edge from the parent child relationship
            filter!(e->e!==ID, dag[parent].Children)
            filter!(e->e!==ID, dag[child].Dependencies)

            # Reconnect the parent and child without edge
            push!(dag[parent].Children, child)
            push!(dag[child].Dependencies, parent)

            # Delete Edge from dag
            delete!(dag, ID)
        end
    end

    # Finally add in all of the non-comm edges
    for (ID, _task) in copy(dag)
        for child in _task.Children
            if !haskey(comm_cost, Edge(ID, child))
                comm_cost[Edge(ID, child)] = 0
                comm_cost[Edge(child, ID)] = 0
            end
        end
    end

    # writeDOT("post", dag)

    if length(dot_file) > 0
        writeDOT(dot_file, dag)
    end

    return dag, comm_cost
end

function daggen_init()
    name = "daggen"
    dir = "daggen"
    touch("daggen_init")

    init_script = """
                #!/bin/bash
                pushd $dir
                make
                mv $(name) ../dag_generator
                rm *.o
                popd"""

    open("daggen_init", "w") do file
        write(file, init_script)
    end

    chmod("daggen_init", 2000)
    run(`./daggen_init`)
    rm("daggen_init")
end

"""
- `outfile`: output file
- `n`: number of tasks
- `mindata`: minimum data size
- `maxdata`: maximum data size
- `minalpha`: minimum Amdahl's law parameter value
- `maxalpha`: maximum Amdahl's law parameter value
- `fat`: dag shape - 1.0=max parallel, 0.0=min parallel
- `density`: number of dependencies
- `regular`: regularity for num tasks per level
- `ccr`: communication(MBytes) to computation(sec) ratio
- `jump`: number of levels spanned by communications
- `dot`: output generated DAG in the DOT format
"""
function daggen(;
    outfile="", num_tasks=-1, mindata=-1, maxdata=-1, minalpha=-1, maxalpha=-1,
    fat=-1, density=-1, regular=-1, ccr=-1, jump=-1, dot=false, verbose=false)

    err = verbose ? "daggen_err.txt" : nothing

    base = `./dag_generator`
    out = (outfile == "") ? `` : ` -o $outfile`
    num = (num_tasks == -1) ? `` : ` -n $(Int64(round(num_tasks/2-sqrt(num_tasks+1))))` # Weird Stuff here
    mind = (mindata == -1) ? `` : ` --mindata $mindata`
    maxd = (maxdata == -1) ? `` : ` --maxdata $maxdata`
    mina = (minalpha == -1) ? `` : ` --minalpha $minalpha`
    maxa = (maxalpha == -1) ? `` : ` --maxalpha $maxalpha`
    fat = (fat == -1) ? `` : ` --fat $fat`
    cr = (ccr == -1) ? `` : ` --ccr $ccr`
    jmp = (jump == -1) ? `` : ` --jump $jump`
    dt = !dot ? `` : ` --dot`

    command = `$base $out $num $mind $maxd $mina $maxa $fat $cr $jmp $dt`
    daggen_init()

    buffer = IOBuffer()
    ps = run(pipeline(`$(command)`, stdout=buffer, stderr=err), wait=true)

    daggen_clean()

    # This processes the IOBuffer as a String array
    tasklist = buffer |> take! |> String |> l->split(l,"\n") |> l->deleteat!(l,[1,2,length(l)])

    num_task = tasklist[1] |> c->split(c, " ") |> n->n[2] |> i->parse(Int, i)

    deleteat!(tasklist, 1) # delete count from tasklist

    return String.(tasklist), num_task
end

function read_daggen(file)
    config = readlines(file)

    # This processes the IOBuffer as a String array
    tasklist = config |> l->deleteat!(l,[1,2])

    num_task = tasklist[1] |> c->split(c, " ") |> n->n[2] |> i->parse(Int, i)

    deleteat!(tasklist, 1) # delete count from tasklist

    return String.(tasklist), num_task
end

function daggen_clean()
    return run(`rm dag_generator`)
end

function writeDOT(Filename, dag)
    open(Filename, "w") do file
        write(file, "digraph DAG_Schedule {\n")

        for (_, task) in sort(dag) # Write self
            color = task.Type == :TRANSFER ? "grey" : "black"

            node = "  T$(task.ID) [size=\"$(task.Cost)\", "
            node *= "overhead=\"$(task.Complexity)\", "
            node *= "color=\"$(color)\"]\n"
            write(file, node)

            for ID in task.Children
                child = dag[ID]
                edge = "  T$(task.ID) -> T$(child.ID) "
                edge *= "[size=\"$(task.Cost)\", "
                edge *= "color=\"$(color)\"]\n"
                write(file, edge)
            end
        end

        write(file, "}\n")
    end
end


text:
#==============================================================================

          Solving the HJB with labor supply using the implicit method

               Translated from matlab code on Ben Moll's website:
                             https://benjaminmoll.com/codes/

==============================================================================#

using LinearAlgebra, SparseArrays, Plots, Roots, LaTeXStrings

γ = 2 # parameter from CRRA utility
φ = .5 # frisch elasticity
ρ = 0.05 # the discount rate
r = 0.03
w = 1.0

z1=.1
z2 =.2
z= [z1 z2]

λ1 = 1.5
λ2 = 1.0

λ = [λ1 λ2]

H = 500 # number of points in the grid space
amin= -0.15
amax= 3.0

a = LinRange(amin,amax,H)
da = (amax-amin)/(H-1)

aa = [a a]
zz = ones(H,1)*z

maxit = 20 # the maximum number of iterations we allow the finite differencing algorithm
crit = 10^(-6) # our critical value
Δ = 1000

dVf=zeros(H,2)
dVb = zeros(H,2)
c= zeros(H,2)

# Create the matrix that captures the systems dynamic from λ process
Aswitch = [-1*sparse(I,H,H)*λ[1] sparse(I,H,H)*λ[1]; sparse(I,H,H)*λ[2] -1*sparse(I,H,H)*λ[2]]

x0 = (w*z1)^(φ*(1-γ)/(1+γ*φ)) #inital value for labor

# create function for labor supply
labor_solve(l, a, z, w, r, γ, φ) = (l .- ((w*z*l .+ r*a).^(-γ*φ))*(w*z)^φ)

# set empty matrix and find all possible values for l
l0 = zeros(500,2)

for i in 1:H
    l_solve1(l) = labor_solve(l, a[i], z1, w, r, γ, φ)
    l01 = find_zero(l_solve1,x0)

    l_solve2(l) = labor_solve(l, a[i], z2, w, r, γ, φ)
    l02 = find_zero(l_solve2,x0)

    l0[i,:] = [l01 l02]
end

v0=zeros(H,2)
v0[:,1] = (w*z[1].*l0[1,1] .+ r.*a).^(1-γ)/(1-γ)/ρ
v0[:,2] = (w*z[2].*l0[1,2] .+ r.*a).^(1-γ)/(1-γ)/ρ
plot(a,v0[:,2])

lmin = l0[1,:]
lmax = l0[H,:]

global v=v0

# set up spaces to save variables of interest
dist=[]; V_n=[]; g_r=[]; dV_r=[]; V_r=[]; c_r=[]; adot=[]

# Loop that finds the equilibrium interest rate
for n in 1:maxit
    global V=v
    push!(V_n,V)
    # Forward differencing
    dVf[1:H-1,:]=(V[2:H,:]-V[1:H-1,:])/da
    dVf[H,:] = (w*z'.*lmax .+ r.*amax).^(-γ) # impose state constraint at the max, just in case
    # Backward differencing
    dVb[2:H,:]=(V[2:H,:]-V[1:H-1,:])/da
    dVb[1,:]= (w*z'.*lmin .+ r.*amin).^(-γ)

    # Find consumption and savings with the forward difference
    cf = dVf.^(-1/γ)
    lf = (dVf.*w.*zz).^φ
    sf = w*zz.*lf +r.*aa -cf

    # Find consumption and savings with the backward difference
    cb= dVb.^(-1/γ)
    lb = (dVb.*w.*zz).^φ
    sb = w*zz.*lb +r.*aa -cb

    # Find consumption and savings at steady state
    c0 = w*zz.*l0 + r.*aa
    dV0 = c0.^(-γ)

    # Now implement the upwind scheme in order to select the best differencing method
    Ib = (sb .< 0)
    If = (sf .> 0)#.*(1 .-Ib) this term was included in Matlab code, messes up julia code for some reason??? Probably rounding differences?
    I0 = (1 .- If - Ib)

    # State constraint at amin is automatically implemented
    global c = cf.*If + cb.*Ib + c0.*I0
    global l = lf.*If + lb.*Ib + l0.*I0
    u = c.^(1-γ)/(1-γ) - l.^(1+1/φ)/(1+1/φ)

    # Construct matrix for the evolution of the system
    global X = -Ib.*sb/da
    global Y = -If.*sf/da + Ib.*sb/da
    global Z = If.*sf/da

    A1 = spdiagm(0=>Y[:,1], -1 => X[2:H,1], 1=> Z[1:H-1,1])
    A2 = spdiagm(0=>Y[:,2], -1 => X[2:H,2], 1=> Z[1:H-1,2])
    AA = [A1 spzeros(H,H); spzeros(H,H) A2]

    global A = AA + Aswitch

    B = (ρ +1/Δ)*sparse(I,2*H,2*H) - A

    u_stacked = [u[:,1]; u[:,2]]
    V_stacked = [V[:,1]; V[:,2]]

    b = u_stacked + V_stacked/Δ
    V_stacked = B\b # Solves the system of equations

    V = [V_stacked[1:H] V_stacked[H+1:2*H]]

    V_change = V - v

    global v = V

    push!(dist, findmax(abs.(V_change))[1])

    if dist[n] < crit
        println("Value Function converged, Iteration=")
        println(n)
        break
    end
end

adot = w.*zz.*l + r.*aa - c

#plot savings
plot(a, adot[:,1], label = L"s_1(a)", xlabel="Wealth, a", ylabel=L"\mathrm{Savings,\hspace{1ex} } s_i(a)")
plot!(a, adot[:,2], label = L"s_2(a)")
plot!(a, zeros(H,1), line=:dash, color=:black, label="")

#plot value function for j=1,2
plot(a, v[:,1], label = L"v_1(a)", xlabel="Wealth, a", ylabel=L"\mathrm{Value Function,\hspace{1ex} } v_i(a)")
plot!(a, v[:,2], label = L"v_2(a)", legend=:bottomright)

#plot labor for j=1,2
plot(a, l[:,1], label = L"l_1(a)", xlabel="Wealth, a", ylabel=L"\mathrm{Labor Supply,\hspace{1ex} } l_i(a)")
plot!(a, l[:,2], label = L"l_2(a)")

#plot consumption for j=1,2
plot(a, c[:,1], label = L"c_1(a)", xlabel="Wealth, a", ylabel=L"\mathrm{Consumption,\hspace{1ex} } c_i(a)")
plot!(a, c[:,2], label = L"c_2(a)", legend=:bottomright)