text:
"""
    stepforward!(prob)

Step forward `prob` one time step.
"""
function stepforward!(prob::Problem)
  stepforward!(prob.sol, prob.clock, prob.timestepper, prob.eqn, prob.vars, prob.params, prob.grid)
end

"""
    stepforward!(prob, nsteps::Int)

Step forward `prob` for `nsteps`.
"""
function stepforward!(prob::Problem, nsteps::Int)
  for step = 1:nsteps
    stepforward!(prob)
  end

  return nothing
end

"""
    stepforward!(prob::Problem, diags, nsteps::Int)

Step forward `prob` for `nsteps`, incrementing `diags` along the way. `diags` may be a
single `Diagnostic` or a `Vector` of `Diagnostic`s.
"""
function stepforward!(prob::Problem, diags, nsteps::Int)
  for step = 1:nsteps
    stepforward!(prob)
    increment!(diags)
  end

  return nothing
end

const fullyexplicitsteppers= [
  :ForwardEuler,
  :RK4,
  :AB3,
  :FilteredForwardEuler,
  :FilteredRK4,
  :FilteredAB3
]

isexplicit(stepper) = any(Symbol(stepper) .== fullyexplicitsteppers)

"""
    TimeStepper(stepper, equation, dt=nothing, dev=CPU(); kw...)

Instantiate the Time`stepper` for `equation` with timestep `dt` and
on the `dev`ice. The `kw` are passed to the timestepper constructor.
"""
function TimeStepper(stepper, equation, dt=nothing, dev::Device=CPU(); kw...)
  fullsteppername = Symbol(stepper, :TimeStepper)

  # Create expression that instantiates the time-stepper, depending on whether
  # timestepper is explicit or not.
  expr = isexplicit(stepper) ? Expr(:call, fullsteppername, equation, dev) :
                               Expr(:call, fullsteppername, equation, dt, dev)

  # Add keyword arguments
  length(kw) > 0 && push!(expr.args, Tuple(Expr(:kw, p.first, p.second) for p in kw)...)

  return eval(expr)
end

# The following time-steppers are implemented below
#
#   * Forward Euler
#   * Filtered Forward Euler
#   * RK4
#   * Filtered RK4
#   * ETDRK4
#   * Filtered ETDRK4
#   * AB3
#   * Filtered AB3
#
# Explicit time-steppers are constructed with the signature
#   ts = ExplicitTimeStepper(eq::Equation)
#
# Implicit time-steppers are constructed with the signature
#   ts = ImplicitTimeStepper(eq::Equation, dt)

# --
# Forward Euler
# --

"""
    ForwardEulerTimeStepper{T} <: AbstractTimeStepper{T}

Type for Forward Euler timestepper.
"""
struct ForwardEulerTimeStepper{T} <: AbstractTimeStepper{T}
  N::T # Explicit linear and nonlinear terms
  ForwardEulerTimeStepper(N::T) where T = new{T}(0N)
end

"""
    ForwardEulerTimeStepper(eq::Equation, dev::Device=CPU())

Construct a Forward Euler timestepper.
"""
ForwardEulerTimeStepper(eq::Equation, dev::Device=CPU()) =
  ForwardEulerTimeStepper(devzeros(dev, eq.T, eq.dims))

function stepforward!(sol, clock, ts::ForwardEulerTimeStepper, eq, vars, params, grid)
  eq.calcN!(ts.N, sol, clock.t, clock, vars, params, grid)
  @. sol += clock.dt * (eq.L * sol + ts.N)
  clock.t += clock.dt
  clock.step += 1

  return nothing
end

"""
    FilteredForwardEulerTimeStepper(eq, dev; filterkwargs...)

Type for Forward Euler timestepper with spectral filtering.
"""
struct FilteredForwardEulerTimeStepper{T,Tf} <: AbstractTimeStepper{T}
       N :: T
  filter :: Tf
end

"""
    FilteredForwardEulerTimeStepper(eq, dev; filterkwargs...)

Construct a forward Euler timestepper with spectral filtering.
"""
function FilteredForwardEulerTimeStepper(equation::Equation, dev::Device=CPU(); filterkwargs...)
  filter = makefilter(equation; filterkwargs...)

  return FilteredForwardEulerTimeStepper(devzeros(dev, equation.T, equation.dims), filter)
end

function stepforward!(sol, clock, ts::FilteredForwardEulerTimeStepper, eq, vars, params, grid)
  eq.calcN!(ts.N, sol, clock.t, clock, vars, params, grid)
  @. sol = ts.filter * (sol + clock.dt * (ts.N + eq.L * sol))
  clock.t += clock.dt
  clock.step += 1

  return nothing
end

# --
# RK4
# --

"""
    RK4TimeStepper{T} <: AbstractTimeStepper{T}

Type for 4th-order Runge-Kutta time stepper for the equation `eq`.
"""
struct RK4TimeStepper{T} <: AbstractTimeStepper{T}
  sol₁ :: T
  RHS₁ :: T
  RHS₂ :: T
  RHS₃ :: T
  RHS₄ :: T
end

"""
    RK4TimeStepper(eq::Equation, dev::Device=CPU())

Construct a 4th-order Runge-Kutta time stepper for the equation `eq`.
"""
function RK4TimeStepper(eq::Equation, dev::Device=CPU())
  @devzeros typeof(dev) eq.T eq.dims N sol₁ RHS₁ RHS₂ RHS₃ RHS₄

  return RK4TimeStepper(sol₁, RHS₁, RHS₂, RHS₃, RHS₄)
end

"""
    FilteredRK4TimeStepper{T,Tf} <: AbstractTimeStepper{T}

Type for 4th-order Runge-Kutta time stepper with spectral filtering for the equation `eq`.
"""
struct FilteredRK4TimeStepper{T,Tf} <: AbstractTimeStepper{T}
    sol₁ :: T
    RHS₁ :: T
    RHS₂ :: T
    RHS₃ :: T
    RHS₄ :: T
  filter :: Tf
end

"""
    FilteredRK4TimeStepper(eq::Equation, dev::Device=CPU(); filterkwargs...)

Construct a 4th-order Runge-Kutta time stepper with spectral filtering for the equation `eq`.
"""
function FilteredRK4TimeStepper(eq::Equation, dev::Device=CPU(); filterkwargs...)
  ts = RK4TimeStepper(eq, dev)
  filter = makefilter(eq; filterkwargs...)

  return FilteredRK4TimeStepper(getfield.(Ref(ts), fieldnames(typeof(ts)))..., filter)
end

function addlinearterm!(RHS, L, sol)
  @. RHS += L*sol

  return nothing
end

function substepsol!(newsol, sol, RHS, dt)
  @. newsol = sol + dt*RHS

  return nothing
end

function RK4substeps!(sol, clock, ts, eq, vars, params, grid, t, dt)
  # Substep 1
  eq.calcN!(ts.RHS₁, sol, t, clock, vars, params, grid)
  addlinearterm!(ts.RHS₁, eq.L, sol)

  # Substep 2
  substepsol!(ts.sol₁, sol, ts.RHS₁, dt/2)
  eq.calcN!(ts.RHS₂, ts.sol₁, t+dt/2, clock, vars, params, grid)
  addlinearterm!(ts.RHS₂, eq.L, ts.sol₁)

  # Substep 3
  substepsol!(ts.sol₁, sol, ts.RHS₂, dt/2)
  eq.calcN!(ts.RHS₃, ts.sol₁, t+dt/2, clock, vars, params, grid)
  addlinearterm!(ts.RHS₃, eq.L, ts.sol₁)

  # Substep 4
  substepsol!(ts.sol₁, sol, ts.RHS₃, dt)
  eq.calcN!(ts.RHS₄, ts.sol₁, t+dt, clock, vars, params, grid)
  addlinearterm!(ts.RHS₄, eq.L, ts.sol₁)

  return nothing
end

function RK4update!(sol, RHS₁, RHS₂, RHS₃, RHS₄, dt)
  @. sol += dt*(RHS₁ / 6 + RHS₂ / 3  + RHS₃ / 3 + RHS₄ / 6)

  return nothing
end

function RK4update!(sol, RHS₁, RHS₂, RHS₃, RHS₄, filter, dt)
  @. sol = filter * (sol + dt*(RHS₁ / 6 + RHS₂ / 3  + RHS₃ / 3 + RHS₄ / 6))

  return nothing
end

function stepforward!(sol, clock, ts::RK4TimeStepper, eq, vars, params, grid)
  RK4substeps!(sol, clock, ts, eq, vars, params, grid, clock.t, clock.dt)
  RK4update!(sol, ts.RHS₁, ts.RHS₂, ts.RHS₃, ts.RHS₄, clock.dt)
  clock.t += clock.dt
  clock.step += 1

  return nothing
end

function stepforward!(sol, clock, ts::FilteredRK4TimeStepper, eq, vars, params, grid)
  RK4substeps!(sol, clock, ts, eq, vars, params, grid, clock.t, clock.dt)
  RK4update!(sol, ts.RHS₁, ts.RHS₂, ts.RHS₃, ts.RHS₄, ts.filter, clock.dt)
  clock.t += clock.dt
  clock.step += 1

  return nothing
end

# ------
# ETDRK4
# ------

"""
    ETDRK4TimeStepper{T,TL} <: AbstractTimeStepper{T}

Type for 4th-order exponential-time-differencing Runge-Kutta time stepper.
"""
struct ETDRK4TimeStepper{T,TL} <: AbstractTimeStepper{T}
  # ETDRK4 coefficents
        ζ :: TL
        α :: TL
        β :: TL
        Γ :: TL
  expLdt  :: TL
  exp½Ldt :: TL
     sol₁ :: T
     sol₂ :: T
       N₁ :: T
       N₂ :: T
       N₃ :: T
       N₄ :: T
end

"""
    ETDRK4TimeStepper(eq::Equation, dt, dev::Device=CPU())

Construct a 4th-order exponential-time-differencing Runge-Kutta time stepper.
"""
function ETDRK4TimeStepper(eq::Equation, dt, dev::Device=CPU())
  dt = fltype(eq.T)(dt) # ensure dt is correct type.
  expLdt, exp½Ldt = getexpLs(dt, eq)
  ζ, α, β, Γ = getetdcoeffs(dt, eq.L)
  @devzeros typeof(dev) eq.T eq.dims sol₁ sol₂ N₁ N₂ N₃ N₄

  return ETDRK4TimeStepper(ζ, α, β, Γ, expLdt, exp½Ldt, sol₁, sol₂, N₁, N₂, N₃, N₄)
end

"""
    FilteredETDRK4TimeStepper{T,TL,Tf} <: AbstractTimeStepper{T}

Type for 4th-order exponential-time-differencing Runge-Kutta time stepper with spectral filtering.
"""
struct FilteredETDRK4TimeStepper{T,TL,Tf} <: AbstractTimeStepper{T}
  # ETDRK4 coefficents:
        ζ :: TL
        α :: TL
        β :: TL
        Γ :: TL
  expLdt  :: TL
  exp½Ldt :: TL
     sol₁ :: T
     sol₂ :: T
       N₁ :: T
       N₂ :: T
       N₃ :: T
       N₄ :: T
   filter :: Tf
end

"""
    FilteredETDRK4TimeStepper(equation, dt; filterkwargs...)

Construct a 4th-order exponential-time-differencing Runge-Kutta time stepper with spectral filtering.
"""
function FilteredETDRK4TimeStepper(equation::Equation, dt, dev::Device=CPU(); filterkwargs...)
  timestepper = ETDRK4TimeStepper(equation, dt, dev)
  filter = makefilter(equation; filterkwargs...)

  return FilteredETDRK4TimeStepper(getfield.(Ref(timestepper), fieldnames(typeof(timestepper)))..., filter)
end

function ETDRK4update!(sol, expLdt, α, β, Γ, N₁, N₂, N₃, N₄)
  @. sol = (expLdt * sol +  α * N₁
                         + 2β * (N₂ + N₃)
                         +  Γ * N₄ )

  return nothing
end

function ETDRK4update!(sol, ts, filter)
  @. sol = filter * (ts.expLdt * sol +     ts.α * ts.N₁
                                     + 2 * ts.β * (ts.N₂ + ts.N₃)
                                     +     ts.Γ * ts.N₄ )

  return nothing
end

function ETDRK4substep12!(sol₁, exp½Ldt, sol, ζ, N)
  @. sol₁ = exp½Ldt * sol + ζ * N

  return nothing
end

function ETDRK4substep3!(sol₂, exp½Ldt, sol₁, ζ, N₁, N₃)
  @. sol₂ = exp½Ldt * sol₁ + ζ * (2N₃ - N₁)

  return nothing
end

function ETDRK4substeps!(sol, clock, ts, eq, vars, params, grid)
  # Substep 1
  eq.calcN!(ts.N₁, sol, clock.t, clock, vars, params, grid)
  ETDRK4substep12!(ts.sol₁, ts.exp½Ldt, sol, ts.ζ, ts.N₁)
  @. ts.sol₁ = ts.exp½Ldt * sol + ts.ζ * ts.N₁

  # Substep 2
  t2 = clock.t + clock.dt/2
  eq.calcN!(ts.N₂, ts.sol₁, t2, clock, vars, params, grid)
  ETDRK4substep12!(ts.sol₂, ts.exp½Ldt, sol, ts.ζ, ts.N₂)

  # Substep 3
  eq.calcN!(ts.N₃, ts.sol₂, t2, clock, vars, params, grid)
  ETDRK4substep3!(ts.sol₂, ts.exp½Ldt, ts.sol₁, ts.ζ, ts.N₁, ts.N₃)

  # Substep 4
  t3 = clock.t + clock.dt
  eq.calcN!(ts.N₄, ts.sol₂, t3, clock, vars, params, grid)

  return nothing
end

function stepforward!(sol, clock, ts::ETDRK4TimeStepper, eq, vars, params, grid)
  ETDRK4substeps!(sol, clock, ts, eq, vars, params, grid)
  ETDRK4update!(sol, ts.expLdt, ts.α, ts.β, ts.Γ, ts.N₁, ts.N₂, ts.N₃, ts.N₄)
  clock.t += clock.dt
  clock.step += 1

  return nothing
end

function stepforward!(sol, clock, ts::FilteredETDRK4TimeStepper, eq, vars, params, grid)
  ETDRK4substeps!(sol, clock, ts, eq, vars, params, grid)
  ETDRK4update!(sol, ts, ts.filter) # update
  clock.t += clock.dt
  clock.step += 1

  return nothing
end

# ---
# AB3
# ---

const ab3h1 = 23/12
const ab3h2 = 16/12
const ab3h3 = 5/12

"""
    AB3TimeStepper{T} <: AbstractTimeStepper{T}

Type for 3rd order Adams-Bashforth time-stepper.
"""
struct AB3TimeStepper{T} <: AbstractTimeStepper{T}
  RHS::T
  RHS₋₁::T
  RHS₋₂::T
end

"""
    FilteredAB3TimeStepper(equation, dev; filterkwargs...)

Construct a 3rd order Adams-Bashforth time stepper.
"""
function AB3TimeStepper(equation::Equation, dev::Device=CPU())
  @devzeros typeof(dev) equation.T equation.dims RHS RHS₋₁ RHS₋₂

  return AB3TimeStepper(RHS, RHS₋₁, RHS₋₂)
end

"""
    FilteredAB3TimeStepper{T} <: AbstractTimeStepper{T}

Type for 3rd order Adams-Bashforth time-stepper with spectral filtering.
"""
struct FilteredAB3TimeStepper{T, Tf} <: AbstractTimeStepper{T}
   RHS   :: T
   RHS₋₁ :: T
   RHS₋₂ :: T
  filter :: Tf
end

"""
    FilteredAB3TimeStepper(equation, dev; filterkwargs...)

Construct a 3rd order Adams-Bashforth time stepper with spectral filtering.
"""
function FilteredAB3TimeStepper(equation::Equation, dev::Device=CPU(); filterkwargs...)
  timestepper = AB3TimeStepper(equation, dev)
  filter = makefilter(equation; filterkwargs...)

  return FilteredAB3TimeStepper(getfield.(Ref(timestepper), fieldnames(typeof(timestepper)))..., filter)
end

function AB3update!(sol, ts, clock)
  if clock.step < 3  # forward Euler steps to initialize AB3
    @. sol += clock.dt * ts.RHS    # Update
  else   # Otherwise, stepforward with 3rd order Adams Bashforth:
    @. sol += clock.dt * (ab3h1 * ts.RHS - ab3h2 * ts.RHS₋₁ + ab3h3 * ts.RHS₋₂)
  end

  return nothing
end

function AB3update!(sol, ts::FilteredAB3TimeStepper, clock)
  if clock.step < 3  # forward Euler steps to initialize AB3
    @. sol = ts.filter * (sol + clock.dt * ts.RHS)    # Update
  else   # Otherwise, stepforward with 3rd order Adams Bashforth:
    @. sol = ts.filter * (sol + clock.dt * (ab3h1 * ts.RHS - ab3h2 * ts.RHS₋₁ + ab3h3 * ts.RHS₋₂))
  end

  return nothing
end

function stepforward!(sol, clock, ts::AB3TimeStepper, eq, vars, params, grid)
  eq.calcN!(ts.RHS, sol, clock.t, clock, vars, params, grid)
  addlinearterm!(ts.RHS, eq.L, sol)

  AB3update!(sol, ts, clock)
  clock.t += clock.dt
  clock.step += 1

  @. ts.RHS₋₂ = ts.RHS₋₁          # Store
  @. ts.RHS₋₁ = ts.RHS            # ... previous values of RHS

  return nothing
end

function stepforward!(sol, clock, ts::FilteredAB3TimeStepper, eq, vars, params, grid)
  eq.calcN!(ts.RHS, sol, clock.t, clock, vars, params, grid)
  addlinearterm!(ts.RHS, eq.L, sol)

  AB3update!(sol, ts, clock)
  clock.t += clock.dt
  clock.step += 1

  @. ts.RHS₋₂ = ts.RHS₋₁          # Store
  @. ts.RHS₋₁ = ts.RHS            # ... previous values of RHS

  return nothing
end

# --
# Timestepper utils
# --

function getexpLs(dt, equation)
  expLdt  = @. exp(dt * equation.L)
  exp½Ldt = @. exp(dt * equation.L/2)

  return expLdt, exp½Ldt
end

"""
    getetdcoeffs(dt, L; ncirc=32, rcirc=1)

Calculate ETDRK4 coefficients associated with the (diagonal) linear coefficient
L by integrating over a small circle in complex space.
"""
function getetdcoeffs(dt, L; ncirc=32, rcirc=1)
  shape = Tuple(cat(ncirc, ones(Int, ndims(L)), dims=1))

  circ = zeros(Complex{Float64}, shape) # use double precision for this calculation
  circ .= rcirc * exp.(2π * im/ncirc * (0.5:1:(ncirc-0.5)))
  circ = permutedims(circ, ndims(circ):-1:1)

  zc = dt * L .+ circ
  M = ndims(L) + 1

  # Four coefficients: ζ, α, β, Γ
  ζc = @.               ( exp(zc/2)-1 ) / zc
  αc = @. ( -4 - zc + exp(zc) * (4 - 3zc + zc^2) ) / zc^3
  βc = @.    ( 2  + zc + exp(zc) * (-2 + zc) ) / zc^3
  Γc = @. ( -4 - 3zc - zc^2 + exp(zc) * (4 - zc) ) / zc^3

  ζ = dt * dropdims(mean(ζc, dims=M), dims=M)
  α = dt * dropdims(mean(αc, dims=M), dims=M)
  β = dt * dropdims(mean(βc, dims=M), dims=M)
  Γ = dt * dropdims(mean(Γc, dims=M), dims=M)

  if eltype(L) <: Real # this is conservative, but unclear if necessary
    ζ = real.(ζ)
    α = real.(α)
    β = real.(β)
    Γ = real.(Γ)
  end

  return ζ, α, β, Γ
end

"""
    step_until!(prob, stop_time)

Step forward `prob` until `stop_time`. Cannot be used with ETDRK4 time steppers.

See also: [`stepforward!`](@ref)
"""
step_until!(prob, stop_time) = step_until!(prob, prob.timestepper, stop_time)

step_until!(prob, ::Union{ETDRK4TimeStepper, FilteredETDRK4TimeStepper}, stop_time) =
  error("step_until! requires fully explicit time stepper; does not work with ETDRK4")

function step_until!(prob, timestepper, stop_time)
  # Throw an error if stop_time is not greater than the current problem time
  stop_time > prob.clock.t || error("stop_time must be greater than prob.clock.t")

  # Extract current time step
  dt = prob.clock.dt

  # Step forward until just before stop_time
  time_interval = stop_time - prob.clock.t
  nsteps = floor(Int, time_interval / dt)
  stepforward!(prob, nsteps)

  # Take one final small step so that prob.clock.t = stop_time
  t_remaining = time_interval - prob.clock.t
  prob.clock.dt = t_remaining
  stepforward!(prob)

  # Restore previous time-step
  prob.clock.dt = dt

  return nothing
end


text:
function check_model_estimated(len::Int)
    if len == 0
        println("Score Driven Model does not have unknowns.")
        return true
    end
    return false
end

function NaN2zero!(m::Matrix{T}, i::Integer) where T
    for j in axes(m, 2)
        if isnan(m[i, j])
            m[i, j] = zero(T)
        end
    end
    return
end

function big_threshold!(m::Matrix{T}, threshold::F, i::Integer) where {T, F}
    for j in axes(m, 2)
        if m[i, j] >= threshold
            m[i, j] = threshold
        end
        if m[i, j] <= -threshold
            m[i, j] = -threshold
        end
    end
    return
end

function small_threshold!(m::Matrix{T}, threshold::F, i::Integer) where {T, F}
    for j in axes(m, 2)
        if m[i, j] <= threshold && m[i, j] >= 0
            m[i, j] = threshold
        end
        if m[i, j] >= -threshold && m[i, j] <= 0
            m[i, j] = -threshold
        end
    end
    return
end

function sample_observation(dist::Distribution)
    return rand(dist)
end

function find_unknowns(v::Vector{T}) where T
    return findall(isnan, v)
end

function find_unknowns(m::Matrix{T}) where T
    return findall(isnan, vec(m))
end

ext: jl


text:
using CellListMapBenchmarks
using Test

@testset "CellListMapBenchmarks.jl" begin
    # Write your tests here.
end


text:
export CG_generic

"""
A non-linear conjugate gradient algorithm.
Can use differents conjugate gradient formulae:
  - Fletcher & Reeves
  - Hager & Zhang
  - Pollack & Ribière
  - Hestenes & Stiefel

  Inputs:
  - An AbstractNLPModel (our problem)
  - An NLPStopping (stopping fonctions/criterion of our problem)
  - (opt) verbose = true will show function value, norm of gradient and the
    inner product of the gradient and the descent direction
  - (opt) hessian_rep, different way to implement the Hessian. The default
    way is a LinearOperator

  Outputs:
  - An NLPAtX, so all the information of the last iterate
  - A boolean true if we reahced an optimal solution, false otherwise
"""
function CG_generic(nlp        :: AbstractNLPModel,
                    nlp_stop   :: NLPStopping;
                    verbose    :: Bool = false,
                    linesearch :: Function = armijo_ls,
                    CG_formula :: Function = formula_HZ,
                    scaling    :: Bool = true,
                    kwargs...)
    nlp_at_x = nlp_stop.current_state
    x = copy(nlp.meta.x0)
    n = nlp.meta.nvar

    xt = Array{Float64}(undef, n)
    ∇ft = Array{Float64}(undef, n)

    f = obj(nlp, x)
    ∇f = grad(nlp, x)

    iter = 0

    OK = update_and_start!(nlp_stop, x = x, fx = f, gx = ∇f, g0 = ∇f)
    ∇fNorm = norm(nlp_at_x.gx)

    verbose && @printf("%4s  %8s  %7s  %8s \n", " iter", "f", "‖∇f‖", "α")
    verbose && @printf("%5d  %9.2e  %8.1e", iter, nlp_at_x.fx, ∇fNorm)

    β = 0.0
    # d = zeros(∇f)
    d = zero(nlp_at_x.gx)
    scale = 1.0

    h = LineModel(nlp, x, d * scale)

    while !OK
        d = - ∇f + β * d
        slope = ∇f⋅d
        if slope > 0.0  # restart with negative gradient
          #stalled_ascent_dir = true
          d = - ∇f
          slope =  ∇f⋅d
        end

        h = redirect!(h, x, d * scale)
        ls_at_t = LSAtT(0.0, h₀ = nlp_at_x.fx, g₀ = slope)
        stop_ls = LS_Stopping(h, (x, y) -> armijo(x, y; kwargs...), ls_at_t)
        verbose && println(" ")
        ls_at_t, good_step_size = linesearch(h, stop_ls)

        ls_at_t.x *= scale
        xt = nlp_at_x.x + ls_at_t.x * d
        ft = obj(nlp, xt);
        ∇ft = grad(nlp, xt);

        ∇ft = grad!(nlp, xt, ∇ft)
        # Move on.
        s = xt - x
        y = ∇ft - ∇f
        β = 0.0
        if (∇ft⋅∇f) < 0.2 * (∇ft⋅∇ft)   # Powell restart
            β = CG_formula(∇f, ∇ft, s, d)
        end
        if scaling
            scale = (y⋅s) / (y⋅y)
        end
        if scale <= 0.0
            #println(" scale = ",scale)
            #println(" ∇f⋅s = ",∇f⋅s,  " ∇ft⋅s = ",∇ft⋅s)
            scale = 1.0
        end
        x = xt
        f = ft
        dxt = xt - x
        BLAS.blascopy!(n, ∇ft, 1, ∇f, 1)

        ∇fNorm = norm(∇ft, Inf)
        iter = iter + 1
        verbose && @printf("%4d  %9.2e  %8.1e %8.2e", iter, nlp_at_x.fx, ∇fNorm, ls_at_t.x)

        OK = update_and_stop!(nlp_stop, x = xt, fx = ft, gx = ∇ft)
    end
    verbose && @printf("\n")

    return nlp_stop, nlp_stop.meta.optimal
end


text:
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.11.2
#   kernelspec:
#     display_name: Julia 1.8.0-DEV
#     language: julia
#     name: julia-1.8
# ---

# %%
struct Foo
    var"local"::Float64
end
@show foo = Foo(1.2)
@show foo.local;

# %%
:(struct Foo
    var"local"::Float64
end) |> Base.remove_linenums! |> print

# %%

max_stars_repo_name: genkuroki/public

max_issues_repo_name: genkuroki/public

max_forks_repo_name: genkuroki/public


text:
#
## String methods
#

function +(a::Sidef_Types_String_String, b::Sidef_Types_String_String)
    Sidef_Types_String_String(a.value * b.value)
end

function *(a::Sidef_Types_String_String, b::Sidef_Types_Number_Number)
    Sidef_Types_String_String(repeat(a.value, Int(b.value)))
end

function eq(a::Sidef_Types_String_String, b::Sidef_Types_String_String)
    (a.value == b.value) ? TRUE : FALSE
end

function eq(::Sidef_Types_String_String, ::Sidef_Object)
    FALSE
end

function ne(a::Sidef_Types_String_String, b::Sidef_Types_String_String)
    (a.value != b.value) ? TRUE : FALSE
end

function ne(::Sidef_Types_String_String, ::Sidef_Object)
    TRUE
end

function string(s::Sidef_Types_String_String)
    s.value
end

max_stars_repo_name: trizen/Sidef-Deparse-Julia

max_issues_repo_name: trizen/Sidef-Deparse-Julia

max_forks_repo_name: trizen/Sidef-Deparse-Julia


text:
using LinearAlgebra
using PolynomialBasis
using ImplicitDomainQuadrature
using Revise
using CutCell
include("../../../test/useful_routines.jl")

function bulk_modulus(l, m)
    return l + 2m / 3
end

function lame_lambda(k, m)
    return k - 2m / 3
end

function solve_for_displacement(
    basis,
    cellquads,
    interfacequads,
    facequads,
    cutmesh,
    levelset,
    levelsetcoeffs,
    stiffness,
    theta0,
    penalty,
)

    bilinearforms = CutCell.BilinearForms(basis, cellquads, stiffness, cutmesh)
    interfacecondition = CutCell.incoherent_interface_condition(
        basis,
        interfacequads,
        stiffness,
        cutmesh,
        penalty,
    )

    sysmatrix = CutCell.SystemMatrix()
    sysrhs = CutCell.SystemRHS()

    CutCell.assemble_bilinear_form!(sysmatrix, bilinearforms, cutmesh)
    CutCell.assemble_interface_condition!(sysmatrix, interfacecondition, cutmesh)
    CutCell.assemble_bulk_transformation_linear_form!(
        sysrhs,
        transfstress,
        basis,
        cellquads,
        cutmesh,
    )
    CutCell.assemble_incoherent_interface_transformation_rhs!(
        sysrhs,
        transfstress,
        basis,
        interfacequads,
        cutmesh,
    )

    matrix = CutCell.make_sparse(sysmatrix, cutmesh)
    rhs = CutCell.rhs(sysrhs, cutmesh)

    topleftnodeid = CutCell.nodes_per_mesh_side(cutmesh.mesh)[2]
    CutCell.apply_dirichlet_bc!(matrix, rhs, [1, topleftnodeid], 1, 0.0, 2)
    CutCell.apply_dirichlet_bc!(matrix, rhs, [1], 2, 0.0, 2)

    nodaldisplacement = matrix \ rhs

    return nodaldisplacement
end

function spatial_closest_points(refclosestpoints, refclosestcellids, mesh)
    dim, npts = size(refclosestpoints)
    spclosestpoints = zeros(dim, npts)
    for i = 1:npts
        cellmap = CutCell.cell_map(mesh, refclosestcellids[i])
        spclosestpoints[:, i] .= cellmap(refclosestpoints[:, i])
    end
    return spclosestpoints
end

function normal_stress_component(stressvector, normal)
    snn =
        normal[1] * stressvector[1] * normal[1] +
        2.0 * normal[1] * stressvector[3] * normal[2] +
        normal[2] * stressvector[2] * normal[2]
end

function normal_stress_component_over_points(stresses,normals)
    numstresscomponents,npts = size(stresses)
    @assert size(normals) == (2,npts)

    normalstresscomp = zeros(npts)
    for i = 1:npts
        normalstresscomp[i] = normal_stress_component(stresses[:,i],normals[:,i])
    end
    return normalstresscomp
end

function stress_inner_product_over_points(stresses)
    return stresses[1,:].^2 + stresses[2,:].^2 + 2.0*stresses[3,:].^2 + stresses[4,:].^2
end

function angular_position(points)
    cpoints = points[1, :] + im * points[2, :]
    return rad2deg.(angle.(cpoints))
end

K1, K2 = 247.0e9, 192.0e9    # Pa
mu1, mu2 = 126.0e9, 87.0e9   # Pa
lambda1 = lame_lambda(K1, mu1)
lambda2 = lame_lambda(K2, mu2)
stiffness = CutCell.HookeStiffness(lambda1, mu1, lambda2, mu2)

rho1 = 3.93e3   # Kg/m^3
rho2 = 3.68e3   # Kg/m^3

V1 = inv(rho1)
V2 = inv(rho2)
theta0 = -0.067
ΔG = -1022.0    # J/mol
molarmass = 0.147   # Kg/mol
diffG0 = ΔG/molarmass # J/Kg

width = 1.0e-3
penaltyfactor = 1e2

polyorder = 3
numqp = required_quadrature_order(polyorder) + 2
nelmts = 33
center = [width / 2, width / 2]
inradius = width / 4
outradius = width

transfstrain = CutCell.plane_transformation_strain(theta0)
transfstress = CutCell.plane_strain_transformation_stress(lambda1, mu1, theta0)

dx = width / nelmts
meanmoduli = 0.5 * (lambda1 + lambda2 + mu1 + mu2)
penalty = penaltyfactor / dx * meanmoduli

basis = TensorProductBasis(2, polyorder)
mesh = CutCell.Mesh([0.0, 0.0], [width, width], [nelmts, nelmts], basis)
levelset = InterpolatingPolynomial(1, basis)
levelsetcoeffs =
    CutCell.levelset_coefficients(x -> -circle_distance_function(x, center, inradius), mesh)

cutmesh = CutCell.CutMesh(levelset, levelsetcoeffs, mesh)
cellquads = CutCell.CellQuadratures(levelset, levelsetcoeffs, cutmesh, numqp)
interfacequads = CutCell.InterfaceQuadratures(levelset, levelsetcoeffs, cutmesh, numqp)
facequads = CutCell.FaceQuadratures(levelset, levelsetcoeffs, cutmesh, numqp)

nodaldisplacement = solve_for_displacement(
    basis,
    cellquads,
    interfacequads,
    facequads,
    cutmesh,
    levelset,
    levelsetcoeffs,
    stiffness,
    theta0,
    penalty,
)

refseedpoints, spatialseedpoints, seedcellids =
    CutCell.seed_zero_levelset_with_interfacequads(interfacequads, cutmesh)
nodalcoordinates = CutCell.nodal_coordinates(cutmesh)

tol = 1e-8
boundingradius = 3.0
refclosestpoints, refclosestcellids, refgradients =
    CutCell.closest_reference_points_on_zero_levelset(
        nodalcoordinates,
        refseedpoints,
        spatialseedpoints,
        seedcellids,
        levelset,
        levelsetcoeffs,
        mesh,
        tol,
        boundingradius,
    )
spclosestpoints = spatial_closest_points(refclosestpoints,refclosestcellids,cutmesh)

relspatialpoints = spclosestpoints .- center
angularposition = angular_position(relspatialpoints)
sortidx = sortperm(angularposition)
angularposition = angularposition[sortidx]

refclosestpoints = refclosestpoints[:,sortidx]
refclosestcellids = refclosestcellids[sortidx]
refgradients = refgradients[:,sortidx]

invjac = CutCell.inverse_jacobian(cutmesh)
normals = diagm(invjac) * refgradients
CutCell.normalize_normals!(normals)

product_stress_at_cp = CutCell.product_stress_at_reference_points(
    refclosestpoints,
    refclosestcellids,
    basis,
    stiffness,
    transfstress,
    theta0,
    nodaldisplacement,
    cutmesh,
)
parent_stress_at_cp = CutCell.parent_stress_at_reference_points(
    refclosestpoints,
    refclosestcellids,
    basis,
    stiffness,
    nodaldisplacement,
    cutmesh,
)

productpressure = CutCell.pressure_at_points(product_stress_at_cp)
parentpressure = CutCell.pressure_at_points(parent_stress_at_cp)

productdevstress = CutCell.deviatoric_stress_at_points(product_stress_at_cp,productpressure)
parentdevstress = CutCell.deviatoric_stress_at_points(parent_stress_at_cp,parentpressure)

productdevnorm = stress_inner_product_over_points(productdevstress)
parentdevnorm = stress_inner_product_over_points(parentdevstress)

productnormaldevstress = normal_stress_component_over_points(productdevstress,normals)
parentnormaldevstress = normal_stress_component_over_points(parentdevstress,normals)

productspecificvolume = V1*(1.0 .- productpressure/K1)
parentspecificvolume = V2*(1.0 .- parentpressure/K2)

productp1 = productspecificvolume .* productpressure
productp2 = V1/(2K1)*(productpressure.^2)
productp3 = -(productspecificvolume .* productnormaldevstress)
productp4 = V1/(4mu1)*productdevnorm
productpotential = productp1+productp2+productp3+productp4

parentp1 = parentspecificvolume .* parentpressure
parentp2 = V2/(2K2)*(parentpressure.^2)
parentp3 = -(parentspecificvolume .* parentnormaldevstress)
parentp4 = V2/(4mu2)*parentdevnorm
parentpotential = parentp1+parentp2+parentp3+parentp4

potentialdifference = diffG0 .+ productpotential - parentpotential

fig,ax = PyPlot.subplots()
ax.plot(angularposition,potentialdifference)
ax.grid()
ax.set_xlabel("Angular position (deg)")
ax.set_ylabel("Potential difference (J/Kg)")
ax.set_title("Potential difference along interface circumference")
fig.tight_layout()
fig.savefig("potential-difference.png")

fig,ax = PyPlot.subplots()
ax.plot(angularposition,parentp1,label="term1")
ax.plot(angularposition,parentp2,label="term2")
ax.plot(angularposition,parentp3,label="term3")
ax.plot(angularposition,parentp4,label="term4")
ax.plot(angularposition,parentpotential,label="potential")
ax.set_xlabel("Angular position (deg)")
ax.set_ylabel("Energy density (J/Kg)")
ax.set_title("Potential on parent side of interface")
ax.set_ylim(-1e6,2e6)
ax.grid()
ax.legend()
fig.tight_layout()
fig.savefig("parent-potential.png")

fig,ax = PyPlot.subplots()
ax.plot(angularposition,productp1,label="term1")
ax.plot(angularposition,productp2,label="term2")
ax.plot(angularposition,productp3,label="term3")
ax.plot(angularposition,productp4,label="term4")
ax.plot(angularposition,productpotential,label="potential")
ax.set_xlabel("Angular position (deg)")
ax.set_ylabel("Energy density (J/Kg)")
ax.set_ylim(-1e6,2e6)
ax.set_title("Potential on product side of interface")
ax.grid()
ax.legend()
fig.tight_layout()
fig.savefig("product-potential.png")


text:
"""
    abstract type AbstractIntervalGrid{T} <: AbstractGrid1d{T}

An AbstractIntervalGrid is a grid that is defined on an interval, i.e. it is connected.
"""
abstract type AbstractIntervalGrid{T} <: AbstractGrid1d{T}
end

isperiodic(::AbstractIntervalGrid) = false
support(grid::AbstractIntervalGrid) = Interval(grid[1], grid[end])
size(grid::AbstractIntervalGrid) = (grid.n,)

instantiate(::Type{T}, n::Int, ::Type{ELT})  where {T<:AbstractIntervalGrid,ELT} = T(n,UnitInterval{ELT}())

"""
    abstract type AbstractEquispacedGrid{T} <: AbstractIntervalGrid{T}

An equispaced grid has equispaced points, and therefore it has a step.
"""
abstract type AbstractEquispacedGrid{T} <: AbstractIntervalGrid{T}
end

range(grid::AbstractEquispacedGrid) = grid.range
==(g1::AbstractEquispacedGrid, g2::AbstractEquispacedGrid) =
    range(g1)==range(g2) && support(g1) == support(g2)

size(grid::AbstractEquispacedGrid) = (length(range(grid)),)
step(grid::AbstractEquispacedGrid) = step(range(grid))
unsafe_getindex(grid::AbstractEquispacedGrid, i::Int) = unsafe_getindex(range(grid), i)

"""
    struct EquispacedGrid{T} <: AbstractEquispacedGrid{T}

An equispaced grid with n points on an interval [a,b], including the endpoints.
It has step (b-a)/(n-1).

# Example
```jldocs
julia> EquispacedGrid(4,0,1)
4-element EquispacedGrid{Float64}:
 0.0
 0.3333333333333333
 0.6666666666666666
 1.0
```
"""
struct EquispacedGrid{T} <: AbstractEquispacedGrid{T}
    # Use StepRangeLen for higher precision
    range   :: LinRange{T}

    EquispacedGrid{T}(n::Int, a, b) where {T} = new(LinRange(T(a),T(b),n))
end

name(g::EquispacedGrid) = "Equispaced grid"

"""
    struct PeriodicEquispacedGrid{T} <: AbstractEquispacedGrid{T}

A periodic equispaced grid is an equispaced grid that omits the right endpoint.
It has step (b-a)/n.

# Example
```jldocs
julia> PeriodicEquispacedGrid(4,0,1)
4-element PeriodicEquispacedGrid{Float64}:
 0.0
 0.25
 0.5
 0.75
```
"""
struct PeriodicEquispacedGrid{T} <: AbstractEquispacedGrid{T}
    range   :: LinRange{T}
    a   ::  T
    b   ::  T

    PeriodicEquispacedGrid{T}(n::Int, a, b) where {T} = new(LinRange(T(a),T(b),n+1)[1:end-1], a, b)
end

name(::PeriodicEquispacedGrid) = "Periodic equispaced grid"
support(grid::PeriodicEquispacedGrid) = Interval(grid.a, grid.b)
isperiodic(::PeriodicEquispacedGrid) = true

"""
    struct MidpointEquispacedGrid{T} <: AbstractEquispacedGrid{T}

A MidpointEquispaced grid is an equispaced grid with grid points in the centers of the equispaced
subintervals. In other words, this is a DCT-II grid.
It has step `(b-a)/n`.

# Example
```jldocs
julia> MidpointEquispacedGrid(4,0,1)
4-element MidpointEquispacedGrid{Float64}:
 0.125
 0.375
 0.6249999999999999
 0.875
```
"""
struct MidpointEquispacedGrid{T} <: AbstractEquispacedGrid{T}
    range   ::LinRange{T}
    a   ::  T
    b   ::  T

    MidpointEquispacedGrid{T}(n::Int, a, b) where {T} = new(LinRange(T(a),T(b),2n+1)[2:2:end], a, b)
end

name(g::MidpointEquispacedGrid) = "Equispaced midpoints grid"
support(grid::MidpointEquispacedGrid) = Interval(grid.a, grid.b)
isperiodic(::MidpointEquispacedGrid) = true

"""
    struct FourierGrid{T} <: AbstractEquispacedGrid{T}

A Fourier grid is a periodic equispaced grid on the interval [0,1).

# example
```jldocs
julia> FourierGrid(4)
4-element FourierGrid{Float64}:
 0.0
 0.25
 0.5
 0.75
```
"""
struct FourierGrid{T} <: AbstractEquispacedGrid{T}
    range ::LinRange{T}

    FourierGrid{T}(n::Int) where {T} = new(LinRange(T(0),T(1),n+1)[1:end-1])
end

name(g::FourierGrid) = "Periodic Fourier grid"
support(g::FourierGrid{T}) where {T} = UnitInterval{T}()
isperiodic(::FourierGrid) = true

include("gauss.jl")

support(::ChebyshevTNodes{T}) where T = ChebyshevInterval{T}()
name(g::ChebyshevTNodes) = "ChebyshevT nodes"

name(g::ChebyshevExtremae) = "Chebyshev extremae"
support(::ChebyshevExtremae{T}) where T = ChebyshevInterval{T}()

name(g::ChebyshevUNodes) = "ChebyshevU nodes"
support(::ChebyshevUNodes{T}) where T = ChebyshevInterval{T}()

name(g::LegendreNodes) = "Legendre nodes"
support(::LegendreNodes{T}) where T = ChebyshevInterval{T}()
LegendreNodes{T}(n::Int) where T = gausslegendre(T, n)[1]

name(g::LaguerreNodes) = "Laguerre nodes α=$(g.α)"
support(::LaguerreNodes{T}) where T = HalfLine{T}()
similargrid(grid::LaguerreNodes, T, n::Int) = LaguerreNodes{T}(n, T(grid.α))
LaguerreNodes(n::Int, α::T) where T = gausslaguerre(T, n, α)[1]
LaguerreNodes{T}(n::Int, α::T) where T = gausslaguerre(T, n, α)[1]

name(g::HermiteNodes) = "Hermite nodes"
support(::HermiteNodes{T}) where T = DomainSets.GeometricSpace{T}()
HermiteNodes{T}(n::Int) where T = gausshermite(T, n)[1]

name(g::JacobiNodes) = "Jacobi nodes α=$(g.α), β=$(g.β)"
support(::JacobiNodes{T}) where T = ChebyshevInterval{T}()
similargrid(grid::JacobiNodes, T, n::Int) = JacobiNodes{T}(n, T(grid.α), T(grid.β))
JacobiNodes{T}(n::Int, α::T, β::T) where T = gaussjacobi(T, n, α, β)[1]
JacobiNodes(n::Int, α::T, β::T) where T = gaussjacobi(T, n, α, β)[1]

# Grids with flexible support
for GRID in (:PeriodicEquispacedGrid, :MidpointEquispacedGrid, :EquispacedGrid)
    @eval $GRID(n::Int, d::AbstractInterval) =
        $GRID(n, endpoints(d)...)
    @eval similargrid(grid::$GRID, ::Type{T}, n::Int) where {T} =
        $GRID{T}(n, map(T, endpoints(support(grid)))...)
    @eval rescale(grid::$GRID, a, b) =
        $GRID{promote_type(typeof(a/2),typeof(b/2),eltype(grid))}(length(grid), a, b)
    @eval $GRID(n::Int, a, b) =
        $GRID{promote_type(typeof(a/2),typeof(b/2))}(n, a, b)
    @eval mapped_grid(grid::$GRID, map::AffineMap) =
        $GRID(length(grid), endpoints(map*support(grid))...)
end

# Grids with fixed support and one variable
for GRID in (:FourierGrid, :ChebyshevNodes, :ChebyshevExtremae, :ChebyshevUNodes, :LegendreNodes, :HermiteNodes)
    @eval similargrid(g::$GRID, ::Type{T}, n::Int) where {T} = $GRID{T}(n)
    @eval $GRID(n::Int) = $GRID{Float64}(n)
    @eval $GRID(n::Int, d::AbstractInterval) =
        $GRID(n, endpoints(d)...)
    @eval $GRID(n::Int, a, b) = rescale($GRID{typeof((b-a)/n)}(n), a, b)
end

# extendable grids
_extension_size(::PeriodicEquispacedGrid, n::Int, factor::Int) = factor*n
_extension_size(::FourierGrid, n::Int, factor::Int) = factor*n
_extension_size(::EquispacedGrid, n::Int, factor::Int) = factor*n-1

for GRID in (:PeriodicEquispacedGrid,:FourierGrid,:EquispacedGrid)
    @eval hasextension(::$GRID) = true
    @eval extend(grid::$GRID, factor::Int) =
        resize(grid, _extension_size(grid, length(grid), factor))
end

# function mapped_grid(grid::FourierGrid{T}, map::AffineMap) where T
#     s = map*support(grid)
#     s≈UnitInterval{T}() ?
#         grid : PeriodicEquispacedGrid{T}(length(grid), endpoints(s)...)
# end

function rescale(g::FourierGrid, a, b)
	m = interval_map(endpoints(support(g))..., a, b)
	mapped_grid(g, m)
end

mapped_grid(g::FourierGrid, map::AffineMap) = MappedGrid(g, map)


text:
#script
## Binary  estimates
    # estSS_1GW,estGas_1GW = 0, 0
    # eMidMod_1GW = DynNets.SdErgmDirBin1(degsIO_T, [ zeros(2N) ,0.9  * ones(1), 0.01 * ones(1)], [Int.(ones(N)),Int.(ones(1)),Int.(ones(1))],"")
    # @time estSS_1GW =  DynNets.estimateSnapSeq(eMidMod_1GW)
    # @time estGas_1GW = DynNets.estimate(eMidMod_1GW)
    # eMidMod = DynNets.SdErgmDirBin1(degsIO_T)
    # estSS,estGas,estGasTarg = 0, 0,0
    # estSS =  DynNets.estimateSnapSeq(eMidMod)
    # #estGasTarg = DynNets.estimateTargeting(eMidMod)
    # estGas = DynNets.estimate(eMidMod)
    # # save Bin
    # save_path = save_fold*file_nameStart* "_Bin1_" *file_nameEnd#
    # @save(save_path,estSS,estGas,estGasTarg,estSS_1GW,estGas_1GW,AeMidWeekly_T, YeMidWeekly_T ,weekInd,degsIO_T)
    #
## Weighted Estimates

estSSW1 =  DynNets.estimateSnapSeq(eMidModW1)
estMean = DynNets.estSingSnap(eMidModW1,meanSq(strIO_T,1))

prod(isfinite.(estSSW1[:,2:end]),2)
parest_t=StaticNets.estimate(fooErgmDirW1; strIO = strIO_T[20,:] ,groupsInds = eMidModW1.groupsInds[1])
strIO_T[20,N+1:2N]
unique(parest_t[1])
prod(parest_t[1][N+1:2N] .> 0)
indmin(parest_t[1][N+1:2N])
log( parest_t[1][N+1]  + parest_t[1][1:N] )
StaticNets.bndPar2uBndPar(fooErgmDirW1,parest_t[1])

using StatsFuns
smallVal = 1e-15
log(bigConstVal - smallVal)  - log(bigConstVal) + log1p(smallVal/bigConstVal)

##
save_path = save_fold*file_nameStart* "_W1_" *file_nameEnd#
@load(save_path,estSSW1,estGasW1,estGasTargW1,estSSW1_1GW,estGasW1_1GW,AeMidWeekly_T, YeMidWeekly_T ,weekInd,degsIO_T)

# test

B = ones(1)* 0.15999999999999998
A = ones(1)*0.0003000000000000001
uBndPar = StaticNets.bndPar2uBndPar(StaticNets.fooErgmDirW1, meanSq( DynNets.uBndPar2bndPar_T(eMidModW1,estSSW1),1) )
#StaticNets.bndPar2uBndPar(StaModType(eMidModW1),meanSq(estSSW1,1))
W = uBndPar.*(1-B)
dgpParArr = [W,B,A];dgpParVec = [W;B;A]
filterParTestUbnd,testlike = DynNets.score_driven_filter_or_dgp(eMidModW1,dgpParVec)

filterParTest = DynNets.uBndPar2bndPar_T(eMidModW1,filterParTestUbnd)
 pin = plot(filterParTest[:,1:N][:,meanSq(filterParTest[:,1:N],1).<1e5]);
 pout = plot(filterParTest[:,1+N:2N][:,meanSq(filterParTest[:,1+N:2N],1).<1e5]);
 plot(pin,pout,layout=(2,1),size = (1200,600))

##
targMean =meanSq( estSSW1[1:Tweeks,:],1)# StaticNets.bndPar2uBndPar(StaticNets.fooErgmDirW1, meanSq( DynNets.uBndPar2bndPar_T(eMidModW1,estSSW1),1) )

 estGasTargW1,x,fVal = DynNets.estimate(eMidModW1;targeting = true, meanEstSS = targMean)
 println(estGasTargW1[2],estGasTargW1[3])
 estGasW1,estSSW1_1GW,estGasW1_1GW = 0,0,0

estParVec = [estGasTargW1[1];estGasTargW1[2];estGasTargW1[3]]
 filterParTestUbndEst,testlike = DynNets.score_driven_filter_or_dgp(eMidModW1,estParVec)
 filterParTest = filterParTestUbndEst#DynNets.uBndPar2bndPar_T(eMidModW1,filterParTestUbndEst)#
 pin = plot(filterParTest[:,1:N][:,meanSq(filterParTest[:,1:N],1).<1e5]);
 pout = plot(filterParTest[:,1+N:2N][:,meanSq(filterParTest[:,1+N:2N],1).<1e5]);
 plot(pin,pout,layout=(2,1),size = (1200,600))

#@save(save_path,estSSW1,estGasW1,estGasTargW1,estSSW1_1GW,estGasW1_1GW,AeMidWeekly_T, YeMidWeekly_T ,weekInd,degsIO_T)

##

max_stars_repo_name: domenicodigangi/SDERGM

max_issues_repo_name: domenicodigangi/SDERGM

max_forks_repo_name: domenicodigangi/SDERGM


text:
"""
    μ = center!(X[, μ, obsdim])

or

    μ = center!(D[, colnames, μ])

where `X` is of type Matrix or Vector and `D` of type DataFrame.

Center `X` along `obsdim` around the corresponding entry in the
vector `μ`. If `μ` is not specified then it defaults to the
feature specific means.

For DataFrames, `obsdim` is obsolete and centering is done column wise.
Instead the vector `colnames` allows to specify which columns to center.
If `colnames` is not provided all columns of type T<:Real are centered.

Example:

    X = rand(4, 100)
    D = DataFrame(A=rand(10), B=collect(1:10), C=[string(x) for x in 1:10])

    μ = center!(X, obsdim=2)
    μ = center!(X, ObsDim.First())
    μ = center!(D)
    μ = center!(D, [:A, :B])

"""
function center!(X, μ; obsdim=LearnBase.default_obsdim(X))
    center!(X, μ, convert(ObsDimension, obsdim))
end

function center!(X; obsdim=LearnBase.default_obsdim(X))
    center!(X, convert(ObsDimension, obsdim))
end

function center!(X::AbstractArray{T,N}, μ::AbstractVector, ::ObsDim.Last) where {T,N}
    center!(X, μ, ObsDim.Constant{N}())
end

function center!(X::AbstractArray{T,N}, ::ObsDim.Last) where {T,N}
    center!(X, ObsDim.Constant{N}())
end

function center!(X::AbstractArray{T,N}, obsdim::ObsDim.Constant{M}) where {T,N,M}
    μ = vec(mean(X, dims=M))
    center!(X, μ, obsdim)
end

function center!(X::AbstractVector{T}, ::ObsDim.Constant{1}) where T
    μ = mean(X)
    for i in 1:length(X)
        X[i] = X[i] - μ
    end
    μ
end

function center!(X::AbstractVector, μ::AbstractVector, ::ObsDim.Constant{1})
    @inbounds for i in 1:length(X)
        X[i] = X[i] - μ[i]
    end
    μ
end

function center!(X::AbstractVector, μ::AbstractFloat, ::ObsDim.Constant{1})
    @inbounds for i in 1:length(X)
        X[i] = X[i] - μ
    end
    μ
end

function center!(X::AbstractMatrix, μ::AbstractVector, ::ObsDim.Constant{1})
    nObs, nVars = size(X)
    for iVar in 1:nVars
        @inbounds for iObs in 1:nObs
            X[iObs, iVar] = X[iObs, iVar] - μ[iVar]
        end
    end
    μ
end

function center!(X::AbstractMatrix, μ::AbstractVector, ::ObsDim.Constant{2})
    nVars, nObs = size(X)
    for iObs in 1:nObs
        @inbounds for iVar in 1:nVars
            X[iVar, iObs] = X[iVar, iObs] - μ[iVar]
        end
    end
    μ
end

function center!(D::AbstractDataFrame)
    μ_vec = Float64[]

    flt = Bool[T <: Real for T in eltype.(eachcol(D))]
    for colname in propertynames(D)[flt]
        μ = mean(D[!, colname])
        center!(D, colname, μ)
        push!(μ_vec, μ)
    end
    μ_vec
end

function center!(D::AbstractDataFrame, colnames::AbstractVector{Symbol})
    μ_vec = Float64[]
    for colname in colnames
        if eltype(D[!, colname]) <: Real
            μ = mean(D[!, colname])
            if ismissing(μ)
                @warn("Column \"$colname\" contains missing values, skipping rescaling of this column!")
                continue
            end
            center!(D, colname, μ)
            push!(μ_vec, μ)
        else
            @warn("Skipping \"$colname\", centering only valid for columns of type T <: Real.")
        end
    end
    μ_vec
end

function center!(D::AbstractDataFrame, colnames::AbstractVector{Symbol}, μ::AbstractVector)
    for (icol, colname) in enumerate(colnames)
        if eltype(D[!, colname]) <: Real
            center!(D, colname, μ[icol])
        else
            @warn("Skipping \"$colname\", centering only valid for columns of type T <: Real.")
        end
    end
    μ
end

function center!(D::AbstractDataFrame, colname::Symbol, μ)
    if any(ismissing, D[!, colname])
        @warn("Column \"$colname\" contains missing values, skipping centering on this column!")
    else
        newcol::Vector{Float64} = convert(Vector{Float64}, D[!, colname])
        nobs = length(newcol)
        @inbounds for i in eachindex(newcol)
            newcol[i] -= μ
        end
        D[!, colname] = newcol
    end
    μ
end

"""
    μ, σ = rescale!(X[, μ, σ, obsdim])

or

    μ, σ = rescale!(D[, colnames, μ, σ])

where `X` is of type Matrix or Vector and `D` of type DataFrame.

Center `X` along `obsdim` around the corresponding entry in the
vector `μ` and then rescale each feature using the corresponding
entry in the vector `σ`.

For DataFrames, `obsdim` is obsolete and centering is done column wise.
The vector `colnames` allows to specify which columns to center.
If `colnames` is not provided all columns of type T<:Real are centered.

Example:

    X = rand(4, 100)
    D = DataFrame(A=rand(10), B=collect(1:10), C=[string(x) for x in 1:10])

    μ, σ = rescale!(X, obsdim=2)
    μ, σ = rescale!(X, ObsDim.First())
    μ, σ = rescale!(D)
    μ, σ = rescale!(D, [:A, :B])

"""
function rescale!(X, μ, σ; obsdim=LearnBase.default_obsdim(X))
    rescale!(X, μ, σ, convert(ObsDimension, obsdim))
end

function rescale!(X::AbstractArray{T,N}, μ, σ, ::ObsDim.Last) where {T,N}
    rescale!(X, μ, σ, ObsDim.Constant{N}())
end

function rescale!(X; obsdim=LearnBase.default_obsdim(X))
    rescale!(X, convert(ObsDimension, obsdim))
end

function rescale!(X::AbstractArray{T,N}, ::ObsDim.Last) where {T,N}
    rescale!(X, ObsDim.Constant{N}())
end

function rescale!(X::AbstractArray{T,N}, obsdim::ObsDim.Constant{M}) where {T,N,M}
    μ = vec(mean(X, dims=M))
    σ = vec(std(X, dims=M))
    rescale!(X, μ, σ, obsdim)
end

function rescale!(X::AbstractVector, ::ObsDim.Constant{1})
    μ = mean(X)
    σ = std(X)
    @inbounds for i in 1:length(X)
        X[i] = (X[i] - μ) / σ
    end
    μ, σ
end

function rescale!(X::AbstractMatrix, μ::AbstractVector, σ::AbstractVector, ::ObsDim.Constant{2})
    σ[σ .== 0] .= 1
    nVars, nObs = size(X)
    for iObs in 1:nObs
        @inbounds for iVar in 1:nVars
            X[iVar, iObs] = (X[iVar, iObs] - μ[iVar]) / σ[iVar]
        end
    end
    μ, σ
end

function rescale!(X::AbstractMatrix, μ::AbstractVector, σ::AbstractVector, ::ObsDim.Constant{1})
    σ[σ .== 0] .= 1
    nObs, nVars = size(X)
    for iVar in 1:nVars
        @inbounds for iObs in 1:nObs
            X[iObs, iVar] = (X[iObs, iVar] - μ[iVar]) / σ[iVar]
        end
    end
    μ, σ
end

function rescale!(X::AbstractVector, μ::AbstractVector, σ::AbstractVector, ::ObsDim.Constant{1})
    @inbounds for i in 1:length(X)
        X[i] = (X[i] - μ[i]) / σ[i]
    end
    μ, σ
end

function rescale!(X::AbstractVector, μ::AbstractFloat, σ::AbstractFloat, ::ObsDim.Constant{1})
    @inbounds for i in 1:length(X)
        X[i] = (X[i] - μ) / σ
    end
    μ, σ
end

function rescale!(D::AbstractDataFrame)
    μ_vec = Float64[]
    σ_vec = Float64[]

    flt = Bool[T <: Real for T in eltype.(eachcol(D))]
    for colname in propertynames(D)[flt]
        μ = mean(D[!, colname])
        σ = std(D[!, colname])
        rescale!(D, colname, μ, σ)
        push!(μ_vec, μ)
        push!(σ_vec, σ)
    end
    μ_vec, σ_vec
end

function rescale!(D::AbstractDataFrame, colnames::Vector{Symbol})
    μ_vec = Float64[]
    σ_vec = Float64[]
    for colname in colnames
        if eltype(D[!, colname]) <: Real
            μ = mean(D[!, colname])
            σ = std(D[!, colname])
            if ismissing(μ)
                @warn("Column \"$colname\" contains missing values, skipping rescaling of this column!")
                continue
            end
            rescale!(D, colname, μ, σ)
            push!(μ_vec, μ)
            push!(σ_vec, σ)
        else
            @warn("Skipping \"$colname\", rescaling only valid for columns of type T <: Real.")
        end
    end
    μ_vec, σ_vec
end

function rescale!(D::AbstractDataFrame, colnames::Vector{Symbol}, μ::AbstractVector, σ::AbstractVector)
    for (icol, colname) in enumerate(colnames)
        if eltype(D[!, colname]) <: Real
            rescale!(D, colname, μ[icol], σ[icol])
        else
            @warn("Skipping \"$colname\", rescaling only valid for columns of type T <: Real.")
        end
    end
    μ, σ
end

function rescale!(D::AbstractDataFrame, colname::Symbol, μ, σ)
    if any(ismissing, D[!, colname])
        @warn("Column \"$colname\" contains missing values, skipping rescaling of this column!")
    else
        σ_div = σ == 0 ? one(σ) : σ
        newcol::Vector{Float64} = convert(Vector{Float64}, D[!, colname])
        nobs = length(newcol)
        @inbounds for i in eachindex(newcol)
            newcol[i] = (newcol[i] - μ) / σ_div
        end
        D[!, colname] = newcol
    end
    μ, σ
end

struct FeatureNormalizer
    offset::Vector{Float64}
    scale::Vector{Float64}

    function FeatureNormalizer(offset::Vector{Float64}, scale::Vector{Float64})
        @assert length(offset) == length(scale)
        new(offset, scale)
    end
end

function FeatureNormalizer(X::AbstractMatrix{T}) where T<:Real
    FeatureNormalizer(vec(mean(X, dims=2)), vec(std(X, dims=2)))
end

function StatsBase.fit(::Type{FeatureNormalizer}, X::AbstractMatrix{T}) where T<:Real
    FeatureNormalizer(X)
end

function StatsBase.predict!(cs::FeatureNormalizer, X::AbstractMatrix{T}) where T<:Real
    @assert length(cs.offset) == size(X, 1)
    rescale!(X, cs.offset, cs.scale)
    X
end

function StatsBase.predict(cs::FeatureNormalizer, X::AbstractMatrix{T}) where T<:AbstractFloat
    Xnew = copy(X)
    StatsBase.predict!(cs, Xnew)
end

function StatsBase.predict(cs::FeatureNormalizer, X::AbstractMatrix{T}) where T<:Real
    X = convert(AbstractMatrix{AbstractFloat}, X)
    StatsBase.predict!(cs, X)
end


text:
using InteractUIkit, Observables, Colors, DataStructures
import InteractBase: scope
@static if VERSION < v"0.7.0-DEV.2005"
    using Base.Test
else
    using Test
end

@testset "input widgets" begin
    a = toggle(label = "Agreed")
    @test observe(a)[] == false
    s = Observable(true)
    a = toggle(s, label = "Agreed")
    @test observe(a)[] == true
end


text:
using OrdinaryDiffEq, ModelingToolkit, Test
using GalacticOptim, RecursiveArrayTools, GalacticOptimJL

N = 32
const xyd_brusselator = range(0, stop = 1, length = N)
brusselator_f(x, y, t) = (((x - 0.3)^2 + (y - 0.6)^2) <= 0.1^2) * (t >= 1.1) * 5.0
limit(a, N) = ModelingToolkit.ifelse(a == N + 1, 1, ModelingToolkit.ifelse(a == 0, N, a))
function brusselator_2d_loop(du, u, p, t)
    A, B, alpha, dx = p
    alpha = alpha / dx^2
    @inbounds for I in CartesianIndices((N, N))
        i, j = Tuple(I)
        x, y = xyd_brusselator[I[1]], xyd_brusselator[I[2]]
        ip1, im1, jp1, jm1 = limit(i + 1, N), limit(i - 1, N), limit(j + 1, N),
                             limit(j - 1, N)
        du[i, j, 1] = alpha * (u[im1, j, 1] + u[ip1, j, 1] + u[i, jp1, 1] + u[i, jm1, 1] -
                       4u[i, j, 1]) +
                      B + u[i, j, 1]^2 * u[i, j, 2] - (A + 1) * u[i, j, 1] +
                      brusselator_f(x, y, t)
        du[i, j, 2] = alpha * (u[im1, j, 2] + u[ip1, j, 2] + u[i, jp1, 2] + u[i, jm1, 2] -
                       4u[i, j, 2]) +
                      A * u[i, j, 1] - u[i, j, 1]^2 * u[i, j, 2]
    end
end

# Test with tuple parameters
p = (3.4, 1.0, 10.0, step(xyd_brusselator))

function init_brusselator_2d(xyd)
    N = length(xyd)
    u = zeros(N, N, 2)
    for I in CartesianIndices((N, N))
        x = xyd[I[1]]
        y = xyd[I[2]]
        u[I, 1] = 22 * (y * (1 - y))^(3 / 2)
        u[I, 2] = 27 * (x * (1 - x))^(3 / 2)
    end
    u
end
u0 = init_brusselator_2d(xyd_brusselator)

# Test with 3-tensor inputs
prob_ode_brusselator_2d = ODEProblem(brusselator_2d_loop,
                                     u0, (0.0, 11.5), p)

modelingtoolkitize(prob_ode_brusselator_2d)

## Optimization

rosenbrock(x, p) = (p[1] - x[1])^2 + p[2] * (x[2] - x[1]^2)^2
x0 = zeros(2)
p = [1.0, 100.0]

prob = OptimizationProblem(rosenbrock, x0, p)
sys = modelingtoolkitize(prob) # symbolicitize me captain!

prob = OptimizationProblem(sys, x0, p, grad = true, hess = true)
sol = solve(prob, NelderMead())
@test sol.minimum < 1e-8

sol = solve(prob, BFGS())
@test sol.minimum < 1e-8

sol = solve(prob, Newton())
@test sol.minimum < 1e-8

## SIR System Regression Test

β = 0.01# infection rate
λ_R = 0.05 # inverse of transition time from  infected to recovered
λ_D = 0.83 # inverse of transition time from  infected to dead
i₀ = 0.075 # fraction of initial infected people in every age class
𝒫 = vcat([β, λ_R, λ_D]...)

# regional contact matrix and regional population

## regional contact matrix
regional_all_contact_matrix = [3.45536 0.485314 0.506389 0.123002;
                               0.597721 2.11738 0.911374 0.323385;
                               0.906231 1.35041 1.60756 0.67411;
                               0.237902 0.432631 0.726488 0.979258] # 4x4 contact matrix

## regional population stratified by age
N = [723208, 874150, 1330993, 1411928] # array of 4 elements, each of which representing the absolute amount of population in the corresponding age class.

# Initial conditions
I₀ = repeat([i₀], 4)
S₀ = N .- I₀
R₀ = [0.0 for n in 1:length(N)]
D₀ = [0.0 for n in 1:length(N)]
D_tot₀ = [0.0 for n in 1:length(N)]
ℬ = vcat([S₀, I₀, R₀, D₀, D_tot₀]...)

# Time
final_time = 20
𝒯 = (1.0, final_time);

function SIRD_ac!(du, u, p, t)
    # Parameters to be calibrated
    β, λ_R, λ_D = p

    # initialize this parameter (death probability stratified by age, taken from literature)

    δ₁, δ₂, δ₃, δ₄ = [
        0.003 / 100,
        0.004 / 100,
        (0.015 + 0.030 + 0.064 + 0.213 + 0.718) / (5 * 100),
        (2.384 + 8.466 + 12.497 + 1.117) / (4 * 100),
    ]
    δ = vcat(repeat([δ₁], 1), repeat([δ₂], 1), repeat([δ₃], 1), repeat([δ₄], 4 - 1 - 1 - 1))

    C = regional_all_contact_matrix

    # State variables
    S = @view u[(4 * 0 + 1):(4 * 1)]
    I = @view u[(4 * 1 + 1):(4 * 2)]
    R = @view u[(4 * 2 + 1):(4 * 3)]
    D = @view u[(4 * 3 + 1):(4 * 4)]
    D_tot = @view u[(4 * 4 + 1):(4 * 5)]

    # Differentials
    dS = @view du[(4 * 0 + 1):(4 * 1)]
    dI = @view du[(4 * 1 + 1):(4 * 2)]
    dR = @view du[(4 * 2 + 1):(4 * 3)]
    dD = @view du[(4 * 3 + 1):(4 * 4)]
    dD_tot = @view du[(4 * 4 + 1):(4 * 5)]

    # Force of infection
    Λ = β * [sum([C[i, j] * I[j] / N[j] for j in 1:size(C)[1]]) for i in 1:size(C)[2]]

    # System of equations
    @. dS = -Λ * S
    @. dI = Λ * S - ((1 - δ) * λ_R + δ * λ_D) * I
    @. dR = λ_R * (1 - δ) * I
    @. dD = λ_D * δ * I
    @. dD_tot = dD[1] + dD[2] + dD[3] + dD[4]
end;

# create problem and check it works
problem = ODEProblem(SIRD_ac!, ℬ, 𝒯, 𝒫)
@time solution = solve(problem, Tsit5(), saveat = 1:final_time);

problem = ODEProblem(SIRD_ac!, ℬ, 𝒯, 𝒫)
sys = modelingtoolkitize(problem)
fast_problem = ODEProblem(sys, ℬ, 𝒯, 𝒫)
@time solution = solve(fast_problem, Tsit5(), saveat = 1:final_time)

## Issue #778

r0 = [1131.340, -2282.343, 6672.423]
v0 = [-5.64305, 4.30333, 2.42879]
Δt = 86400.0 * 365
μ = 398600.4418
rv0 = ArrayPartition(r0, v0)

f = function (dy, y, μ, t)
    r = sqrt(sum(y[1, :] .^ 2))
    dy[1, :] = y[2, :]
    dy[2, :] = -μ .* y[1, :] / r^3
end

prob = ODEProblem(f, rv0, (0.0, Δt), μ)
modelingtoolkitize(prob)

# Index reduction and mass matrix handling
using LinearAlgebra
function pendulum!(du, u, p, t)
    x, dx, y, dy, T = u
    g, L = p
    du[1] = dx
    du[2] = T * x
    du[3] = dy
    du[4] = T * y - g
    du[5] = x^2 + y^2 - L^2
    return nothing
end
pendulum_fun! = ODEFunction(pendulum!, mass_matrix = Diagonal([1, 1, 1, 1, 0]))
u0 = [1.0, 0, 0, 0, 0]
p = [9.8, 1]
tspan = (0, 10.0)
pendulum_prob = ODEProblem(pendulum_fun!, u0, tspan, p)
pendulum_sys_org = modelingtoolkitize(pendulum_prob)
sts = states(pendulum_sys_org)
pendulum_sys = dae_index_lowering(pendulum_sys_org)
prob = ODEProblem(pendulum_sys, Pair[], tspan)
sol = solve(prob, Rodas4())
l2 = sol[sts[1]] .^ 2 + sol[sts[3]] .^ 2
@test all(l -> abs(sqrt(l) - 1) < 0.05, l2)

ff911 = (du, u, p, t) -> begin
    du[1] = u[2] + 1.0
    du[2] = u[1] - 1.0
end
prob = ODEProblem(ff911, zeros(2), (0, 1.0))
@test_nowarn modelingtoolkitize(prob)

k(x, p, t) = p * x
x0 = 1.0
p = 0.98
tspan = (0.0, 1.0)
prob = ODEProblem(k, x0, tspan, p)
sys = modelingtoolkitize(prob)

k(x, p, t) = 0.98 * x
x0 = 1.0
tspan = (0.0, 1.0)
prob = ODEProblem(k, x0, tspan)
sys = modelingtoolkitize(prob)

## https://github.com/SciML/ModelingToolkit.jl/issues/1054
using LabelledArrays
using ModelingToolkit

# ODE model: simple SIR model with seasonally forced contact rate
function SIR!(du, u, p, t)

    # states
    (S, I, R) = u[1:3]
    N = S + I + R

    # params
    β = p.β
    η = p.η
    φ = p.φ
    ω = 1.0 / p.ω
    μ = p.μ
    σ = p.σ

    # FOI
    βeff = β * (1.0 + η * cos(2.0 * π * (t - φ) / 365.0))
    λ = βeff * I / N

    # change in states
    du[1] = (μ * N - λ * S - μ * S + ω * R)
    du[2] = (λ * S - σ * I - μ * I)
    du[3] = (σ * I - μ * R - ω * R)
    du[4] = (σ * I) # cumulative incidence
end

# Solver settings
tmin = 0.0
tmax = 10.0 * 365.0
tspan = (tmin, tmax)

# Initiate ODE problem
theta_fix = [1.0 / (80 * 365)]
theta_est = [0.28, 0.07, 1.0 / 365.0, 1.0, 1.0 / 5.0]
p = @LArray [theta_est; theta_fix] (:β, :η, :ω, :φ, :σ, :μ)
u0 = @LArray [9998.0, 1.0, 1.0, 1.0] (:S, :I, :R, :C)

# Initiate ODE problem
problem = ODEProblem(SIR!, u0, tspan, p)
sys = modelingtoolkitize(problem)

@parameters t
@test all(isequal.(parameters(sys), getproperty.(@variables(β, η, ω, φ, σ, μ), :val)))
@test all(isequal.(Symbol.(states(sys)), Symbol.(@variables(S(t), I(t), R(t), C(t)))))

# https://github.com/SciML/ModelingToolkit.jl/issues/1158

function ode_prob(du, u, p::NamedTuple, t)
    du[1] = u[1] + p.α * u[2]
    du[2] = u[2] + p.β * u[1]
end
params = (α = 1, β = 1)
prob = ODEProblem(ode_prob, [1 1], (0, 1), params)
sys = modelingtoolkitize(prob)
@test nameof.(parameters(sys)) == [:α, :β]

function ode_prob(du, u, p::Tuple, t)
    α, β = p
    du[1] = u[1] + α * u[2]
    du[2] = u[2] + β * u[1]
end

params = (1, 1)
prob = ODEProblem(ode_prob, [1 1], (0, 1), params)
sys = modelingtoolkitize(prob)
@test nameof.(parameters(sys)) == [:α₁, :α₂]


text:
standartize(x; dims=2) = reshape(x, dims, :)

loss_x(ux, p) = L2dist(K_all(standartize(ux), p.wp_x, p.vs_x) .* p.w_x, p.k_x .* p.w_x)
loss_y(uy, p) = L2dist(K_all(standartize(uy), p.wp_y, p.vs_y) .* p.w_y, p.k_y .* p.w_y)
loss_y_cross(uy, ux, p) = L2dist(K_all(standartize(ux), standartize(uy), p.wp_xy, p.vs_x_cross, p.vs_y_cross) .* p.w_y_cross, p.k_y_cross .* p.w_y_cross)


text:
function write_csv(names, data, filename)
    open(joinpath(@__DIR__, filename), "w") do f
        write(f, join(names, ",") * "\n")
        for tup in data
            write(f, join(tup, ",") * "\n")
        end
    end
    return nothing
end

function normal_sample(μ, Σ)
    n = length(μ)
    x = μ + sqrt(Σ) * randn(n)
    return x
end


text:
function _inverse_retract_nlsolve(
    M::Manifold,
    p,
    q,
    retraction,
    X0,
    project_tangent,
    project_point,
    nlsolve_kwargs;
    kwargs...,
)
    function f!(F, X)
        project_tangent && project!(M, X, p, X)
        retract!(M, F, p, project(M, p, X), retraction; kwargs...)
        project_point && project!(M, q, q)
        F .-= q
        return F
    end
    res = NLsolve.nlsolve(f!, X0; nlsolve_kwargs...)
    return res
end


text:
module QuadraticRelaxations

import JuMP, MathOptInterface
import LinearAlgebra, SparseArrays
import PiecewiseLinearOpt

const MOI = MathOptInterface
const MOIU = MOI.Utilities
const VI = MOI.VariableIndex
const SV = MOI.SingleVariable
const SAT = MOI.ScalarAffineTerm{Float64}
const SAF = MOI.ScalarAffineFunction{Float64}
const SQT = MOI.ScalarQuadraticTerm{Float64}
const SQF = MOI.ScalarQuadraticFunction{Float64}
const VAT = MOI.VectorAffineTerm{Float64}
const VAF = MOI.VectorAffineFunction{Float64}
const VQT = MOI.VectorQuadraticTerm{Float64}
const VQF = MOI.VectorQuadraticFunction{Float64}

abstract type AbstractRelaxation end

abstract type AbstractShiftMethod end

include("moi_util.jl")
include("reformulation.jl")
include("shift.jl")

include("relaxations/exact.jl")
include("relaxations/secant.jl")
include("relaxations/cda.jl")
include("relaxations/neural_net.jl")
include("relaxations/standard_nmdt.jl")
include("relaxations/tightened_nmdt.jl")

end # module


text:
import NearestNeighbors
import NearestNeighbors: KDTree

export KDTree

###########################################################################################
# Standard API
###########################################################################################
function Neighborhood.searchstructure(::Type{KDTree}, data, metric; kwargs...)
    return KDTree(data, metric; kwargs...)
end

datatype(::Type{<:KDTree{V}}) where V = V
getmetric(tree::KDTree) = tree.metric

function Neighborhood.search(tree::KDTree, query, t::NeighborNumber, skip=alwaysfalse; sortds=true)
    return NearestNeighbors.knn(tree, query, t.k, sortds, skip)
end

function Neighborhood.search(tree::KDTree, query, t::WithinRange, skip=alwaysfalse; sortds=true)
    idxs = NearestNeighbors.inrange(tree, query, t.r)
    skip ≠ alwaysfalse && filter!(!skip, idxs)
    ds = _NN_get_ds(tree, query, idxs)
    if sortds # sort according to distances
        sp = sortperm(ds)
        sort!(ds)
        idxs = idxs[sp]
    end
    return idxs, ds
end

function _NN_get_ds(tree::KDTree, query, idxs)
    if tree.reordered
        ds = [
            evaluate(tree.metric, query, tree.data[
                findfirst(isequal(i), tree.indices)
            ]) for i in idxs]
    else
        ds = [evaluate(tree.metric, query, tree.data[i]) for i in idxs]
    end
end

# Performance method when distances are not required (can't sort then)
function Neighborhood.isearch(tree::KDTree, query, t::WithinRange, skip=alwaysfalse; sortds=false)
    idxs = NearestNeighbors.inrange(tree, query, t.r, sortds)
    skip ≠ alwaysfalse && filter!(!skip, idxs)
    return idxs
end

###########################################################################################
# Bulk and skip predicates
###########################################################################################
function Neighborhood.bulksearch(tree::KDTree, queries, t::NeighborNumber; sortds=true)
    return NearestNeighbors.knn(tree, queries, t.k, sortds)
end

# no easy skip version for knn bulk search
function Neighborhood.bulksearch(tree::KDTree, queries, t::NeighborNumber, skip; sortds=true)
    k, N = t.k, length(queries)
    dists = [Vector{eltype(queries[1])}(undef, k) for _ in 1:N]
    idxs = [Vector{Int}(undef, k) for _ in 1:N]
    for j in 1:N
        # Notice that this `_skip` definition matches our API definition!
        _skip = i -> skip(i, j)
        @inbounds NearestNeighbors.knn_point!(tree, queries[j], sortds, dists[j], idxs[j], _skip)
    end
    return idxs, dists
end

# but there is an easy skip version for inrange bulk isearch!
function Neighborhood.bulkisearch(tree::KDTree, queries, t::WithinRange, skip=alwaysfalse)
    vec_of_idxs = NearestNeighbors.inrange(tree, queries, t.r)
    skip ≠ alwaysfalse && vecskipfilter!(vec_of_idxs, skip)
    return vec_of_idxs
end

function Neighborhood.bulksearch(tree::KDTree, queries, t::WithinRange, skip=alwaysfalse; sortds=true)
    vec_of_idxs = NearestNeighbors.inrange(tree, queries, t.r)
    vec_of_ds = [ _NN_get_ds(tree, queries[j], vec_of_idxs[j]) for j in 1:length(queries)]
    skip ≠ alwaysfalse && vecskipfilter!(vec_of_idxs, vec_of_ds, skip)
    if sortds # sort according to distances
        for i in 1:length(queries)
            @inbounds ds = vec_of_ds[i]
            length(ds) ≤ 1 && continue
            sp = sortperm(ds)
            sort!(ds)
            vec_of_idxs[i] = vec_of_idxs[i][sp]
        end
    end
    return vec_of_idxs, vec_of_ds
end

function vecskipfilter!(vec_of_idxs, skip)
    for j in 1:length(vec_of_idxs)
        @inbounds idxs = vec_of_idxs[j]
        filter!(i -> !skip(i, j), idxs)
    end
    return vec_of_idxs
end
function vecskipfilter!(vec_of_idxs, vec_of_ds, skip)
    for j in 1:length(vec_of_idxs)
        @inbounds idxs = vec_of_idxs[j]
        todelete = [i for i in 1:length(idxs) if skip(idxs[i], j)]
        deleteat!(idxs, todelete)
        deleteat!(vec_of_ds[j], todelete)
    end
    return vec_of_idxs
end


text:
nconstraints(f::Union{SVF, SAF}, s) = 1
nconstraints(f::VVF, s) = length(f.variables)

function _allocate_constraint(m::SOItoMOIBridge, f, s)
    ci = CI{typeof(f), typeof(s)}(m.nconstrs)
    n = nconstraints(f, s)
    # Fails on Julia v0.6
    #m.constrmap[ci] = m.nconstrs .+ (1:n)
    m.constrmap[ci] = (m.nconstrs + 1):(m.nconstrs + n)
    m.nconstrs += n
    return ci
end
function MOIU.allocate_constraint(m::SOItoMOIBridge, f::SAF, s::SupportedSets)
    _allocate_constraint(m::SOItoMOIBridge, f, s)
end

function loadcoefficients!(m::SOItoMOIBridge, cs::UnitRange,
                           f::MOI.ScalarAffineFunction, s)
    f = MOIU.canonical(f) # sum terms with same variables and same outputindex
    @assert length(cs) == 1
    c = first(cs)
    rhs = MOIU.getconstant(s) - MOI._constant(f)
    for t in f.terms
        if !iszero(t.coefficient)
            for (blk, i, j, coef, shift) in varmap(m, t.variable_index)
                if !iszero(blk)
                    @assert !iszero(coef)
                    setconstraintcoefficient!(m.sdoptimizer, t.coefficient*coef, c, blk, i, j)
                end
                rhs -= t.coefficient * shift
            end
        end
    end
    setconstraintconstant!(m.sdoptimizer, rhs, c)
end

function MOIU.load_constraint(m::SOItoMOIBridge, ci::CI, f::SAF, s::SupportedSets)
    setconstant!(m, ci, s)
    cs = m.constrmap[ci]
    @assert !isempty(cs)
    loadcoefficients!(m, cs, f, s)
end


text:
using Base.Test
using StatsFuns
using PMTK4Stats

a = 3
b = 9
icdf(p::Real) = betainvcdf(a, b, p)
lb, ub = hdiFromIcdf(icdf)

@test_approx_eq_eps .04 lb 1e-2
@test_approx_eq_eps .48 ub 1e-2


text:
using PyPlot
using GeoDatasets

res = 'c'

# coastline
for (lon,lat) in GeoDatasets.gshhg(res,[1,5])
    plot(lon,lat,"-",color="k")
end

# rivers
for (lon,lat) in GeoDatasets.WDBII_river(res,[1,2])
    plot(lon,lat,"b-",alpha = 0.5, linewidth = 0.5)
end

# borders
for (lon,lat) in GeoDatasets.WDBII_border(res,1)
    plot(lon,lat,"g-",alpha = 0.8, linewidth = 0.3)
end

max_stars_repo_name: JuliaGeo/GeoDatasets

max_issues_repo_name: JuliaGeo/GeoDatasets

max_forks_repo_name: JuliaGeo/GeoDatasets


text:
# Julia wrapper for header: /usr/local/include/libavfilter/avfiltergraph.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

export
    avfilter_graph_alloc,
    avfilter_graph_get_filter,
    avfilter_graph_add_filter,
    avfilter_graph_create_filter,
    avfilter_graph_config,
    avfilter_graph_free,
    avfilter_graph_parse

function avfilter_graph_alloc()
    ccall((:avfilter_graph_alloc,libavfilter),Ptr{AVFilterGraph},())
end

function avfilter_graph_get_filter(graph,name)
    ccall((:avfilter_graph_get_filter,libavfilter),Ptr{AVFilterContext},(Ptr{AVFilterGraph},Ptr{UInt8}),graph,name)
end

function avfilter_graph_add_filter(graphctx,filter)
    ccall((:avfilter_graph_add_filter,libavfilter),Cint,(Ptr{AVFilterGraph},Ptr{AVFilterContext}),graphctx,filter)
end

function avfilter_graph_create_filter(filt_ctx,filt,name,args,opaque,graph_ctx)
    ccall((:avfilter_graph_create_filter,libavfilter),Cint,(Ptr{Ptr{AVFilterContext}},Ptr{AVFilter},Ptr{UInt8},Ptr{UInt8},Ptr{Void},Ptr{AVFilterGraph}),filt_ctx,filt,name,args,opaque,graph_ctx)
end

function avfilter_graph_config(graphctx,log_ctx)
    ccall((:avfilter_graph_config,libavfilter),Cint,(Ptr{AVFilterGraph},Ptr{Void}),graphctx,log_ctx)
end

function avfilter_graph_free(graph)
    ccall((:avfilter_graph_free,libavfilter),Void,(Ptr{Ptr{AVFilterGraph}},),graph)
end

function avfilter_graph_parse(graph,filters,inputs,outputs,log_ctx)
    ccall((:avfilter_graph_parse,libavfilter),Cint,(Ptr{AVFilterGraph},Ptr{UInt8},Ptr{AVFilterInOut},Ptr{AVFilterInOut},Ptr{Void}),graph,filters,inputs,outputs,log_ctx)
end

max_stars_repo_name: asafasban/AAVideoIO

max_issues_repo_name: asafasban/AAVideoIO

max_forks_repo_name: asafasban/AAVideoIO


text:
module day04
using Pipe: @pipe

include("inputs.jl")
ex = inputs.ex_d4()
invalid = inputs.invalid_d4()
valid = inputs.valid_d4()
input = inputs.input_d4()

function count_valid_pp(input)
    passports = [Dict{String, String}(split.(blob, ":")) for blob in input]
    required = ["byr", "iyr", "eyr", "hgt", "hcl", "ecl", "pid"]
    counts = [0, 0]

    # Part 1 test
    for p in passports
        counts[1] += all([haskey(p, req) for req in required]) ? 1 : 0
    end

    # Part 2 test
    eyecolors = ["amb", "blu", "brn", "gry", "grn", "hzl", "oth"]
    for p in passports
        height = match(r"^(\d+)(\w+)$", get(p, "hgt", ""))
        counts[2] += (
            all([haskey(p, req) for req in required]) &&
            1920 <= parse(Int, get(p, "byr", "-1")) <= 2002 &&
            2010 <= parse(Int, get(p, "iyr", "-1")) <= 2020 &&
            2020 <= parse(Int, get(p, "eyr", "-1")) <= 2030 &&
            height isa RegexMatch && (
                (height[2] == "cm" && (150 <= parse(Int, height[1]) <= 193)) ||
                (height[2] == "in" && (59  <= parse(Int, height[1]) <= 76))) &&
            occursin(r"^#[a-f0-9]{6}$", get(p, "hcl", "-1")) &&
            get(p, "ecl", "-1") ∈ eyecolors &&
            occursin(r"^\d{9}$", get(p, "pid", "-1")) &&
            true ? 1 : 0)

    end

    return counts
end

@assert count_valid_pp(ex)[1] == 2
@info "Day 4, Part 1 answer: $(count_valid_pp(input)[1])"

@assert count_valid_pp(invalid)[2] == 0
@assert count_valid_pp(valid)[2] == 4
@info "Day 4, Part 2 answer: $(count_valid_pp(input)[2])"

end


text:
module DiariesController

using Genie.Renderer
using Genie.Router
using SearchLight

using Diaries
using Dates

using ViewHelper

function index()
    date = try
        year, month = parse.(Int, [@params(:year), @params(:month)])
        Date(year, month)
    catch
        Date(now())
    end
    fdom, ldom = firstdayofmonth(date), lastdayofmonth(date)
    # TODO: to Scope
    diaries_by_day = get_diaries_in_month(fdom, ldom)

    html!(:diaries, :index,
          diaries_by_day = diaries_by_day,
          dayrange       = fdom:Day(1):ldom,
          title          = Dates.format(date, "yyyy年m月"),
          current_year   = year( date),
          current_month  = month(date),
          now_year       = year(Date(now())),
          now_month      = month(Date(now())),
    )
end

function show()
    date = Date(@params(:date))
    diary = SearchLight.find_one_by(Diary, SQLWhereExpression("date = ?", date)) |> get
    html!(:diaries, :show, diary = diary, year = year(date), month = month(date))
end

function new()
    date = Date(@params(:date))
    diary = Diary(date = date)
    html!(:diaries, :new, diary = diary, year = year(date), month = month(date))
end

function create()
    date = Date(@params(:date))
    Diary(date = @params(:date), weather = @params(:weather), note = @params(:note)) |> save &&
        redirect_to(link_to(:show_diary, date = date))
end

function edit()
    date = Date(@params(:date))
    diary = SearchLight.find_one_by(Diary, SQLWhereExpression("date = ?", date)) |> get
    html!(:diaries, :edit, diary = diary, year = year(date), month = month(date))
end

function update()
    date = Date(@params(:date))
    diary = SearchLight.find_one!!(Diary, @params(:diary_id))
    diary.date    = date
    diary.weather = @params(:weather)
    diary.note    = @params(:note)
    save(diary) && redirect_to(link_to(:show_diary, date = date))
end

function get_diaries_in_month(fdom, ldom)
    diaries = SearchLight.find_by(Diary, SQLWhereExpression("date >= ? AND date <= ?", [fdom, ldom]), order = SQLOrder("date"))
    reduce((dict, d) -> merge(dict, Dict(d.date => d)), diaries, init = Dict{Date, Diary}())
end

end

max_stars_repo_name: ky-s/DiaryApp

max_issues_repo_name: ky-s/DiaryApp

max_forks_repo_name: ky-s/DiaryApp


text:
"""
	struct ExtractKeyAsField{S,V} <: AbstractExtractor
		key::S
		item::V
	end

extracts all items in `vec` and in `other` and return them as a ProductNode.
"""
struct ExtractKeyAsField{S,V} <: BagExtractor
	key::S
	item::V
end

extract_empty_bag_item(e::ExtractKeyAsField, store_input) =
	ProductNode((key = e.key(extractempty; store_input), item = e.item(extractempty; store_input)))[1:0]

(s::ExtractKeyAsField)(v::MissingOrNothing; store_input=false) = extract_missing_bag(s, v; store_input)

(e::ExtractKeyAsField)(v::ExtractEmpty; store_input=false) =
	make_empty_bag(ProductNode((key = e.key(v; store_input), item = e.item(v; store_input))), v)

function (e::ExtractKeyAsField)(vs::Dict; store_input=false)
	isempty(vs) && return extract_missing_bag(e, vs; store_input)
	items = mapreduce(catobs, collect(vs)) do (k,v)
		ProductNode((key = e.key(k; store_input), item = e.item(v; store_input)))
	end
	_make_bag_node(items, [1:length(vs)], [vs], store_input)
end

Base.hash(e::ExtractKeyAsField, h::UInt) = hash((e.key, e.item), h)
Base.:(==)(e1::ExtractKeyAsField, e2::ExtractKeyAsField) = e1.key == e2.key && e1.item == e2.item


text:
module JTAGManager

export JTAG

# external dependencies
using ProgressMeter

# stdlib dependencies
using Sockets

# Messages - for interfacing the the TCL server running in either Quartus
# SystemConsole or STP.
jtag_open(id = 0) = "jtag_open $id"
jtag_close() = "jtag_close"
jtag_write(start, data) = "jtag_write $start $(join(data, " "))"
jtag_read(start, bytes) = "jtag_read $start $bytes"

"""
    JTAG(ip = ip"127.0.0.1", port = 2540) :: JTAG

Stand-in type for controlling communication to a TCL server running in Quartus
SystemConsole. Provides remote access to the JTAG read and write mechanisms.

API: `open`, `close`, `write`, `read`
"""
mutable struct JTAG{T <: IPAddr}
    ip :: T
    port :: Int
    # Remember if the JTAG Connection has been opened yet.
    # Default to false on construction.
    isopen :: Bool
end

JTAG(ip::IPAddr = ip"127.0.0.1", port::Int = 2540) = JTAG(ip, port, false)

Sockets.connect(jtag::JTAG) = connect(jtag.ip, jtag.port)

"""
    open(jtag::JTAG)

Open the SystemConsole connection to an FPGA. Should not be commonly called
directly.
"""
function Base.open(jtag::JTAG)
    socket = connect(jtag)
    println(socket, jtag_open())
    close(socket)
    # Set open flag.
    jtag.isopen = true
    return nothing
end
safeopen(jtag::JTAG) = (jtag.isopen || open(jtag); return nothing)

"""
    close(jtag::JTAG)

Close the SystemConsole connection to an FPGA.
"""
function Base.close(jtag::JTAG)
    socket = connect(jtag)
    println(socket, jtag_close())
    close(socket)
    # Clear open flag.
    jtag.isopen = false
end

"""
    write(jtag::JTAG, address, data; bytes_per_transfer = 2^18)

Write `data` to the connected device starting at `addresses`. If `data` is a
large vector, break it up into sizes of `bytes_per_transfer`.
"""
Base.write(jtag::JTAG, address, data::Integer; kwargs...) = write(jtag, address, [data]; kwargs...)
function Base.write(
            jtag::JTAG,
            address,
            data;
            bytes_per_transfer = 2 ^ 18,
       )

    # Establish connection if it hasn't been established yet.
    safeopen(jtag)

    # Open the socket to the TCL server
    socket = connect(jtag)

    # Compute the number of transfers
    ntransfers = ceil(Int, length(data) / bytes_per_transfer)

    # Transfer data in chunks
    @showprogress for i in 1:ntransfers
        # Compute start address for this transfer
        start = address + bytes_per_transfer * (i-1)

        # Get the index range of data for this transfer
        idx_low = bytes_per_transfer * (i-1) + 1
        idx_high = min(bytes_per_transfer * i, length(data))

        transfer_data = view(data, idx_low:idx_high)

        # Send the data over TCP
        println(socket, jtag_write(start, transfer_data))
    end
    close(socket)
end

"""
    read(jtag::JTAG, address, bytes; bytes_per_transfer = 2^18) :: Vector{UInt8}

Read `bytes` number of bytes from the JTAG device starting at `address`.
"""
function Base.read(jtag::JTAG, address, bytes; bytes_per_transfer = 2^18)

    # Do the standard setup
    safeopen(jtag)
    socket = connect(jtag)
    ntransfers = ceil(Int, length(bytes) / bytes_per_transfer)

    # Data to store result.
    data = UInt8[]
    sizehint!(data, bytes)

    @showprogress for i in 1:ntransfers
        start = address + bytes_per_transfer * (i-1)

        println(socket, jtag_read(start, bytes))
        back_string = readline(socket)

        # Split the returned string on spaces, parse back as UInt8
        split_string = split(back_string, " ")
        append!(data, parse.(UInt8, split_string))
    end

    close(socket)

    return data
end

end # module


text:
module PSPipePlanckRender

end

max_stars_repo_name: xzackli/PSPipePlanckRender

max_issues_repo_name: xzackli/PSPipePlanckRender

max_forks_repo_name: xzackli/PSPipePlanckRender


text:
abstract type LSHIndex end

function _get_partition(threshold, max_n_hashes)
    min_error = Inf
    opt = (0, 0)
    bestfp, bestfn = (0,0)
    for b in 1 : max_n_hashes
        max_r = Int(floor(max_n_hashes / b))
        for r in 1 : max_r
            fp = _integrate(s -> 1 - (1 - s^r)^b, 0.0, threshold)
            fn = _integrate(s -> (1 - s^r)^b, threshold, 1.0)
            error = fp + fn
            if error < min_error
                min_error = error
                opt = (b, r)
                bestfp, bestfn = fp, fn
            end
        end
    end
    @show bestfp, bestfn
    return opt
end

function _integrate(f, a, b)
    p = 0.001
    area = 0.0
    x = a
    while x < b
        area = area + f(x+0.5*p)*p
        x = x + p
    end
    return area
end

"""
Find all pairs of similar strings in index.
"""
function similar_pairs(mh::LSHIndex)::Vector{Tuple{Int, Int}}
    sims = Set{Tuple{Int, Int}}()
    for table in mh.tables
        for similars in values(table)
            for i in 1:length(similars), j in i + 1 : length(similars)
                push!(sims, (similars[i],similars[j]))
            end
        end
    end
    collect(sims)
end


text:
abstract type AbstractEnergyModel end

abstract type AbstractData end

abstract type Expression end

"Can be `add`ed to a JuMP model, must have an `objects` dictionary and a `class`"
abstract type Component end

abstract type DeviceFormulation end

"Connected to at least one `Bus`. Additionally to `addto!` provides `p`, `cost` and `busattributes`"
abstract type Device{DF<:DeviceFormulation} <: Component end

abstract type ModelType end

# Expressions

abstract type Emission <: Expression end
struct Cost <: Expression end
struct CO2 <: Emission end
struct NodalActivePower <: Expression end
struct NodalReactivePower <: Expression end

max_stars_repo_name: coroa/EnergyModels

max_issues_repo_name: coroa/EnergyModels

max_forks_repo_name: coroa/EnergyModels


text:
module TabularTMazeHeatmapSet

using Curiosity
using MinimalRLCore
const TTMU = Curiosity.TabularTMazeUtils
include("../experiments/tabular_tmaze.jl")
using Random
using StatsBase
using GVFHordes
using Statistics
using JLD2

StatsBase.sample(p::GVFHordes.GVFParamFuncs.FunctionalPolicy, s, actions) =
    sample(Weights([p.func(s, a) for a in actions]))

function gen_dataset()
    parsed = TabularTMazeExperiment.default_args()
    parsed["cumulant_schedule"] = "Constant"
    parsed["cumulant"] = 1.0

    cumulant_schedule = TTMU.get_cumulant_schedule(parsed)

    exploring_starts = parsed["exploring_starts"]
    env = TabularTMaze(exploring_starts, cumulant_schedule)
    start_states = Curiosity.valid_state_mask()

    function del!(arr, ind_to_del)
        arr = arr[1:end .!= ind_to_del]
        return arr
    end
    #NOTE: This lets the agent start in terminal states. This is not normally possible.
    # Work backwards with index so you dont have to reindex due to shifting
    ind_to_delete = [21, 17, 5, 1]
    for i in ind_to_delete
        start_states = del!(start_states, i)
    end
    observations = []
    for s in start_states
        MinimalRLCore.reset!(env,s)
        push!(observations, MinimalRLCore.get_state(env))
    end

    num_actions = 4

    total_obs = repeat(observations, inner=num_actions)
    total_actions = repeat(collect(1:4),outer=length(observations))

    return total_obs, total_actions
end

function save_data(obs,actions)
    ValueSet = Dict()
    ValueSet["actions"] = actions
    ValueSet["states"] = obs
    @save "./src/data/TTMazeValueSet.jld2" ValueSet
end

end

max_stars_repo_name: MatthewMcLeod/curiosity

max_issues_repo_name: MatthewMcLeod/curiosity

max_forks_repo_name: MatthewMcLeod/curiosity


text:
using ApplicationBuilder
using Pkg

examples_blink = joinpath(@__DIR__, "..", "..", "examples", "sdl.jl")

# Allow this file to be called either as a standalone file to build the above
# example, or from runtests.jl using a globally-defined builddir.
@isdefined(builddir) || (builddir="builddir")

using SimpleDirectMediaLayer
SDL2 = SimpleDirectMediaLayer

sdlPkg = Pkg.dir("SimpleDirectMediaLayer")
homebrew = Pkg.dir("Homebrew")

# Copy and modify all of the required libraries for SDL2 to run. Only need
# to do this manual step once.
# Note: this step is required because the SDL libs are linked to reference
# eachother, so some manual modification must be done. Perhaps this could be
# automated by ApplicationBuilder as well in the future.
libs = joinpath(builddir, "sdl_libs")
mkpath(libs)
function cp_lib(l)
    name = basename(l)
    cp(l, joinpath(libs, name), follow_symlinks=true, force=true)
    l = joinpath(libs, name)
    run(`install_name_tool -id "$name" $l`)
    try
        external_deps = readlines(pipeline(`otool -L $l`,
                    `grep .julia`,  # filter julia lib deps
                    `sed 's/(.*)$//'`)) # remove trailing parens
        for line in external_deps
            path = strip(line)
            depname = basename(path)
            depname = "$(match(r"(\w+)", depname)[1]).dylib" # strip version
            cmd = `install_name_tool -change "$path" "@rpath/$depname" $l`
            println(cmd)
            run(cmd)
        end
    catch
    end
end
cp_lib(SDL2.libSDL2)
cp_lib(SDL2.libSDL2_ttf)
cp_lib(SDL2.libSDL2_mixer)
cp_lib(joinpath(homebrew, "deps/usr/lib/libmodplug.dylib"))
cp_lib(joinpath(homebrew, "deps/usr/lib/libvorbisfile.dylib"))
cp_lib(joinpath(homebrew, "deps/usr/lib/libvorbis.dylib"))
cp_lib(joinpath(homebrew, "deps/usr/lib/libfreetype.dylib"))
cp_lib(joinpath(homebrew, "deps/usr/lib/libogg.dylib"))
cp_lib(joinpath(homebrew, "deps/usr/lib/libpng16.dylib"))

build_app_bundle(examples_blink;
    verbose = true,
    resources = [joinpath(sdlPkg,
                         "assets","fonts","FiraCode","ttf","FiraCode-Regular.ttf"),
                ],
    libraries = ["$libs/*"],
    appname="HelloSDL2", builddir=builddir)


text:
using SchoolBusRouting, JLD, DataFrames, CSV

# experiment parameters
conditions = ["one", "many", "LBH"]
maxiter = 1
maxtime = 3600.
fractions = [0.05, 0.1, 0.15, 0.2, 0.25]
schools = [50, 100, 150, 200]

# output arrays
exp_id = Int[]; exp_name = AbstractString[]; nbuses = Int[]; iter = Int[];

counter = 0
for i = 1:maxiter, f in fractions, s in schools, condition in conditions
	srand(i)
	counter += 1
	data = SchoolBusRouting.syntheticproblem(nschools=s,
	                            nyards=1,
	                            districtsize = 30000.,
	                            dwelltime = 10.*60.,
	                            studentsperschool = 100,
	                            schoolrange=(3000., 8000.),
	                            studentspread=200.,
	                            studentclustering=0.66,
	                            maxwalking=1610*0.5,
	                            d2dfraction=f)
	println("***********************")
	@show data
	data.params.max_time_on_bus = maxtime
	data.params.velocity = 35/3.6 # 35 kph
	scenariolist = []
	params = []
	if condition == "many"
		# many scenarios
		λs = [1e2, 5e2, 5e3, 1e4, 5e4, 1e5];
		params = [SchoolBusRouting.ScenarioParameters(maxRouteTimeLower=maxtime-2000,
										 maxRouteTimeUpper=maxtime+2000,
										 nGreedy=10, λ=λ,
										 nIterations=10) for λ in λs]
		srand(i)
		scenariolist = SchoolBusRouting.computescenarios(data, params, OutputFlag=0, MIPGap=1e-4, TimeLimit=90);
	elseif condition == "one"
		params = [SchoolBusRouting.ScenarioParameters(maxRouteTimeLower=maxtime-2000,
										 maxRouteTimeUpper=maxtime+2000,
										 nGreedy=20, λ=λ,
										 nIterations=50) for λ in [1e8]]
		srand(i)
		scenariolist = SchoolBusRouting.computescenarios(data, params, OutputFlag=0, MIPGap=1e-4, TimeLimit=90);
	end
	if condition == "LBH"
		scenarios, routes, usedscenario, buses = SchoolBusRouting.lbh(data)
		SchoolBusRouting.loadsolution!(data, scenarios, routes, usedscenario, buses)
	else
		SchoolBusRouting.loadroutingscenarios!(data, scenariolist);
		SchoolBusRouting.routeBuses!(data, OutputFlag=0, TimeLimit=3600);
	end
	SchoolBusRouting.test(data)
	println("$counter - $condition: $(length(data.buses)) buses")
	CSV.write("$(Pkg.dir("SchoolBusRouting"))/output/tableS2/solution-$counter.csv", SchoolBusRouting.output(data))
	push!(nbuses, length(data.buses))
	push!(exp_id, counter)
	push!(exp_name, condition)
	push!(iter, i)
	JLD.save("$(Pkg.dir("SchoolBusRouting"))/output/tableS2/solution-$counter.jld", "data", data)
end
df = DataFrame()
df[:exp_id] = exp_id; df[:exp_name] = exp_name; df[:buses] = nbuses; df[:iter] = iter;
CSV.write("$(Pkg.dir("SchoolBusRouting"))/output/pnas-tableS2.csv", df)

max_stars_repo_name: jihoyeo/SchoolBusRouting

max_issues_repo_name: jihoyeo/SchoolBusRouting

max_forks_repo_name: jihoyeo/SchoolBusRouting


text:
module test_octo_repo

using Test
using Octo.Adapters.SQL

Repo.connect(adapter=SQL)

Repo.query("select * from albums")
Repo.query(SubString("select * from albums"))
Repo.execute("select * from albums")
Repo.execute(SubString("select * from albums"))

end # module test_octo_repo


text:
using PastaQ
using ITensors
using Test
using Random

@testset "trotter circuit: real dynamics" begin
  N = 2
  B = 0.2

  sites = qubits(N)

  ampo = OpSum()
  # loop over the pauli operators
  for j in 1:N-1
    ampo += -1.0,"Z",j,"Z",j+1
    ampo += -B,"X",j
  end
  ampo += -B,"X",N
  Hmpo = MPO(ampo, sites)

  H = OpSum()
  # loop over the pauli operators
  for j in 1:N-1
    H += -1.0,"ZZ",(j,j+1)
    H += -B,"X",j
  end
  H += -B,"X",N

  t = 1.0
  ψ₀ = productstate(sites)
  ψtest = noprime(exp(-im * prod(Hmpo) * t) * prod(ψ₀))

  # 1. Specify final time and trotter step
  δt = 0.001
  t = 1.0
  circuit = trottercircuit(H; δt = δt, t = t)
  @test length(circuit) == 1000
  ψ = runcircuit(sites, circuit)
  @test prod(ψ) ≈ ψtest atol = 1e-5

  # 2. Specify  a time sequence
  ts = 0:0.001:1.0
  circuit = trottercircuit(H; ts = ts)
  @test length(circuit) == 1000
  ψ = runcircuit(sites, circuit)
  @test prod(ψ) ≈ ψtest atol = 1e-5
  ts = collect(0:0.001:1.0)
  circuit = trottercircuit(H; ts = ts)
  @test length(circuit) == 1000
  ψ = runcircuit(sites, circuit)
  @test prod(ψ) ≈ ψtest atol = 1e-5

  # 3. Specify a non-uniform time-sequence
  ts = vcat(collect(0:0.001:0.5), collect(0.502:0.002:1.0))
  circuit = trottercircuit(H; ts = ts)
  @test length(circuit) == length(ts)-1
  ψ = runcircuit(sites, circuit)
  @test prod(ψ) ≈ ψtest atol = 1e-5
end

@testset "trotter circuit: imaginary dynamics" begin
  N = 2
  B = 0.2

  sites = qubits(N)

  ampo = OpSum()
  # loop over the pauli operators
  for j in 1:N-1
    ampo += -1.0,"Z",j,"Z",j+1
    ampo += -B,"X",j
  end
  ampo += -B,"X",N
  Hmpo = MPO(ampo, sites)

  H = OpSum()
  # loop over the pauli operators
  for j in 1:N-1
    H += -1.0,"ZZ",(j,j+1)
    H += -B,"X",j
  end
  H += -B,"X",N

  τ = 1.0
  ψ₀ = productstate(sites)
  ψtest = noprime(exp(-prod(Hmpo) * τ) * prod(ψ₀))

  ## 1. Specify final time and trotter step
  δτ = 0.001
  circuit = trottercircuit(H; δτ = δτ, τ = τ)
  @test length(circuit) == 1000
  ψ = runcircuit(sites, circuit)
  @test prod(ψ) ≈ ψtest atol = 1e-5

  # 2. Specify  a time sequence
  τs = 0:0.001:1.0
  circuit = trottercircuit(H; τs = τs)
  @test length(circuit) == 1000
  ψ = runcircuit(sites, circuit)
  @test prod(ψ) ≈ ψtest atol = 1e-5
  τs = collect(0:0.001:1.0)
  circuit = trottercircuit(H; τs = τs)
  @test length(circuit) == 1000
  ψ = runcircuit(sites, circuit)
  @test prod(ψ) ≈ ψtest atol = 1e-5

  # 3. Specify a non-uniform time-sequence
  τs = vcat(collect(0:0.001:0.5), collect(0.502:0.002:1.0))
  circuit = trottercircuit(H; τs = τs)
  @test length(circuit) == length(τs)-1
  ψ = runcircuit(sites, circuit)
  @test prod(ψ) ≈ ψtest atol = 1e-5
end

#@testset "trotter layer 1" begin
#  Random.seed!(1234)
#  N = 20
#  H = OpSum()
#  nterms = 30
#  for j in 1:nterms
#    h = rand()
#    localop = rand(["X","Y","Z"])
#    site = rand(1:N)
#    H += h, localop, site
#
#    g = rand()
#    localop = rand(["CX","CZ"])
#    siteA = rand(1:N)
#    siteB = rand(deleteat!(collect(1:N),siteA))
#    H += g, localop, siteA,siteB
#  end
#
#  layer = PastaQ.trotterlayer(H, 0.1)
#  twoqubitgates = layer[1:nterms]
#  @test length(layer) == 2*nterms
#  @test all([g[2] isa Tuple for g in twoqubitgates])
#  onequbitgates = layer[nterms+1:2*nterms]
#  @test all([g[2] isa Int for g in onequbitgates])
#end
#
#@testset "trotter layer 2" begin
#  Random.seed!(1234)
#  N = 20
#  H = OpSum()
#  nterms = 30
#  for j in 1:nterms
#    h = rand()
#    localop = rand(["X","Y","Z"])
#    site = rand(1:N)
#    H += h, localop, site
#
#    g = rand()
#    localop = rand(["CX","CZ"])
#    siteA = rand(1:N)
#    siteB = rand(deleteat!(collect(1:N),siteA))
#    H += g, localop, siteA,siteB
#  end
#  layer = PastaQ.trotterlayer(H; δt = 0.1)
#  twoqubitgates = layer[1:nterms]
#  for (j,g) in enumerate(twoqubitgates)
#    @test g == layer[end+1-j]
#    @test g[2] isa Tuple
#  end
#  onequbitgates = layer[nterms+1:3*nterms]
#  for g in onequbitgates
#    @test g[2] isa Int
#  end
#end


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/Materials.jl/blob/master/LICENSE

using FEMMaterials, Tensors, Test

analysis, problem, element, bc_elements, ip = get_one_element_material_analysis(:Chaboche)
update!(element, "E", 200.0e3)
update!(element, "nu", 0.3)
update!(element, "R0", 100.0)
update!(element, "Kn", 10.0)
update!(element, "nn", 20.0)
update!(element, "C1", 0.0)
update!(element, "D1", 100.0)
update!(element, "C2", 0.0)
update!(element, "D2", 1000.0)
update!(element, "Q", 0.0)
update!(element, "b", 0.1)

times = [0.0]
loads = [0.0]
dt = 0.5
E = 200.0e3
nu = 0.3
G = 0.5*E/(1+nu)
syield = 100.0
# vm = sqrt(3)*G*ga | ea = ga
ea = 2*syield/(sqrt(3)*G)
# Go to elastic border
push!(times, times[end]+dt)
push!(loads, loads[end] + ea*dt)
 # Proceed to plastic flow
push!(times, times[end]+dt)
push!(loads, loads[end] + ea*dt)
 # Reverse direction
push!(times, times[end]+dt)
push!(loads, loads[end] - ea*dt)
 # Continue and pass yield criterion
push!(times, times[end]+dt)
push!(loads, loads[end] - 2*ea*dt)
loading = ShearStrainLoading(times, loads)
update_bc_elements!(bc_elements, loading)

analysis.properties.t1 = maximum(times)
run!(analysis)
stresses = zeros(length(times), 6)
for (i,t) in zip(1:length(times), times)
    stresses[i,:] = tovoigt(ip("stress", t))
    @test isapprox(stresses[i,[1,2,3,4,6]], zeros(5); atol=1e-6)
end
s31 = stresses[:,5]
eeqs = [ip("cumeq", t) for t in times]

@test isapprox(s31[2], syield/sqrt(3.0))
@test isapprox(s31[3]*sqrt(3.0), syield + 10.0*((eeqs[3]-eeqs[2])/dt)^(1.0/20.0); rtol=1e-2)
@test isapprox(s31[4], s31[3]-G*ea*dt)
@test isapprox(s31[5]*sqrt(3.0), -(syield + 10.0*((eeqs[5]-eeqs[4])/dt)^(1.0/20.0)); rtol=1e-2)


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

module LoadBalancerClient

using Random
using Dates
using Swagger
import Swagger: field_name, property_type, hasproperty, validate_property, SwaggerApi, SwaggerModel
import Base: convert, propertynames

include("modelincludes.jl")

include("api_LoadBalancersApi.jl")

# export models
export convert, ApplicationGatewayBackendAddress
export convert, ApplicationGatewayBackendAddressPool
export convert, ApplicationGatewayBackendAddressPoolPropertiesFormat
export convert, ApplicationSecurityGroup
export convert, ApplicationSecurityGroupPropertiesFormat
export convert, BackendAddressPool
export convert, BackendAddressPool2
export convert, BackendAddressPoolPropertiesFormat
export convert, BackendAddressPoolPropertiesFormat2
export convert, CloudError
export convert, CloudErrorBody
export convert, CustomDnsConfigPropertiesFormat
export convert, DdosSettings
export convert, Delegation
export convert, FlowLog
export convert, FlowLogFormatParameters
export convert, FlowLogPropertiesFormat
export convert, FrontendIPConfiguration
export convert, FrontendIPConfiguration2
export convert, FrontendIPConfigurationPropertiesFormat
export convert, FrontendIPConfigurationPropertiesFormat2
export convert, IPAllocationMethod
export convert, IPConfiguration
export convert, IPConfigurationProfile
export convert, IPConfigurationProfilePropertiesFormat
export convert, IPConfigurationPropertiesFormat
export convert, IPVersion
export convert, InboundNatPool
export convert, InboundNatPoolPropertiesFormat
export convert, InboundNatRule
export convert, InboundNatRule2
export convert, InboundNatRuleListResult
export convert, InboundNatRulePropertiesFormat
export convert, InboundNatRulePropertiesFormat2
export convert, IpTag
export convert, LoadBalancer
export convert, LoadBalancerBackendAddress
export convert, LoadBalancerBackendAddress2
export convert, LoadBalancerBackendAddressPoolListResult
export convert, LoadBalancerBackendAddressPropertiesFormat
export convert, LoadBalancerBackendAddressPropertiesFormat2
export convert, LoadBalancerFrontendIPConfigurationListResult
export convert, LoadBalancerListResult
export convert, LoadBalancerLoadBalancingRuleListResult
export convert, LoadBalancerOutboundRuleListResult
export convert, LoadBalancerProbeListResult
export convert, LoadBalancerPropertiesFormat
export convert, LoadBalancerSku
export convert, LoadBalancingRule
export convert, LoadBalancingRulePropertiesFormat
export convert, NetworkInterface
export convert, NetworkInterfaceDnsSettings
export convert, NetworkInterfaceIPConfiguration
export convert, NetworkInterfaceIPConfigurationPrivateLinkConnectionProperties
export convert, NetworkInterfaceIPConfigurationPropertiesFormat
export convert, NetworkInterfaceListResult
export convert, NetworkInterfacePropertiesFormat
export convert, NetworkInterfaceTapConfiguration
export convert, NetworkInterfaceTapConfigurationPropertiesFormat
export convert, NetworkSecurityGroup
export convert, NetworkSecurityGroupPropertiesFormat
export convert, OutboundRule
export convert, OutboundRulePropertiesFormat
export convert, PrivateEndpoint
export convert, PrivateEndpointProperties
export convert, PrivateLinkServiceConnection
export convert, PrivateLinkServiceConnectionProperties
export convert, PrivateLinkServiceConnectionState
export convert, Probe
export convert, ProbePropertiesFormat
export convert, ProvisioningState
export convert, PublicIPAddress
export convert, PublicIPAddressDnsSettings
export convert, PublicIPAddressPropertiesFormat
export convert, PublicIPAddressSku
export convert, Resource
export convert, ResourceNavigationLink
export convert, ResourceNavigationLinkFormat
export convert, RetentionPolicyParameters
export convert, Route
export convert, RouteNextHopType
export convert, RoutePropertiesFormat
export convert, RouteTable
export convert, RouteTablePropertiesFormat
export convert, SecurityRule
export convert, SecurityRuleAccess
export convert, SecurityRuleDirection
export convert, SecurityRulePropertiesFormat
export convert, ServiceAssociationLink
export convert, ServiceAssociationLinkPropertiesFormat
export convert, ServiceDelegationPropertiesFormat
export convert, ServiceEndpointPolicy
export convert, ServiceEndpointPolicyDefinition
export convert, ServiceEndpointPolicyDefinitionPropertiesFormat
export convert, ServiceEndpointPolicyPropertiesFormat
export convert, ServiceEndpointPropertiesFormat
export convert, SubResource
export convert, Subnet
export convert, SubnetPropertiesFormat
export convert, TagsObject
export convert, TrafficAnalyticsConfigurationProperties
export convert, TrafficAnalyticsProperties
export convert, TransportProtocol
export convert, TransportProtocol2
export convert, VirtualNetworkTap
export convert, VirtualNetworkTapPropertiesFormat

# export operations
export convert, LoadBalancersApi

export check_required, field_name, property_type, hasproperty, propertynames, validate_property, convert

end


text:
function UnivariatePointValuePair(arg0::jdouble, arg1::jdouble)
    return UnivariatePointValuePair((jdouble, jdouble), arg0, arg1)
end

function get_point(obj::UnivariatePointValuePair)
    return jcall(obj, "getPoint", jdouble, ())
end

function get_value(obj::UnivariatePointValuePair)
    return jcall(obj, "getValue", jdouble, ())
end