text:
################################################################################
# Pooling in image dimension (width and height)
################################################################################
function max_pooling_forward{T}(input::Array{T}, output::Array{T}, mask::Array{Csize_t}, layer)
  width, height, channels, num = size(input)
  pooled_width = size(output, 1)
  pooled_height = size(output, 2)

  for n = 1:num
    for c = 1:channels
      for ph = 1:pooled_height
        for pw = 1:pooled_width
          hstart = (ph-1)*layer.stride[2] - layer.pad[2] + 1
          wstart = (pw-1)*layer.stride[1] - layer.pad[1] + 1
          hend = min(hstart + layer.kernel[2] - 1, height)
          wend = min(wstart + layer.kernel[1] - 1, width)
          hstart = max(hstart, 1)
          wstart = max(wstart, 1)

          maxval = -Inf
          maxw = 0
          maxh = 0
          for w = wstart:wend
            for h = hstart:hend
              @inbounds val = input[w,h,c,n]
              if val > maxval
                maxval = val
                maxw = w
                maxh = h
              end
            end
          end
          @inbounds output[pw,ph,c,n] = maxval
          @inbounds mask[pw,ph,c,n] = (maxh-1) * width + maxw-1
        end
      end
    end
  end
end

function mean_pooling_forward{T}(input::Array{T}, output::Array{T}, layer)
  width, height, channels, num = size(input)
  pooled_width = size(output, 1)
  pooled_height = size(output, 2)
  kernel_size = layer.kernel[1] * layer.kernel[2]

  for n = 1:num
    for c = 1:channels
      for ph = 1:pooled_height
        for pw = 1:pooled_width
          hstart = (ph-1)*layer.stride[2] - layer.pad[2] + 1
          wstart = (pw-1)*layer.stride[1] - layer.pad[1] + 1
          hend = min(hstart + layer.kernel[2] - 1, height)
          wend = min(wstart + layer.kernel[1] - 1, width)
          hstart = max(hstart, 1)
          wstart = max(wstart, 1)

          the_sum = 0.0
          for w = wstart:wend
            for h = hstart:hend
              @inbounds the_sum += input[w,h,c,n]
            end
          end
          @inbounds output[pw,ph,c,n] = the_sum / kernel_size
        end
      end
    end
  end
end

function max_pooling_backward{T}(input::Array{T}, output::Array{T}, mask::Array{Csize_t}, layer)
  width, height, channels, num = size(input)
  pooled_width = size(output, 1)
  pooled_height = size(output, 2)

  fill!(input, 0)
  for n = 1:num
    for c = 1:channels
      for ph = 1:pooled_height
        for pw = 1:pooled_width
          index = mask[pw,ph,c,n]
          idx_w = (index % width) + 1
          idx_h = div(index, width) + 1
          @inbounds input[idx_w, idx_h, c, n] += output[pw,ph,c,n]
        end
      end
    end
  end
end

function mean_pooling_backward{T}(input::Array{T}, output::Array{T}, layer)
  width, height, channels, num = size(input)
  pooled_width = size(output, 1)
  pooled_height = size(output, 2)
  kernel_size = layer.kernel[1] * layer.kernel[2]

  fill!(input, 0)
  for n = 1:num
    for c = 1:channels
      for ph = 1:pooled_height
        for pw = 1:pooled_width
          hstart = (ph-1)*layer.stride[2] - layer.pad[2] + 1
          wstart = (pw-1)*layer.stride[1] - layer.pad[1] + 1
          hend = min(hstart + layer.kernel[2] - 1, height)
          wend = min(wstart + layer.kernel[1] - 1, width)
          hstart = max(1, hstart)
          wstart = max(1, wstart)

          @inbounds val = output[pw,ph,c,n] / kernel_size
          for w = wstart:wend
            for h = hstart:hend
              @inbounds input[w,h,c,n] += val
            end
          end
        end
      end
    end
  end
end


text:
module Day2

using DelimitedFiles

mutable struct Position
    x::Int
    depth::Int
    aim::Int
end

update_aiming!(position, command) = begin
    parts = split(command)
    dir = parts[1]
    n = parse(Int, parts[2])
    if dir == "forward"
        position.x += n
        position.depth += n * position.aim
    elseif dir == "down"
        position.aim += n
    elseif dir == "up"
        position.aim -= n
    end
end

update!(position, command) = begin
    parts = split(command)
    dir = parts[1]
    n = parse(Int, parts[2])
    if dir == "forward"
        position.x += n
    elseif dir == "down"
        position.depth += n
    elseif dir == "up"
        position.depth -= n
    end
end

run_part1(commands) = begin
    pos = Position(0, 0, 0)
    for command in commands
        update!(pos, command)
    end
    @show pos.x
    @show pos.depth
    return pos.x * pos.depth
end

run_part2(commands) = begin
    pos = Position(0, 0, 0)
    for command in commands
        update_aiming!(pos, command)
    end
    return pos.x * pos.depth
end

test = ["forward 5"
        "down 5"
        "forward 8"
        "up 3"
        "down 8"
        "forward 2"];

@assert run_part1(test) == 150;
@assert run_part2(test) == 900;

input = readdlm("resources/day2_input.txt", '\t', String, '\n') |> vec;

@show (input |> run_part1)
@show (input |> run_part2)

end

ext: jl


text:
"""
"""
immutable Cursor
    graph::Node
end


text:
function ssfm_solve(num::Numerics, coeffs::Coefficients)
  time_steps = Int(floor(num.T / num.dt))
  space_steps = Int(floor(num.S / num.ds))
  time = LinRange(0, num.T, time_steps)
  space = LinRange(-num.S / 2, num.S / 2, space_steps)

  # print("Simulation diagnostics:\n")
  # @printf("Time  : from %6.3i ms, to  %6.3i ms; %6d steps \n", time[1] * 1e3, time[time_steps] * 1e3, time_steps)
  # @printf("Space : from %6.3f mm, to  %6.3f mm; %6d steps \n", space[1] * 1e3, space[space_steps] * 1e3, space_steps)

  # Spatial frequency range computation
  k = 2 * pi * LinRange(-1 / (2 * num.ds), 1 / (2 * num.ds), space_steps)
  k = fftshift(k)
  ψ = zeros(ComplexF64, space_steps, time_steps)
  ψ_spect = zeros(ComplexF64, space_steps, time_steps)

  waveform = coeffs.initial.(space)
  # fig0 = plot(space * 1e3,
  #   abs.(waveform) .^ 2,
  #   title="initial configuration",
  #   xlabel="space [mm]",
  #   show=true,
  #   reuse=false,
  #   label="initial state")

  # fig1 = plot(space * 1e3,
  #   abs.(coeffs.β.(space)),
  #   show=true,
  #   reuse=false,
  #   label="potential")

  ## [SPACE INDEX, TIME INDEX]
  ψ[:, 1] = waveform

  fwd_disp = exp.(num.dt / 2 .* coeffs.α * k .^ 2)
  fwd_curvature = exp.(num.dt / 2 .* coeffs.β.(space))
  # display("Gamma at peak")
  # display(coeffs.γ.(ψ[Int(floor(space_steps/2)),1]))
  #@showprogress "Propagating the field... "
  for n = 1:time_steps-1
    # if 1+2*app.as*(app.N - 1) * abs.(ψ[Int(floor(space_steps/2)), n])^2 > 1
    #   display("Collapse detected")
    #   return
    # end
    ψ_spect[:, n] = fft(ψ[:, n])
    ψ_spect[:, n] = ψ_spect[:, n] .* fwd_disp
    ψ[:, n] = ifft(ψ_spect[:, n])
    ψ[:, n+1] = ψ[:, n] .* exp.(num.dt / 2 .* coeffs.γ.(ψ[:, n])) .* fwd_curvature  ## this is an Euler step
  end
  ψ_spect[:, time_steps] = fft(ψ[:, time_steps])

  return time, space, ψ, ψ_spect
end

function ssfm_propagate(num::Numerics, coeffs::Coefficients)
  time_steps = Int(floor(num.T / num.dt))
  space_steps = Int(floor(num.S / num.ds))
  time = LinRange(0, num.T, time_steps)
  space = LinRange(-num.S / 2, num.S / 2, space_steps)

  # Spatial frequency range computation
  k = 2 * pi * LinRange(-1 / (2 * num.ds), 1 / (2 * num.ds), space_steps)
  k = fftshift(k)
  ψ = zeros(ComplexF64, space_steps)
  ψ_spect = zeros(ComplexF64, space_steps)

  waveform = coeffs.initial.(space)

  ## [SPACE INDEX, TIME INDEX]
  ψ = waveform

  fwd_disp = exp.(num.dt / 2 .* coeffs.α * k .^ 2)
  fwd_curvature = exp.(num.dt / 2 .* coeffs.β.(space))
  max_amplitude = maximum(abs.(ψ).^2)
  for n = 1:time_steps-1

    ψ_spect = fft(ψ)
    ψ_spect = ψ_spect .* fwd_disp
    ψ = ifft(ψ_spect)
    ψ= ψ .* exp.(num.dt / 2 .* coeffs.γ.(ψ)) .* fwd_curvature  ## this is an Euler step
    if max_amplitude < maximum(abs.(ψ).^2)
      max_amplitude = maximum(abs.(ψ).^2)
    end
  end
  ψ_spect = fft(ψ)

  return time, space, ψ, ψ_spect, max_amplitude
end

function ground_state_solve(num::Numerics, coeffs::Coefficients)
  time_steps = Int(floor(num.T / num.dt))
  space_steps = Int(floor(num.S / num.ds))
  time = LinRange(0, num.T, time_steps)
  space = LinRange(-num.S / 2, num.S / 2, space_steps)

  # print("Simulation diagnostics:\n")
  # @printf("Time  : from %6.3i ms, to  %6.3i ms; %6d steps \n", time[1] * 1e3, time[time_steps] * 1e3, time_steps)
  # @printf("Space : from %6.3f mm, to  %6.3f mm; %6d steps \n", space[1] * 1e3, space[space_steps] * 1e3, space_steps)

  # Spatial frequency range computation
  k = 2 * pi * LinRange(-1 / (2 * num.ds), 1 / (2 * num.ds), space_steps)
  k = fftshift(k)
  #fig = Plots.plot(k, show=true)
  ψ = zeros(ComplexF64, space_steps)
  ψ_spect = zeros(ComplexF64, space_steps)

  #Implement gaussian initial state
  wave(s::Float64) = sqrt(1 / (sqrt(2 * pi) * num.S / 3)) * exp.(-(s) .^ 2 / (4 * num.S / 3^2))
  integral, error = quadgk(s -> abs(wave(s))^2, -num.S / 2, +num.S / 2)

  waveform = wave.(space) / sqrt(integral)

  ## [SPACE INDEX, TIME INDEX]
  ψ[:] = waveform
  fwd_disp = exp.(-im * num.dt / 2 .* coeffs.α * k .^ 2)
  fwd_curvature = exp.(-im * num.dt / 2 .* coeffs.β.(space))
  # display("Gamma at peak")
  # display(coeffs.γ.(ψ[Int(floor(space_steps/2)),1]))
  target_Δ = 1e-10
  Δ = 1
  let prev, current
    prev = 0
    while Δ > target_Δ
      # if 1+2*app.as*(app.N - 1) * abs.(ψ[Int(floor(space_steps/2)), n])^2 > 1
      #   display("Collapse detected")
      #   return
      # end
      ψ_spect = fft(ψ)
      ψ_spect = ψ_spect .* fwd_disp
      ψ = ifft(ψ_spect)
      ψ = ψ .* exp.(im * num.dt / 2 .* coeffs.γ.(ψ)) .* fwd_curvature  ## this is an Euler step
      # renormalization
      #isplay(sum(abs.(ψ) .^ 2 * num.ds))
      current = sum(abs.(ψ) .^ 2 * num.ds)
      Δ = abs(prev - current)
      ψ = ψ / sqrt(current)
      prev = current
    end
  end
  ψ_spect = fft(ψ)
  return time, space, ψ, ψ_spect
end

function plot_dynamics(time, space, ψ, ψ_spect)
  pyplot()
  # Natural orientation choice for plot
  fig2 = heatmap(time * 1e3,
    space * 1e3,
    abs.(ψ[:, :]) .^ 2,
    show=true,
    title="|ψ|^2",
    ylabel="space [mm]",
    xlabel="time [ms]",
    colorrange=(0, 1),
    reuse=false,
    size=(800, 600))

  # fig2_bis = heatmap(time * 1e3,
  #   space * 1e3,
  #   angle.(ψ[:, :]),
  #   show=true,
  #   title="|ψ|^2",
  #   ylabel="space [mm]",
  #   xlabel="time [ms]",
  #   colorrange=(0, 1),
  #   reuse=false,
  #   size=(800, 600))

  fig3 = plot(space * 1e3,
    abs.(ψ[:, end]) .^ 2,
    title="evolved |ψ|^2",
    xlabel="space [mm]",
    reuse=false,
    label="t=t_max")
  plot!(space * 1e3,
    abs.(ψ[:, 1]) .^ 2,
    show=true,
    reuse=false,
    label="t=0",
    size=(800, 600))
end

function plot_ground_state(time, space, ψ, ψ_spect)
  pyplot()
  # Natural orientation choice for plot
  fig3 = plot(space * 1e3,
    abs.(ψ) .^ 2,
    title="stationary |ψ|^2",
    xlabel="space [mm]",
    reuse=false,
    show=true,
    label="t=t_max")
end

max_stars_repo_name: lorenzifrancesco/soliton-BEC

max_issues_repo_name: lorenzifrancesco/soliton-BEC

max_forks_repo_name: lorenzifrancesco/soliton-BEC


text:
using Base.Meta

iscall(e) = false
iscall(e::Expr) = isexpr(e, :call)

getfunsym(e::Expr) = iscall(e) ? e.args[1] : e.head
getfunsym(e) = e

setfunsym!(e::Expr, s) = iscall(e) ? (e.args[1] = s) : (e.head = s)
setfunsym!(e, s) = s

getfunargs(e::Expr) = e.args[(iscall(e) ? 2 : 1):end]

getfunargs(e) = []

function setfunargs!(e::Expr, args::Vector)
    e.args[(iscall(e) ? 2 : 1):end] = args
end
setfunargs!(e, args) = []

istree(e::Expr) = true
istree(a) = false


text:
module LeastsquaresWrapper

using JavaCall

const AbsoluteDate = @jimport org.orekit.time.AbsoluteDate
const AbstractIntegratedPropagator = @jimport org.orekit.propagation.integration.AbstractIntegratedPropagator
const BatchLSEstimator = @jimport org.orekit.estimation.leastsquares.BatchLSEstimator
const BatchLSModel = @jimport org.orekit.estimation.leastsquares.BatchLSModel
const BatchLSODModel = @jimport org.orekit.estimation.leastsquares.BatchLSODModel
const BatchLSObserver = @jimport org.orekit.estimation.leastsquares.BatchLSObserver
const Class = @jimport java.lang.Class
const ConvergenceChecker = @jimport org.hipparchus.optim.ConvergenceChecker
const DSSTBatchLSModel = @jimport org.orekit.estimation.leastsquares.DSSTBatchLSModel
const DSSTPropagator = @jimport org.orekit.propagation.semianalytical.dsst.DSSTPropagator
const EstimatedMeasurement = @jimport org.orekit.estimation.measurements.EstimatedMeasurement
const EstimationsProvider = @jimport org.orekit.estimation.measurements.EstimationsProvider
const Incrementor = @jimport org.hipparchus.util.Incrementor
const IntegratedPropagatorBuilder = @jimport org.orekit.propagation.conversion.IntegratedPropagatorBuilder
const JString = @jimport java.lang.String
const LeastSquaresOptimizer = @jimport org.hipparchus.optim.nonlinear.vector.leastsquares.LeastSquaresOptimizer
const LeastSquaresOptimizer_Optimum = @jimport org.hipparchus.optim.nonlinear.vector.leastsquares.LeastSquaresOptimizer$Optimum
const LeastSquaresProblem_Evaluation = @jimport org.hipparchus.optim.nonlinear.vector.leastsquares.LeastSquaresProblem$Evaluation
const List = @jimport java.util.List
const Map = @jimport java.util.Map
const MeasurementHandler = @jimport org.orekit.estimation.leastsquares.MeasurementHandler
const ModelObserver = @jimport org.orekit.estimation.leastsquares.ModelObserver
const MultivariateJacobianFunction = @jimport org.hipparchus.optim.nonlinear.vector.leastsquares.MultivariateJacobianFunction
const NumericalPropagator = @jimport org.orekit.propagation.numerical.NumericalPropagator
const Object = @jimport java.lang.Object
const ObservedMeasurement = @jimport org.orekit.estimation.measurements.ObservedMeasurement
const Orbit = @jimport org.orekit.orbits.Orbit
const Pair = @jimport org.hipparchus.util.Pair
const ParameterDriversList = @jimport org.orekit.utils.ParameterDriversList
const PreCompensation = @jimport org.orekit.estimation.leastsquares.PreCompensation
const PropagationType = @jimport org.orekit.propagation.PropagationType
const RealMatrix = @jimport org.hipparchus.linear.RealMatrix
const RealVector = @jimport org.hipparchus.linear.RealVector

export BatchLSEstimator
export BatchLSModel
export BatchLSODModel
export BatchLSObserver
export DSSTBatchLSModel
export MeasurementHandler
export ModelObserver
export PreCompensation
export add_measurement
export create_propagators
export equals
export estimate
export evaluation_performed
export fetch_evaluated_measurement
export get_class
export get_date
export get_evaluations_count
export get_iterations_count
export get_last_estimations
export get_measurement
export get_measurements_parameters_drivers
export get_optimum
export get_orbital_parameters_drivers
export get_physical_covariances
export get_propagator_parameters_drivers
export get_selected_propagation_drivers_for_builder
export handle_step
export hash_code
export init
export is_forward_propagation
export model_called
export notify
export notify_all
export set_convergence_checker
export set_evaluations_counter
export set_iterations_counter
export set_max_evaluations
export set_max_iterations
export set_observer
export set_parameters_convergence_threshold
export to_string
export value
export wait

include("batch_ls_estimator.jl")
include("batch_ls_model.jl")
include("batch_ls_observer.jl")
include("batch_lsod_model.jl")
include("dsst_batch_ls_model.jl")
include("measurement_handler.jl")
include("model_observer.jl")
include("pre_compensation.jl")

end


text:
module TestMaterialize

include("preamble.jl")
using BangBang: air

@testset begin
    @test materialize!!(nothing, air.([0, 1] .+ 2))::Vector{Int} == [2, 3]
    @test materialize!!([NaN, NaN], air.([0, 1] .+ 2))::Vector{Float64} == [2, 3]
    @test materialize!!(nothing, air.((0, 1) .+ 2))::Tuple === (2, 3)
    @test materialize!!(nothing, air.(SVector(0, 1) .+ 2))::SVector === SVector(2, 3)
end

end  # module


text:
module CubedSphereDiscreteModelsTests
 using GridapGeosciences
 using Gridap
 using Test

 include("ConvergenceAnalysisTools.jl")

 """
 Performs a convergence study of the discrete cubed sphere surface
 using a Gauss quadrature of degree degree for numerical integration.
 order is the order of the polynomial space used to approximate the
 surface of the cubed sphere.

 Returns a tuple with three entries.
  [1] Discrete cubed sphere surfaces for the associated values of n.
  [2] Relative errors among the approximate surface and the exact one.
  [3] Slope of the log(h)-log(err) relative error convergence curve.
 """
 function convergence_discrete_cubed_sphere_surface(n_values,order,degree)
   approx_surfaces=Float64[]
   rel_errors=Float64[]
   exact_surface=4π
   current=1
   for n in n_values
      approx_surface=compute_discrete_cubed_sphere_surface(n,order,degree)
      append!(approx_surfaces,approx_surface)
      append!(rel_errors,abs(approx_surface-exact_surface)/exact_surface)
      current=current*10
   end
   println(n_values)
   hs=[2.0/n for n in n_values]
   return approx_surfaces,rel_errors,slope(hs,rel_errors)
 end

 function compute_discrete_cubed_sphere_surface(n,order,degree)
   model = CubedSphereDiscreteModel(n,order)
   Ω = Triangulation(model)
   dΩ = Measure(Ω,degree)
   f(x)=1
   sum(∫(f)dΩ)
 end

 _,_,s=convergence_discrete_cubed_sphere_surface(generate_n_values(2),1,0)
 @test round(s,digits=3) ≈ 1.986
 _,_,s=convergence_discrete_cubed_sphere_surface(generate_n_values(2),2,4)
 @test round(s,digits=3) ≈ 4.023

 model=CubedSphereDiscreteModel(10,1;radius=rₑ)
 ψₖ = get_cell_map(model)
 @test norm(evaluate(ψₖ[rand(1:num_cells(model))],Point(1,1))) ≈ rₑ

 model=CubedSphereDiscreteModel(10;radius=rₑ)
 ψₖ = get_cell_map(model)
 @test norm(evaluate(ψₖ[rand(1:num_cells(model))],Point(rand(),rand()))) ≈ rₑ

 function detJtJ(Jt::Gridap.Fields.MultiValue{Tuple{D1,D2}}) where {D1,D2}
  J = transpose(Jt)
  sqrt(det(Jt⋅J))
 end

 function JtJ(Jt::Gridap.Fields.MultiValue{Tuple{D1,D2}}) where {D1,D2}
  J = transpose(Jt)
  Jt⋅J
 end

 function detJJt(Jt::Gridap.Fields.MultiValue{Tuple{D1,D2}}) where {D1,D2}
  J = transpose(Jt)
  sqrt(abs(det(J⋅Jt)))
 end

 function JJt(Jt::Gridap.Fields.MultiValue{Tuple{D1,D2}}) where {D1,D2}
  J = transpose(Jt)
  J⋅Jt
 end

 Ω = Triangulation(model)
 dΩ = Measure(Ω,2)

 vdetJtJ   = lazy_map(Operation(detJtJ),lazy_map(Broadcasting(∇),get_cell_map(model)))
 cfdetJtJ  = Gridap.CellData.GenericCellField(vdetJtJ,Ω,ReferenceDomain())
 p=get_cell_points(dΩ.quad)
 detJtJp=cfdetJtJ(p)

 println(detJtJp[1][1])
 println(detJtJp[1][2])
 println(detJtJp[1][3])
 println(detJtJp[1][4])

 vJtJ   = lazy_map(Operation(JtJ),lazy_map(Broadcasting(∇),get_cell_map(model)))
 cfJtJ  = Gridap.CellData.GenericCellField(vJtJ,Ω,ReferenceDomain())
 p=get_cell_points(dΩ.quad)
 JtJp=cfJtJ(p)

 println(JtJp[1][1])
 println(JtJp[1][2])
 println(JtJp[1][3])
 println(JtJp[1][4])

 vdetJJt   = lazy_map(Operation(detJJt),lazy_map(Broadcasting(∇),get_cell_map(model)))
 cfdetJJt  = Gridap.CellData.GenericCellField(vdetJJt,Ω,ReferenceDomain())
 p=get_cell_points(dΩ.quad)
 detJJtp=cfdetJJt(p)

 println(detJJtp[1][1])
 println(detJJtp[1][2])
 println(detJJtp[1][3])
 println(detJJtp[1][4])

 vJJt   = lazy_map(Operation(JJt),lazy_map(Broadcasting(∇),get_cell_map(model)))
 cfJJt  = Gridap.CellData.GenericCellField(vJJt,Ω,ReferenceDomain())
 p=get_cell_points(dΩ.quad)
 JJtp=cfJJt(p)

 println(JJtp[1][1])
 println(JJtp[1][2])
 println(JJtp[1][3])
 println(JJtp[1][4])

 vJt = lazy_map(Broadcasting(∇),get_cell_map(model))
 cfJt  = Gridap.CellData.GenericCellField(vJt,Ω,ReferenceDomain())
 p=get_cell_points(dΩ.quad)
 Jtp=cfJt(p)

 println(Jtp[1][1])
 println(Jtp[1][2])
 println(Jtp[1][3])
 println(Jtp[1][4])

end


text:
include("Tests/Tests.jl")

using Test, JuMP, GLPK
factory = with_optimizer(GLPK.Optimizer)
config = MOI.Test.TestConfig(atol=1e-5, rtol=1e-5, query=false)
@testset "Linear" begin
    Tests.linear_test(factory, config)
end


text:
abstract type AbstractReservoirDriver end

"""
    create_states(
        reservoir_driver::AbstractReservoirDriver, train_data, reservoir_matrix,input_matrix
    )

Return the trained ESN states according to the given driver.
"""
function create_states(
    reservoir_driver::AbstractReservoirDriver,
    train_data,
    washout,
    reservoir_matrix,
    input_matrix,
    bias_vector
)
    train_len = size(train_data, 2)-washout
    res_size = size(reservoir_matrix, 1)

    states = Adapt.adapt(typeof(train_data), zeros(res_size, train_len))
    tmp_array = allocate_tmp(reservoir_driver, typeof(train_data), res_size)
    _state = Adapt.adapt(typeof(train_data), zeros(res_size,1))

    for i in 1:washout
        yv = @view train_data[:,i]
        _state = next_state!(
            _state, reservoir_driver, _state, yv, reservoir_matrix,
            input_matrix, bias_vector, tmp_array
        )
    end

    for j in 1:train_len
        yv = @view train_data[:, washout+j]
        _state = next_state!(
            _state, reservoir_driver, _state, yv,
            reservoir_matrix, input_matrix, bias_vector, tmp_array
        )
        states[:,j] = _state
    end

    return states
end

function create_states(
    reservoir_driver::AbstractReservoirDriver,
    train_data,
    washout,
    reservoir_matrix::Vector,
    input_matrix,
    bias_vector
)
    train_len = size(train_data, 2)-washout
    res_size = sum([size(reservoir_matrix[i], 1) for i in 1:length(reservoir_matrix)])

    states = Adapt.adapt(typeof(train_data), zeros(res_size, train_len))
    tmp_array = allocate_tmp(reservoir_driver, typeof(train_data), res_size)
    _state = Adapt.adapt(typeof(train_data), zeros(res_size))

    for i in 1:washout
        for j in 1:length(reservoir_matrix)
            _inter_state = next_state!(
                _inter_state, reservoir_driver, _inter_state, train_data[:, i],
                reservoir_matrix, input_matrix, bias_vector, tmp_array
            )
        end
        _state = next_state!(
            _state, reservoir_driver, _state, train_data[:, i],
            reservoir_matrix, input_matrix, bias_vector, tmp_array
        )
    end

    for j in 1:train_len
        _state = next_state!(
            _state, reservoir_driver, _state, train_data[:, washout+j],
            reservoir_matrix, input_matrix, bias_vector, tmp_array
        )
        states[:,j] = _state
    end

    return states
end

#standard RNN driver
struct RNN{F,T} <: AbstractReservoirDriver
    activation_function::F
    leaky_coefficient::T
end

"""
    RNN(activation_function, leaky_coefficient)
    RNN(;activation_function=tanh, leaky_coefficient=1.0)

Returns a Recurrent Neural Network initializer for the ESN. This is the default choice.
"""
function RNN(;activation_function= NNlib.fast_act(tanh), leaky_coefficient=1.0)
    RNN(activation_function, leaky_coefficient)
end

function reservoir_driver_params(rnn::RNN, args...)
    rnn
end

function next_state!(out, rnn::RNN, x, y, W, W_in, b, tmp_array)
    mul!(tmp_array[1], W, x)
    mul!(tmp_array[2], W_in, y)
    @. tmp_array[1] = rnn.activation_function(tmp_array[1] + tmp_array[2] + b) * rnn.leaky_coefficient
    @. out = (1-rnn.leaky_coefficient)*x + tmp_array[1]
end

function next_state!(out, rnn::RNN, x, y, W::Vector, W_in, b, tmp_array)
    esn_depth = length(W)
    res_sizes = vcat(0, [get_ressize(W[i]) for i=1:esn_depth])
    inner_states = [x[1+sum(res_sizes[1:i]):sum(res_sizes[1:i+1])]
        for i in 1:esn_depth]
    inner_inputs = vcat([y], inner_states[1:end-1])

    for i in 1:esn_depth
        inner_states[i] = (1-rnn.leaky_coefficient).*inner_states[i] +
            rnn.leaky_coefficient*rnn.activation_function.((W[i]*inner_states[i]).+
            (W_in[i]*inner_inputs[i]).+reduce(vcat,b[i]))
    end
    return reduce(vcat, inner_states)
end

function allocate_tmp(::RNN, tmp_type, res_size)
    return [Adapt.adapt(tmp_type, zeros(res_size, 1)) for i in 1:2]
end

#multiple RNN driver
struct MRNN{F,T,R} <: AbstractReservoirDriver
    activation_function::F
    leaky_coefficient::T
    scaling_factor::R
end

"""
    MRNN(activation_function, leaky_coefficient, scaling_factor)
    MRNN(;activation_function=[tanh, sigmoid], leaky_coefficient=1.0,
        scaling_factor=fill(leaky_coefficient, length(activation_function)))

Returns a Multiple RNN initializer, where multiple function are combined in a linear
combination with chosen parameters ```scaling_factor```. The ```activation_function```
and ```scaling_factor``` arguments must vectors of the same size. Multiple combinations
are possible, the implementation is based upon a double activation function idea,
found in [1].

[1] Lun, Shu-Xian, et al. "_A novel model of leaky integrator echo state network for
time-series prediction._" Neurocomputing 159 (2015): 58-66.

"""
function MRNN(
    ;
    activation_function=[tanh, sigmoid],
    leaky_coefficient=1.0,
    scaling_factor=fill(leaky_coefficient, length(activation_function))
)
    @assert length(activation_function) == length(scaling_factor)
    return MRNN(activation_function, leaky_coefficient, scaling_factor)
end

function reservoir_driver_params(mrnn::MRNN, args...)
    return mrnn
end

function next_state!(out, mrnn::MRNN, x, y, W, W_in, b, tmp_array)
    @. out = (1-mrnn.leaky_coefficient)*x
    for i in 1:length(mrnn.scaling_factor)
        mul!(tmp_array[1], W, x)
        mul!(tmp_array[2], W_in, y)
        @. out += mrnn.activation_function[i](tmp_array[1] + tmp_array[2] + b) * mrnn.scaling_factor[i]
    end

    return out
end

#=
function next_state!(out, mrnn::MRNN, x, y, W, W_in, b, tmp_array)
    rnn_next_state = (1-mrnn.leaky_coefficient).*x
    for i=1:length(mrnn.scaling_factor)
        rnn_next_state += mrnn.scaling_factor[i]*mrnn.activation_function[i].((W*x).+(W_in*y).+b)
    end
    rnn_next_state
end
=#

function allocate_tmp(::MRNN, tmp_type, res_size)
    return [Adapt.adapt(tmp_type, zeros(res_size, 1)) for i in 1:2]
end

#GRU-based driver
struct GRU{F,L,R,V,B} #not an abstractreservoirdriver
    activation_function::F
    inner_layer::L
    reservoir::R
    bias::B
    variant::V
end

#https://arxiv.org/abs/1701.05923# variations of gru
"""
    FullyGated()

Returns a standard Gated Recurrent Unit ESN initializer, as described in [1].

[1] Cho, Kyunghyun, et al. “_Learning phrase representations using RNN encoder-decoder
for statistical machine translation._”
arXiv preprint arXiv:1406.1078 (2014).
"""
struct FullyGated <: AbstractGRUVariant end

"""
    Minimal()

Returns a minimal GRU ESN initializer as described in [1].

[1] Zhou, Guo-Bing, et al. "_Minimal gated unit for recurrent neural networks._"
International Journal of Automation and Computing 13.3 (2016): 226-234.
"""
struct Minimal <: AbstractGRUVariant end

#layer_init and activation_function must be vectors
"""
    GRU(;activation_function=[NNlib.sigmoid, NNlib.sigmoid, tanh],
        inner_layer = fill(DenseLayer(), 2),
        reservoir = fill(RandSparseReservoir(), 2),
        bias = fill(DenseLayer(), 2),
        variant = FullyGated())

Returns a Gated Recurrent Unit [1] reservoir driver.

[1] Cho, Kyunghyun, et al. “_Learning phrase representations using RNN encoder-decoder for
statistical machine translation._” arXiv preprint arXiv:1406.1078 (2014).
"""
function GRU(
    ;
    activation_function=[NNlib.sigmoid, NNlib.sigmoid, tanh],
    inner_layer = fill(DenseLayer(), 2),
    reservoir = fill(RandSparseReservoir(0), 2),
    bias = fill(DenseLayer(), 2),
    variant = FullyGated()
)
    return GRU(activation_function, inner_layer, reservoir, bias, variant)
end

#the actual params are only available inside ESN(), so a different driver is needed
struct GRUParams{F,V,S,I,N,SF,IF,NF} <: AbstractReservoirDriver
    activation_function::F
    variant::V
    Wz_in::S
    Wz::I
    bz::N
    Wr_in::SF
    Wr::IF
    br::NF
end

#vreation of the actual driver
function reservoir_driver_params(gru::GRU, res_size, in_size)
    gru_params = create_gru_layers(gru, gru.variant, res_size, in_size)
    return gru_params
end

#dispatch on the differenct gru variations
function create_gru_layers(gru, variant::FullyGated, res_size, in_size)

    Wz_in = create_layer(gru.inner_layer[1], res_size, in_size)
    Wz = create_reservoir(gru.reservoir[1], res_size)
    bz = create_layer(gru.bias[1], res_size, 1)

    Wr_in = create_layer(gru.inner_layer[2], res_size, in_size)
    Wr = create_reservoir(gru.reservoir[2], res_size)
    br = create_layer(gru.bias[2], res_size, 1)

    return GRUParams(gru.activation_function, variant, Wz_in, Wz, bz, Wr_in, Wr, br)
end

#check this one, not sure
function create_gru_layers(gru, variant::Minimal, res_size, in_size)

    Wz_in = create_layer(gru.inner_layer, res_size, in_size)
    Wz = create_reservoir(gru.reservoir, res_size)
    bz = create_layer(gru.bias, res_size, 1)

    Wr_in = nothing
    Wr = nothing
    br = nothing

    return GRUParams(gru.activation_function, variant, Wz_in, Wz, bz, Wr_in, Wr, br)
end

#in case the user wants to use this driver
function reservoir_driver_params(gru::GRUParams, args...)
    return gru
end

#dispatch on the important function: next_state
function next_state!(out, gru::GRUParams, x, y, W, W_in, b, tmp_array)
    gru_next_state = obtain_gru_state!(out, gru.variant, gru, x, y, W, W_in, b, tmp_array)
    return gru_next_state
end

function allocate_tmp(::GRUParams, tmp_type, res_size)
    return [Adapt.adapt(tmp_type, zeros(res_size, 1)) for i in 1:9]
end

#W=U, W_in=W in papers. x=h, and y=x. I know, it's confusing. ( on the left our notation)
#fully gated gru
function obtain_gru_state!(out, variant::FullyGated, gru, x, y, W, W_in, b, tmp_array)

    mul!(tmp_array[1], gru.Wz_in, y)
    mul!(tmp_array[2], gru.Wz, x)
    @. tmp_array[3] = gru.activation_function[1](tmp_array[1] + tmp_array[2] + gru.bz)

    mul!(tmp_array[4], gru.Wr_in, y)
    mul!(tmp_array[5], gru.Wr, x)
    @. tmp_array[6] = gru.activation_function[2](tmp_array[4] + tmp_array[5] + gru.br)

    mul!(tmp_array[7], W_in, y)
    mul!(tmp_array[8], W, tmp_array[6].*x)
    @. tmp_array[9] = gru.activation_function[3](tmp_array[7] + tmp_array[8] + b)
    @. out = (1 - tmp_array[3]) * x + tmp_array[3] * tmp_array[9]
end

#minimal
function obtain_gru_state!(out, variant::Minimal, gru, x, y, W, W_in, b, tmp_array)

    mul!(tmp_array[1], gru.Wz_in, y)
    mul!(tmp_array[2], gru.Wz, x)
    @. tmp_array[3] = gru.activation_function[1](tmp_array[1] + tmp_array[2] + gru.bz)

    mul!(tmp_array[4], W_in, y)
    mul!(tmp_array[5], W, tmp_array[3].*x)
    @. tmp_array[6] =  gru.activation_function[2](tmp_array[4] + tmp_array[5] + b)

    @. out = (1 - tmp_array[3]) * x + tmp_array[3] * tmp_array[6]
end

max_stars_repo_name: MartinuzziFrancesco/EchoStateNetwork

max_issues_repo_name: MartinuzziFrancesco/EchoStateNetwork

max_forks_repo_name: MartinuzziFrancesco/EchoStateNetwork


text:
include("../utils.jl")

@testset "GLEScript" begin
    g = G.GS()
    G.gsave(g)
    G.grestore(g)
    @test isapproxstr(String(take!(g)), "gsave grestore")

    "foo" |> g
    @test String(take!(g)) // "foo"

    g2 = G.GS()
    "foo" |> (g, g2)
    @test String(take!(g)) // String(take!(g2))

    "bar" |> g2
    g2 |> g
    @test String(take!(g)) // "bar"
end


text:
module SubMod1
function test_me(x)
    return x^2
end
end

max_stars_repo_name: JeffreySarnoff/techytok-examples

max_issues_repo_name: JeffreySarnoff/techytok-examples

max_forks_repo_name: JeffreySarnoff/techytok-examples


text:
# Before running this, please make sure to activate and instantiate the environment
# corresponding to [this `Project.toml`](https://raw.githubusercontent.com/alan-turing-institute/DataScienceTutorials.jl/master/Project.toml) and [this `Manifest.toml`](https://raw.githubusercontent.com/alan-turing-institute/DataScienceTutorials.jl/master/Manifest.toml)
# so that you get an environment which matches the one used to generate the tutorials:
#
# ```julia
# cd("DataScienceTutorials") # cd to folder with the *.toml
# using Pkg; Pkg.activate("."); Pkg.instantiate()
# ```

# ## Generating dummy data
# Let's start by generating some dummy data with both numerical values and categorical values:
using MLJ
using PrettyPrinting

@load KNNRegressor
# input
X = (age    = [23, 45, 34, 25, 67],
     gender = categorical(['m', 'm', 'f', 'm', 'f']))
# target
height = [178, 194, 165, 173, 168];

# Note that the scientific type of `age` is `Count` here:
scitype(X.age)

# We will want to coerce that to `Continuous` so that it can be given to a regressor that expects such values.
# ## Declaring a pipeline
# A typical workflow for such data is to one-hot-encode the categorical data and then apply some regression model on the data.# Let's say that we want to apply the following steps:# 1. standardize the target variable (`:height`)# 1. one hot encode the categorical data# 1. train a KNN regression model
# The `@pipeline` macro helps you define such a simple (non-branching) pipeline of steps to be applied in order:
pipe = @pipeline(X -> coerce(X, :age=>Continuous),
                OneHotEncoder(),
                KNNRegressor(K=3),
                target = UnivariateStandardizer());

# Note the coercion of the `:age` variable to Continuous since `KNNRegressor` expects `Continuous` input.# Note also the `target` keyword where you can specify a transformation of the target variable.
# Hyperparameters of this pipeline can be accessed (and set) using dot syntax:
pipe.knn_regressor.K = 2
pipe.one_hot_encoder.drop_last = true;

# Evaluation for a pipe can be done with the `evaluate!` method; implicitly it will construct machines that will contain the fitted parameters etc:
evaluate(pipe, X, height, resampling=Holdout(),
         measure=rms) |> pprint

# This file was generated using Literate.jl, https://github.com/fredrikekre/Literate.jl


text:
module Virtuals

using Base.Broadcast: Broadcasted, Extruded
using LinearAlgebra

using Swizzles.ArrayifiedArrays
using Swizzles.Properties

export Virtual, VirtualArray, virtualize, virtual

struct Virtual{T}
    ex
end

function virtual(root, Data::Type)
    ins = instance(Data)
    if ins !== nothing
        return something(ins)
    else
        return Virtual{Data}(root)
    end
end

Base.ndims(::Virtual{T}) where {T} = ndims(T)
Base.ndims(::Virtual{<:Tuple}) = 1
Base.ndims(::Type{Virtual{T}}) where {T} = ndims(T)
Base.ndims(::Type{<:Virtual{<:Tuple}}) = 1
Base.axes(::Virtual{T}) where {T} = ndims(T) == 0 ? () : error() #FIXME return virtual axes!
Base.Broadcast.broadcastable(x::Virtual{<:Union{AbstractArray,Number,Ref,Tuple,Broadcasted}}) = x
Base.eltype(::Virtual{T}) where {T} = eltype(T)

struct VirtualArray{T, N, Data<:AbstractArray{T, N}} <: AbstractArray{T, N}
    ex
end

function virtualtuple(root, Data::Type{<:Tuple})
    data = map(((i, Datum),) -> virtualize(:($root[$i]), Datum), enumerate(Data.parameters))
    return (data...,)
end

function virtualize(root, Data::Type)
    ins = instance(Data)
    if ins !== nothing
        return something(ins)
    else
        return Virtual{Data}(root)
    end
end

function virtualize(root, Data::AbstractArray{T, N}) where {T, N}
    return VirtualArray{T, N, Data}(root)
end

function virtualize(root, ::Type{<:Adjoint{T, Arg}}) where {T, Arg}
    arg = virtualize(:($root.parent), Arg)
    return Adjoint{T, typeof(arg)}(arg)
end

function virtualize(root, ::Type{<:Transpose{T, Arg}}) where {T, Arg}
    arg = virtualize(:($root.parent), Arg)
    return Transpose{T, typeof(arg)}(arg)
end

function virtualize(root, ::Type{<:ArrayifiedArray{T, N, Arg}}) where {T, N, Arg}
    arg = virtualize(:($root.arg), Arg)
    return ArrayifiedArray{T, N, typeof(arg)}(arg)
end

function virtualize(root, Data::Type{<:AbstractArray{T, N}}) where {T, N}
    try
        return Data.instance
    catch Error
        return VirtualArray{T, N, Data}(root)
    end
end

function virtualize(root, ::Type{<:Broadcasted{Style, Axes, F, Args}}) where {Style, Axes, F, Args}
    f = virtual(:($root.f), F)
    args = (map(((i, Arg),) -> virtualize(:($root.args[$i]), Arg), enumerate(Args.parameters))...,)
    axes = (map(((i, Axis),) -> virtualize(:($root.axes[$i]), Axes), enumerate(Axes.parameters))...,)
    return Broadcasted{Style}(f, args, axes)
end

function virtualize(root, ::Type{<:Extruded{Arg, Keeps, Defaults}}) where {Arg, Keeps, Defaults}
    arg = virtualize(:($root.arg), Arg)
    args = (map(((i, K),) -> virtual(:($root.keeps[$i]), K), enumerate(Keeps.parameters))...,)
    axes = (map(((i, D),) -> virtual(:($root.defaults[$i]), D), enumerate(Defaults.parameters))...,)
    return Extruded(arg, keeps, defaults)
end

end


text:
using OptimalTransport

using CUDA
using Distances
using NNlibCUDA

using LinearAlgebra
using Random
using Test

CUDA.allowscalar(false)

Random.seed!(100)

@testset "simple_gpu.jl" begin
    # ensure that a GPU is available
    if !CUDA.functional()
        @warn "skipped GPU tests: no GPU available"
    else
        # source histogram
        m = 200
        μ = normalize!(rand(Float32, m), 1)
        cu_μ = cu(μ)

        # target histogram
        n = 250
        ν = normalize!(rand(Float32, n), 1)
        cu_ν = cu(ν)

        # random cost matrix
        C = pairwise(SqEuclidean(), randn(Float32, 1, m), randn(Float32, 1, n); dims=2)
        cu_C = cu(C)

        # regularization parameter
        ε = 0.01f0

        @testset "sinkhorn" begin
            for alg in (
                (),
                (SinkhornGibbs(),),
                (SinkhornStabilized(),),
                (SinkhornEpsilonScaling(SinkhornGibbs()),),
                (SinkhornEpsilonScaling(SinkhornStabilized()),),
            )
                # compute transport plan and cost on the GPU
                γ = sinkhorn(cu_μ, cu_ν, cu_C, ε, alg...)
                @test γ isa CuArray{Float32,2}
                c = sinkhorn2(cu_μ, cu_ν, cu_C, ε, alg...)
                @test c isa Float32

                # compare with results on the CPU
                γ_cpu = sinkhorn(μ, ν, C, ε, alg...)
                @test convert(Array, γ) ≈ γ_cpu
                @test c ≈ sinkhorn2(μ, ν, C, ε, alg...)

                # batches of histograms
                d = 10
                for (size2_μ, size2_ν) in
                    (((), (d,)), ((1,), (d,)), ((d,), ()), ((d,), (1,)), ((d,), (d,)))
                    # generate uniform histograms
                    μ_batch = repeat(cu_μ, 1, size2_μ...)
                    ν_batch = repeat(cu_ν, 1, size2_ν...)

                    # compute optimal transport plan and check that it is consistent with the
                    # plan for individual histograms
                    γ_all = sinkhorn(μ_batch, ν_batch, cu_C, ε, alg...)
                    @test γ_all isa CuArray{Float32,3}
                    @test size(γ_all) == (m, n, d)
                    @test all(γi ≈ γ_cpu for γi in eachslice(convert(Array, γ_all); dims=3))

                    # compute optimal transport cost and check that it is consistent with the
                    # cost for individual histograms
                    c_all = sinkhorn2(μ_batch, ν_batch, cu_C, ε, alg...)
                    @test c_all isa CuArray{Float32,1}
                    @test size(c_all) == (d,)
                    @test all(ci ≈ c for ci in convert(Array, c_all))
                end
            end
        end

        @testset "sinkhorn_unbalanced" begin
            # source histogram
            μscaled = 1.5f0 .* μ
            cu_μscaled = cu(μscaled)

            # compute transport plan and cost on the GPU
            λ1 = 0.4f0
            λ2 = 0.6f0
            γ = sinkhorn_unbalanced(cu_μscaled, cu_ν, cu_C, λ1, λ2, ε)
            c = sinkhorn_unbalanced2(cu_μscaled, cu_ν, cu_C, λ1, λ2, ε)

            # compare with results on the CPU
            @test convert(Array, γ) ≈ sinkhorn_unbalanced(μscaled, ν, C, λ1, λ2, ε)
            @test c ≈ sinkhorn_unbalanced2(μscaled, ν, C, λ1, λ2, ε)
        end

        @testset "quadreg" begin
            # use a different reg parameter
            ε_quad = 1.0f0
            γ = quadreg(cu_μ, cu_ν, cu_C, ε_quad, QuadraticOTNewton(0.1f0, 0.5f0, 1f-5, 50))
            # compare with results on the CPU
            @test convert(Array, γ) ≈
                  quadreg(μ, ν, C, ε_quad, QuadraticOTNewton(0.1f0, 0.5f0, 1f-5, 50)) atol =
                1f-4 rtol = 1f-4
        end
    end
end


text:
# Julia wrapper for header: /usr/include/scip/scip_param.h
# Automatically generated using Clang.jl wrap_c

function SCIPaddBoolParam(scip, name, desc, valueptr, isadvanced, defaultvalue, paramchgd, paramdata)
    ccall((:SCIPaddBoolParam, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Cstring, Cstring, Ptr{UInt32}, UInt32, UInt32, Ptr{Cvoid}, Ptr{SCIP_PARAMDATA}), scip, name, desc, valueptr, isadvanced, defaultvalue, paramchgd, paramdata)
end

function SCIPaddIntParam(scip, name, desc, valueptr, isadvanced, defaultvalue, minvalue, maxvalue, paramchgd, paramdata)
    ccall((:SCIPaddIntParam, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Cstring, Cstring, Ptr{Cint}, UInt32, Cint, Cint, Cint, Ptr{Cvoid}, Ptr{SCIP_PARAMDATA}), scip, name, desc, valueptr, isadvanced, defaultvalue, minvalue, maxvalue, paramchgd, paramdata)
end

function SCIPaddLongintParam(scip, name, desc, valueptr, isadvanced, defaultvalue, minvalue, maxvalue, paramchgd, paramdata)
    ccall((:SCIPaddLongintParam, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Cstring, Cstring, Ptr{Clonglong}, UInt32, Clonglong, Clonglong, Clonglong, Ptr{Cvoid}, Ptr{SCIP_PARAMDATA}), scip, name, desc, valueptr, isadvanced, defaultvalue, minvalue, maxvalue, paramchgd, paramdata)
end

function SCIPaddRealParam(scip, name, desc, valueptr, isadvanced, defaultvalue, minvalue, maxvalue, paramchgd, paramdata)
    ccall((:SCIPaddRealParam, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Cstring, Cstring, Ptr{Cdouble}, UInt32, Cdouble, Cdouble, Cdouble, Ptr{Cvoid}, Ptr{SCIP_PARAMDATA}), scip, name, desc, valueptr, isadvanced, defaultvalue, minvalue, maxvalue, paramchgd, paramdata)
end

function SCIPaddCharParam(scip, name, desc, valueptr, isadvanced, defaultvalue, allowedvalues, paramchgd, paramdata)
    ccall((:SCIPaddCharParam, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Cstring, Cstring, Cstring, UInt32, UInt8, Cstring, Ptr{Cvoid}, Ptr{SCIP_PARAMDATA}), scip, name, desc, valueptr, isadvanced, defaultvalue, allowedvalues, paramchgd, paramdata)
end

function SCIPaddStringParam(scip, name, desc, valueptr, isadvanced, defaultvalue, paramchgd, paramdata)
    ccall((:SCIPaddStringParam, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Cstring, Cstring, Ptr{Cstring}, UInt32, Cstring, Ptr{Cvoid}, Ptr{SCIP_PARAMDATA}), scip, name, desc, valueptr, isadvanced, defaultvalue, paramchgd, paramdata)
end

function SCIPisParamFixed(scip, name)
    ccall((:SCIPisParamFixed, libscip), UInt32, (Ptr{SCIP_}, Cstring), scip, name)
end

function SCIPgetParam(scip, name)
    ccall((:SCIPgetParam, libscip), Ptr{SCIP_PARAM}, (Ptr{SCIP_}, Cstring), scip, name)
end

function SCIPgetBoolParam(scip, name, value)
    ccall((:SCIPgetBoolParam, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Cstring, Ptr{UInt32}), scip, name, value)
end

function SCIPgetIntParam(scip, name, value)
    ccall((:SCIPgetIntParam, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Cstring, Ptr{Cint}), scip, name, value)
end

function SCIPgetLongintParam(scip, name, value)
    ccall((:SCIPgetLongintParam, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Cstring, Ptr{Clonglong}), scip, name, value)
end

function SCIPgetRealParam(scip, name, value)
    ccall((:SCIPgetRealParam, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Cstring, Ptr{Cdouble}), scip, name, value)
end

function SCIPgetCharParam(scip, name, value)
    ccall((:SCIPgetCharParam, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Cstring, Ptr{Cchar}), scip, name, value)
end

function SCIPgetStringParam(scip, name, value)
    ccall((:SCIPgetStringParam, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Cstring, Ptr{Cstring}), scip, name, value)
end

function SCIPfixParam(scip, name)
    ccall((:SCIPfixParam, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Cstring), scip, name)
end

function SCIPunfixParam(scip, name)
    ccall((:SCIPunfixParam, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Cstring), scip, name)
end

function SCIPsetParam(scip, name, value)
    ccall((:SCIPsetParam, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Cstring, Ptr{Cvoid}), scip, name, value)
end

function SCIPchgBoolParam(scip, param, value)
    ccall((:SCIPchgBoolParam, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_PARAM}, UInt32), scip, param, value)
end

function SCIPsetBoolParam(scip, name, value)
    ccall((:SCIPsetBoolParam, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Cstring, UInt32), scip, name, value)
end

function SCIPisBoolParamValid(scip, param, value)
    ccall((:SCIPisBoolParamValid, libscip), UInt32, (Ptr{SCIP_}, Ptr{SCIP_PARAM}, UInt32), scip, param, value)
end

function SCIPchgIntParam(scip, param, value)
    ccall((:SCIPchgIntParam, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_PARAM}, Cint), scip, param, value)
end

function SCIPsetIntParam(scip, name, value)
    ccall((:SCIPsetIntParam, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Cstring, Cint), scip, name, value)
end

function SCIPisIntParamValid(scip, param, value)
    ccall((:SCIPisIntParamValid, libscip), UInt32, (Ptr{SCIP_}, Ptr{SCIP_PARAM}, Cint), scip, param, value)
end

function SCIPchgLongintParam(scip, param, value)
    ccall((:SCIPchgLongintParam, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_PARAM}, Clonglong), scip, param, value)
end

function SCIPsetLongintParam(scip, name, value)
    ccall((:SCIPsetLongintParam, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Cstring, Clonglong), scip, name, value)
end

function SCIPisLongintParamValid(scip, param, value)
    ccall((:SCIPisLongintParamValid, libscip), UInt32, (Ptr{SCIP_}, Ptr{SCIP_PARAM}, Clonglong), scip, param, value)
end

function SCIPchgRealParam(scip, param, value)
    ccall((:SCIPchgRealParam, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_PARAM}, Cdouble), scip, param, value)
end

function SCIPsetRealParam(scip, name, value)
    ccall((:SCIPsetRealParam, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Cstring, Cdouble), scip, name, value)
end

function SCIPisRealParamValid(scip, param, value)
    ccall((:SCIPisRealParamValid, libscip), UInt32, (Ptr{SCIP_}, Ptr{SCIP_PARAM}, Cdouble), scip, param, value)
end

function SCIPchgCharParam(scip, param, value)
    ccall((:SCIPchgCharParam, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_PARAM}, UInt8), scip, param, value)
end

function SCIPsetCharParam(scip, name, value)
    ccall((:SCIPsetCharParam, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{Cchar}, UInt8), scip, name, value)
end

function SCIPisCharParamValid(scip, param, value)
    ccall((:SCIPisCharParamValid, libscip), UInt32, (Ptr{SCIP_}, Ptr{SCIP_PARAM}, UInt8), scip, param, value)
end

function SCIPchgStringParam(scip, param, value)
    ccall((:SCIPchgStringParam, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_PARAM}, Cstring), scip, param, value)
end

function SCIPsetStringParam(scip, name, value)
    ccall((:SCIPsetStringParam, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Cstring, Cstring), scip, name, value)
end

function SCIPisStringParamValid(scip, param, value)
    ccall((:SCIPisStringParamValid, libscip), UInt32, (Ptr{SCIP_}, Ptr{SCIP_PARAM}, Cstring), scip, param, value)
end

function SCIPreadParams(scip, filename)
    ccall((:SCIPreadParams, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Cstring), scip, filename)
end

function SCIPwriteParam(scip, param, filename, comments, onlychanged)
    ccall((:SCIPwriteParam, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_PARAM}, Cstring, UInt32, UInt32), scip, param, filename, comments, onlychanged)
end

function SCIPwriteParams(scip, filename, comments, onlychanged)
    ccall((:SCIPwriteParams, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Cstring, UInt32, UInt32), scip, filename, comments, onlychanged)
end

function SCIPresetParam(scip, name)
    ccall((:SCIPresetParam, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Cstring), scip, name)
end

function SCIPresetParams(scip)
    ccall((:SCIPresetParams, libscip), SCIP_RETCODE, (Ptr{SCIP_},), scip)
end

function SCIPsetEmphasis(scip, paramemphasis, quiet)
    ccall((:SCIPsetEmphasis, libscip), SCIP_RETCODE, (Ptr{SCIP_}, SCIP_PARAMEMPHASIS, UInt32), scip, paramemphasis, quiet)
end

function SCIPsetSubscipsOff(scip, quiet)
    ccall((:SCIPsetSubscipsOff, libscip), SCIP_RETCODE, (Ptr{SCIP_}, UInt32), scip, quiet)
end

function SCIPsetHeuristics(scip, paramsetting, quiet)
    ccall((:SCIPsetHeuristics, libscip), SCIP_RETCODE, (Ptr{SCIP_}, SCIP_PARAMSETTING, UInt32), scip, paramsetting, quiet)
end

function SCIPsetPresolving(scip, paramsetting, quiet)
    ccall((:SCIPsetPresolving, libscip), SCIP_RETCODE, (Ptr{SCIP_}, SCIP_PARAMSETTING, UInt32), scip, paramsetting, quiet)
end

function SCIPsetSeparating(scip, paramsetting, quiet)
    ccall((:SCIPsetSeparating, libscip), SCIP_RETCODE, (Ptr{SCIP_}, SCIP_PARAMSETTING, UInt32), scip, paramsetting, quiet)
end

function SCIPgetParams(scip)
    ccall((:SCIPgetParams, libscip), Ptr{Ptr{SCIP_PARAM}}, (Ptr{SCIP_},), scip)
end

function SCIPgetNParams(scip)
    ccall((:SCIPgetNParams, libscip), Cint, (Ptr{SCIP_},), scip)
end


text:
# Manage and access multinetwork dimensions

## Dimension data structure

function _initialize_dim()
    Dict{Symbol,Any}(
        :pos    => NamedTuple(), # Position (order) of each dimension
        :prop   => Dict{Symbol,Dict{Int,Dict{String,Any}}}(), # Data relating to the elements of the dimensions
        :meta   => Dict{Symbol,Any}(), # Data relating to dimensions
        :offset => 0 # Offset of nw ids: the id of the first nw is offset+1
        # li (linear indices): generated by add_dimension!()
        # ci (cartesian indices): generated by add_dimension!()
    )
end

"""
    add_dimension!(data, name, properties; metadata)

Add dimension `name` to `data` specifying some `properties` of the dimension ids.

# Arguments
- `data::Dict{String,Any}`: a single-network data dictionary.
- `name::Symbol`: the name to use to refer to the dimension in the code.
- `properties::Dict{Int,Dict{String,Any}}`: properties associated to the dimension ids. The
  outer dictionary is indexed with the ids along the dimension (consecutive `Int`s starting
  from 1). The inner dictionaries, one for each id, store the properties.
- `metadata::Dict{String,Any} = Dict{String,Any}()`: optional metadata describing the
  dimension as a whole.

# Examples
```julia-repl
julia> add_dimension!(data, :scenario, Dict(s => Dict{String,Any}("probability"=>1/4) for s in 1:4))
```

# Extended help
The functions `dim_prop`, `dim_meta` and `dim_length` can be used to access properties,
metadata and length (cardinality) of a dimension. They apply to both data dictionaries
(`Dict{String,Any}`) and powermodels (`PowerModels.AbstractPowerModel`).
"""
function add_dimension!(data::Dict{String,Any}, name::Symbol, properties::Dict{Int,Dict{String,Any}}; metadata::Dict{String,Any}=Dict{String,Any}())
    dim = get!(data, "dim", _initialize_dim())

    if haskey(dim[:pos], name)
        Memento.error(_LOGGER, "A dimension named \"$name\" is already present in data.")
    end
    if Set(keys(properties)) != Set(1:length(properties))
        Memento.error(_LOGGER, "Keys of `properties` Dict must range from 1 to the number of $(name)s.")
    end

    dim[:pos] = (; dim[:pos]..., name => length(dim[:pos])+1)
    dim[:prop][name] = properties
    dim[:meta][name] = metadata
    dim[:li] = LinearIndices(Tuple(1:length(dim[:prop][nm]) for nm in keys(dim[:pos]))).+dim[:offset]
    dim[:ci] = CartesianIndices(dim[:li])

    return dim
end

"""
    add_dimension!(data, name, size; metadata)

Add dimension `name` to `data` specifying the dimension `size`.

# Examples
```julia-repl
julia> add_dimension!(data, :hour, 24)
```
"""
function add_dimension!(data::Dict{String,Any}, name::Symbol, size::Int; metadata::Dict{String,Any}=Dict{String,Any}())
    properties = Dict{Int,Dict{String,Any}}(i => Dict{String,Any}() for i in 1:size)
    add_dimension!(data, name, properties; metadata)
end

"""
    shift_ids!(dim::Dict{Symbol,Any}, offset)
    shift_ids!(data::Dict{String,Any}, offset)

Shift by `offset` the network ids in `dim` or `data`.

The `offset` argument is added to the existing offset.
Return a vector containing the shifted network ids.
`data` must be a single-network `data` dictionary.
"""
function shift_ids! end

function shift_ids!(dim::Dict{Symbol,Any}, offset::Int)
    dim[:offset] += offset
    dim[:li] .+= offset
    return vec(dim[:li])
end

function shift_ids!(data::Dict{String,Any}, offset::Int)
    if _IM.ismultinetwork(data)
        Memento.error(_LOGGER, "`shift_ids!` can only be applied to single-network data dictionaries.")
    end
    shift_ids!(dim(data), offset)
end

"""
    merge_dim!(dim1, dim2, dimension)

Merge `dim1` and `dim2` structures along `dimension`.

The ids of `dim2` must be contiguous to those of `dim1`.
"""
function merge_dim!(dim1::Dict{Symbol,Any}, dim2::Dict{Symbol,Any}, dimension::Symbol)
    dim = Dict{Symbol,Any}()

    if dim1[:pos] != dim2[:pos]
        Memento.error(_LOGGER, "The dimensions to be merged have different names and/or order:\nfirst:  $(dim1[:pos])\nsecond: $(dim1[:pos])")
    end
    dim[:pos] = dim1[:pos]

    if any(dim1[:prop][d] != dim2[:prop][d] for d in delete!(Set(keys(dim[:pos])), dimension))
        diff = join(", $d" for d in delete!(Set(keys(dim[:pos])), dimension) if dim1[:prop][d] != dim2[:prop][d])[2:end]
        Memento.error(_LOGGER, "Different properties found in the following dimension(s):$diff.")
    end
    dim[:prop] = dim1[:prop]
    offset = length(dim1[:prop][dimension])
    for (k,v) in dim2[:prop][dimension]
        dim[:prop][dimension][k+offset] = deepcopy(v)
    end

    if dim1[:meta] != dim2[:meta]
        diff = join(", $d" for d in keys(dim[:pos]) if dim1[:meta][d] != dim2[:meta][d])[2:end]
        Memento.error(_LOGGER, "Different metadata found in the following dimension(s):$diff.")
    end
    dim[:meta] = dim1[:meta]

    if dim2[:li][1] != dim1[:li][end]+1
        Memento.error(_LOGGER, "Multinetworks to be merged must have contiguous ids.")
    end
    dim[:offset] = min(dim1[:offset], dim2[:offset])
    dim[:li] = LinearIndices(Tuple(1:length(dim[:prop][nm]) for nm in keys(dim[:pos]))).+dim[:offset]
    dim[:ci] = CartesianIndices(dim[:li])

    return dim
end

"""
    slice, ids = slice_dim(dim::Dict{Symbol,Any}; kwargs...)

Slice `dim` structure keeping the networks that have the coordinates specified by `kwargs`.

`kwargs` must be of the form `name = <value>`, where `name` is the name of a dimension of
`dim` and `<value>` is an `Int` coordinate of that dimension.

Return `slice`, a sliced `dim` structure whose networks have ids starting at 1, and `ids`, a
vector containing the ids that the networks making up `slice` have in `dim`.
The coordinates of the dimensions at which `dim` is sliced are accessible with
`dim_meta(slice, <name>, "orig_id")`, where `<name>` is the name of one of those dimensions.

# Examples
```julia-repl
julia> slice_dim(dim; hour = 24)
julia> slice_dim(dim; hour = 24, scenario = 3)
```
"""
function slice_dim(dim::Dict{Symbol,Any}; kwargs...)
    slice = Dict{Symbol,Any}()
    slice[:pos] = dim[:pos]
    slice[:prop] = Dict{Symbol,Dict{Int,Dict{String,Any}}}()
    for d in keys(dim[:pos])
        if d ∈ keys(kwargs)
            slice[:prop][d] = Dict{Int,Dict{String,Any}}(1 => deepcopy(dim[:prop][d][kwargs[d]]))
        else
            slice[:prop][d] = deepcopy(dim[:prop][d])
        end
    end
    slice[:meta] = deepcopy(dim[:meta])
    for (d, i) in kwargs
        slice[:meta][d]["orig_id"] = i
    end
    slice[:offset] = 0
    slice[:li] = collect(LinearIndices(Tuple(1:length(slice[:prop][nm]) for nm in keys(slice[:pos]))))
    slice[:ci] = CartesianIndices(slice[:li])

    names = keys(dim[:pos])
    li = dim[:li]
    ids = vec(li[ntuple(i -> get(kwargs, names[i], axes(li,i)), ndims(li))...])

    return slice, ids
end

"""
    require_dim(data, dimensions...)

Verify that the specified `dimensions` are present in `data`; if not, raise an error.
"""
function require_dim(data::Dict{String,Any}, dimensions::Symbol...)
    if !haskey(data, "dim")
        Memento.error(_LOGGER, "Missing `dim` dict in `data`. Use `add_dimension!` to fix.")
    end
    for d in dimensions
        if !haskey(dim(data)[:prop], d)
            Memento.error(_LOGGER, "Missing dimension \"$d\" in `data`. Use `add_dimension!` to fix.")
        end
    end
end

## Access (subsets of) nw ids

"""
    nw_ids(pm::PowerModels.AbstractPowerModel; kwargs...)
    nw_ids(data::Dict{String,Any}; kwargs...)
    nw_ids(dim::Dict{Symbol,Any}; kwargs...)

Sorted vector containing nw ids of `pm`, `data`, or `dim`, optionally filtered by the coordinates of one or more dimensions.

`kwargs` must be of the form `name = <value>` or `name = <interval>` or `name = <subset>`,
where `name` is the name of a dimension of `pm`.

# Examples
```julia-repl
julia> nw_ids(pm)
julia> nw_ids(pm; hour = 24)
julia> nw_ids(pm; hour = 13:24)
julia> nw_ids(pm; hour = [6,12,18,24])
julia> nw_ids(pm; hour = 24, scenario = 3)
```
"""
function nw_ids end

function nw_ids(dim::Dict{Symbol,Any}; kwargs...)::Vector{Int}
    names = keys(dim[:pos])
    li = dim[:li]
    nws = li[ntuple(i -> get(kwargs, names[i], axes(li, i)), ndims(li))...]
    ndims(nws) >= 1 ? vec(nws) : [nws]
end

function nw_ids(data::Dict{String,Any}; kwargs...)::Vector{Int}
    haskey(data, "dim") ? nw_ids(dim(data); kwargs...) : [0]
end

function nw_ids(pm::_PM.AbstractPowerModel; kwargs...)::Vector{Int}
    haskey(pm.ref[:it][_PM.pm_it_sym], :dim) ? nw_ids(dim(pm); kwargs...) : [0]
end

## Compute nw ids given another nw id

"""
    similar_ids(pm::PowerModels.AbstractPowerModel, n::Int; kwargs...)
    similar_ids(data::Dict{String,Any}, n::Int; kwargs...)
    similar_ids(dim::Dict{Symbol,Any}, n::Int; kwargs...)

Sorted vector containing nw ids that have the same coordinates along dimensions as `n` except for dimensions passed in `kwargs`.

`kwargs` must be of the form `name = <value>` or `name = <interval>` or `name = <subset>`,
where `name` is the name of a dimension of `pm`, `data` or `dim`.

# Examples
```julia-repl
julia> similar_ids(pm, 3; hour = 24)
julia> similar_ids(pm, 3; hour = 13:24)
julia> similar_ids(pm, 3; hour = [6,12,18,24])
julia> similar_ids(pm, 3; hour = 24, scenario = 3)
```
"""
function similar_ids end

function similar_ids(dim::Dict{Symbol,Any}, n::Int; kwargs...)::Vector{Int}
    names = keys(dim[:pos])
    offset = dim[:offset]
    li = dim[:li]
    ci_n = dim[:ci][n-offset]
    nws = li[ntuple(i -> get(kwargs, names[i], ci_n[i]), ndims(li))...]
    ndims(nws) >= 1 ? vec(nws) : [nws]
end

similar_ids(data::Dict{String,Any}, n::Int; kwargs...) = similar_ids(dim(data), n; kwargs...)
similar_ids(pm::_PM.AbstractPowerModel, n::Int; kwargs...) = similar_ids(dim(pm), n; kwargs...)

"""
    similar_id(pm::PowerModels.AbstractPowerModel, n::Int; kwargs...)
    similar_id(data::Dict{String,Any}, n::Int; kwargs...)
    similar_id(dim::Dict{Symbol,Any}, n::Int; kwargs...)

Nw id that has the same coordinates along dimensions as `n` except for dimensions passed in `kwargs`.

`kwargs` must be of the form `name = <value>`, where `name` is the name of a dimension of
`pm`, `data` or `dim`.

# Examples
```julia-repl
julia> similar_id(pm, 3; hour = 24)
julia> similar_id(pm, 3; hour = 24, scenario = 3)
```
"""
function similar_id end

function similar_id(dim::Dict{Symbol,Any}, n::Int; kwargs...)::Int
    names = keys(dim[:pos])
    offset = dim[:offset]
    li = dim[:li]
    ci_n = dim[:ci][n-offset]
    li[ntuple(i -> get(kwargs, names[i], ci_n[i])::Int, ndims(li))...]
end

similar_id(data::Dict{String,Any}, n::Int; kwargs...) = similar_id(dim(data), n; kwargs...)
similar_id(pm::_PM.AbstractPowerModel, n::Int; kwargs...) = similar_id(dim(pm), n; kwargs...)

"""
    first_id(pm::PowerModels.AbstractPowerModel, n::Int, dimension::Symbol...)
    first_id(data::Dict{String,Any}, n::Int, dimension::Symbol...)
    first_id(dim::Dict{Symbol,Any}, n::Int, dimension::Symbol...)

Return the first network in `pm`, `data` or `dim` along `dimension` while keeping the other dimensions fixed.
"""
function first_id end

function first_id(dim::Dict{Symbol,Any}, n::Int, dimension::Symbol...)
    names = keys(dim[:pos])
    offset = dim[:offset]
    li = dim[:li]
    ci_n = dim[:ci][n-offset]
    li[ntuple(i -> names[i] in dimension ? 1 : ci_n[i], ndims(li))...]
end

first_id(data::Dict{String,Any}, n::Int, dimension::Symbol...) = first_id(dim(data), n, dimension...)
first_id(pm::_PM.AbstractPowerModel, n::Int, dimension::Symbol...) = first_id(dim(pm), n, dimension...)

"""
    last_id(pm::PowerModels.AbstractPowerModel, n::Int, dimension::Symbol...)
    last_id(data::Dict{String,Any}, n::Int, dimension::Symbol...)
    last_id(dim::Dict{Symbol,Any}, n::Int, dimension::Symbol...)

Return the last network in `pm`, `data` or `dim` along `dimension` while keeping the other dimensions fixed.
"""
function last_id end

function last_id(dim::Dict{Symbol,Any}, n::Int, dimension::Symbol...)
    names = keys(dim[:pos])
    offset = dim[:offset]
    li = dim[:li]
    ci_n = dim[:ci][n-offset]
    li[ntuple(i -> names[i] in dimension ? size(li,i) : ci_n[i], ndims(li))...]
end

last_id(data::Dict{String,Any}, n::Int, dimension::Symbol...) = last_id(dim(data), n, dimension...)
last_id(pm::_PM.AbstractPowerModel, n::Int, dimension::Symbol...) = last_id(dim(pm), n, dimension...)

"""
    prev_id(pm::PowerModels.AbstractPowerModel, n::Int, dimension::Symbol)
    prev_id(data::Dict{String,Any}, n::Int, dimension::Symbol)
    prev_id(dim::Dict{Symbol,Any}, n::Int, dimension::Symbol)

Return the previous network in `pm`, `data` or `dim` along `dimension` while keeping the other dimensions fixed.
"""
function prev_id end

function prev_id(dim::Dict{Symbol,Any}, n::Int, dimension::Symbol)::Int
    pos_d = dim[:pos][dimension]
    offset = dim[:offset]
    li = dim[:li]
    ci_n = dim[:ci][n-offset]
    li[ntuple(i -> i == pos_d ? ci_n[i]-1 : ci_n[i], ndims(li))...]
end

prev_id(data::Dict{String,Any}, n::Int, dimension::Symbol) = prev_id(dim(data), n, dimension)
prev_id(pm::_PM.AbstractPowerModel, n::Int, dimension::Symbol) = prev_id(dim(pm), n, dimension)

"""
    prev_ids(pm::PowerModels.AbstractPowerModel, n::Int, dimension::Symbol)
    prev_ids(data::Dict{String,Any}, n::Int, dimension::Symbol)
    prev_ids(dim::Dict{Symbol,Any}, n::Int, dimension::Symbol)

Return the previous networks in `pm`, `data` or `dim` along `dimension` while keeping the other dimensions fixed.
"""
function prev_ids end

function prev_ids(dim::Dict{Symbol,Any}, n::Int, dimension::Symbol)::Vector{Int}
    pos_d = dim[:pos][dimension]
    offset = dim[:offset]
    li = dim[:li]
    ci_n = Tuple(dim[:ci][n-offset])
    li[CartesianIndex(ci_n[1:pos_d-1]), CartesianIndices((1:ci_n[pos_d]-1,)), CartesianIndex(ci_n[pos_d+1:end])]
end

prev_ids(data::Dict{String,Any}, n::Int, dimension::Symbol) = prev_ids(dim(data), n, dimension)
prev_ids(pm::_PM.AbstractPowerModel, n::Int, dimension::Symbol) = prev_ids(dim(pm), n, dimension)

"""
    next_id(pm::PowerModels.AbstractPowerModel, n::Int, dimension::Symbol)
    next_id(data::Dict{String,Any}, n::Int, dimension::Symbol)
    next_id(dim::Dict{Symbol,Any}, n::Int, dimension::Symbol)

Return the next network in `pm`, `data` or `dim` along `dimension` while keeping the other dimensions fixed.
"""
function next_id end

function next_id(dim::Dict{Symbol,Any}, n::Int, dimension::Symbol)::Int
    pos_d = dim[:pos][dimension]
    offset = dim[:offset]
    li = dim[:li]
    ci_n = dim[:ci][n-offset]
    li[ntuple(i -> i == pos_d ? ci_n[i]+1 : ci_n[i], ndims(li))...]
end

next_id(data::Dict{String,Any}, n::Int, dimension::Symbol) = next_id(dim(data), n, dimension)
next_id(pm::_PM.AbstractPowerModel, n::Int, dimension::Symbol) = next_id(dim(pm), n, dimension)

"""
    next_ids(pm::PowerModels.AbstractPowerModel, n::Int, dimension::Symbol)
    next_ids(data::Dict{String,Any}, n::Int, dimension::Symbol)
    next_ids(dim::Dict{Symbol,Any}, n::Int, dimension::Symbol)

Return the next networks in `pm`, `data` or `dim` along `dimension` while keeping the other dimensions fixed.
"""
function next_ids end

function next_ids(dim::Dict{Symbol,Any}, n::Int, dimension::Symbol)::Vector{Int}
    pos_d = dim[:pos][dimension]
    offset = dim[:offset]
    li = dim[:li]
    ci_n = Tuple(dim[:ci][n-offset])
    li[CartesianIndex(ci_n[1:pos_d-1]), CartesianIndices((ci_n[pos_d]+1:size(li,pos_d),)), CartesianIndex(ci_n[pos_d+1:end])]
end

next_ids(data::Dict{String,Any}, n::Int, dimension::Symbol) = next_ids(dim(data), n, dimension)
next_ids(pm::_PM.AbstractPowerModel, n::Int, dimension::Symbol) = next_ids(dim(pm), n, dimension)

## Query properties of nw ids

"""
    coord(pm::_PM.AbstractPowerModel, n::Int, dimension::Symbol)
    coord(data::Dict{String,Any}, n::Int, dimension::Symbol)
    coord(dim::Dict{Symbol,Any}, n::Int, dimension::Symbol)

Return the coordinate along `dimension` of nw `n` of `pm`, `data` or `dim`.
"""
function coord end

function coord(dim::Dict{Symbol,Any}, n::Int, dimension::Symbol)
    pos = dim[:pos]
    offset = dim[:offset]
    ci_n = dim[:ci][n-offset]
    ci_n[pos[dimension]]
end

coord(data::Dict{String,Any}, args...) = coord(dim(data), args...)
coord(pm::_PM.AbstractPowerModel, args...) = coord(dim(pm), args...)

"""
    is_first_id(pm::PowerModels.AbstractPowerModel, n::Int, dimension::Symbol...)
    is_first_id(data::Dict{String,Any}, n::Int, dimension::Symbol...)
    is_first_id(dim::Dict{Symbol,Any}, n::Int, dimension::Symbol...)

Return whether the network `n` is the first along `dimension` in `pm`, `data` or `dim`.
"""
function is_first_id end

function is_first_id(dim::Dict{Symbol,Any}, n::Int, dimension::Symbol...)
    pos = dim[:pos]
    offset = dim[:offset]
    ci_n = dim[:ci][n-offset]
    all(ci_n[pos[d]] == 1 for d in dimension)
end

is_first_id(data::Dict{String,Any}, args...) = is_first_id(dim(data), args...)
is_first_id(pm::_PM.AbstractPowerModel, args...) = is_first_id(dim(pm), args...)

"""
    is_last_id(pm::PowerModels.AbstractPowerModel, n::Int, dimension::Symbol...)
    is_last_id(data::Dict{String,Any}, n::Int, dimension::Symbol...)
    is_last_id(dim::Dict{Symbol,Any}, n::Int, dimension::Symbol...)

Return whether the network `n` is the last along `dimension` in `pm`, `data` or `dim`.
"""
function is_last_id(dim::Dict{Symbol,Any}, n::Int, dimension::Symbol...)
    pos = dim[:pos]
    offset = dim[:offset]
    li = dim[:li]
    ci_n = dim[:ci][n-offset]
    all(ci_n[pos[d]] == size(li,pos[d]) for d in dimension)
end

is_last_id(data::Dict{String,Any}, args...) = is_last_id(dim(data), args...)
is_last_id(pm::_PM.AbstractPowerModel, args...) = is_last_id(dim(pm), args...)

## Access data relating to dimensions

"""
    dim(data::Dict{String,Any})
    dim(pm::PowerModels.AbstractPowerModel)

Return `dim` data structure.
"""
function dim end
dim(data::Dict{String,Any}) = data["dim"]
dim(pm::_PM.AbstractPowerModel) = pm.ref[:it][_PM.pm_it_sym][:dim]

"""
    dim_prop(dim::Dict{Symbol,Any}[, dimension[, id[, key]]])
    dim_prop(data::Dict{String,Any}[, dimension[, id[, key]]])
    dim_prop(pm::PowerModels.AbstractPowerModel[, dimension[, id[, key]]])

Properties associated to the `id`s of `dimension`.

    dim_prop(dim::Dict{Symbol,Any}, n, dimension[, key])
    dim_prop(data::Dict{String,Any}, n, dimension[, key])
    dim_prop(pm::PowerModels.AbstractPowerModel, n, dimension[, key])

Properties associated to `dimension` of a network `n`.
"""
function dim_prop end
dim_prop(dim::Dict{Symbol,Any}) = dim[:prop]
dim_prop(dim::Dict{Symbol,Any}, dimension::Symbol) = dim[:prop][dimension]
dim_prop(dim::Dict{Symbol,Any}, dimension::Symbol, id::Int) = dim[:prop][dimension][id]
dim_prop(dim::Dict{Symbol,Any}, dimension::Symbol, id::Int, key::String) = dim[:prop][dimension][id][key]

dim_prop(dim::Dict{Symbol,Any}, n::Int, dimension::Symbol) = dim[:prop][dimension][coord(dim,n,dimension)]
dim_prop(dim::Dict{Symbol,Any}, n::Int, dimension::Symbol, key::String) = dim[:prop][dimension][coord(dim,n,dimension)][key]

dim_prop(data::Dict{String,Any}, args...) = dim_prop(dim(data), args...)
dim_prop(pm::_PM.AbstractPowerModel, args...) = dim_prop(dim(pm), args...)

"""
    dim_meta(dim::Dict{Symbol,Any}[, dimension[, key]])
    dim_meta(data::Dict{String,Any}[, dimension[, key]])
    dim_meta(pm::PowerModels.AbstractPowerModel[, dimension[, key]])

Metadata associated to `dimension`.
"""
function dim_meta end
dim_meta(dim::Dict{Symbol,Any}) = dim[:meta]
dim_meta(dim::Dict{Symbol,Any}, dimension::Symbol) = dim[:meta][dimension]
dim_meta(dim::Dict{Symbol,Any}, dimension::Symbol, key::String) = dim[:meta][dimension][key]
dim_meta(data::Dict{String,Any}, args...) = dim_meta(dim(data), args...)
dim_meta(pm::_PM.AbstractPowerModel, args...) = dim_meta(dim(pm), args...)

"""
    dim_length(dim::Dict{Symbol,Any}[, dimension])
    dim_length(data::Dict{String,Any}[, dimension])
    dim_length(pm::PowerModels.AbstractPowerModel[, dimension])

Return the number of networks or, if `dimension` is specified, return its size.
"""
function dim_length end
dim_length(dim::Dict{Symbol,Any}) = length(dim[:li])
dim_length(dim::Dict{Symbol,Any}, dimension::Symbol) = size(dim[:li], dim[:pos][dimension])
dim_length(data::Dict{String,Any}, args...) = dim_length(dim(data), args...)
dim_length(pm::_PM.AbstractPowerModel, args...) = dim_length(dim(pm), args...)


text:
# Autogenerated wrapper script for LoopTools_jll for x86_64-apple-darwin-libgfortran5
export libooptools

using CompilerSupportLibraries_jll
JLLWrappers.@generate_wrapper_header("LoopTools")
JLLWrappers.@declare_library_product(libooptools, "@rpath/libooptools.dylib")
function __init__()
    JLLWrappers.@generate_init_header(CompilerSupportLibraries_jll)
    JLLWrappers.@init_library_product(
        libooptools,
        "lib/libooptools.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using CubicEoS
using CubicEoSDatabase: Data

C₁C₅ = load(
    BrusilovskyEoSMixture,
    names = ("methane", "n-pentane"),
    component_dbs = (Data.martinez(), Data.brusilovsky_comp()),
    mix_eos_db = Data.brusilovsky_mix()
)

# [0.547413, 0.452587] mole fractions from (Mikyska, 2013, Example 2)
# V = 1e-6 from old vtflash
# (5e3 * V) is (5 kmol m⁻³ * V) just an ordinary value
# T = 370 Kelvins at 5 kmol m⁻³ should correspond to two-phase
# T = 440 Kelvins at 5 kmol m⁻³ should correspond to single-phase

V = 1e-6
χ = [0.547413, 0.452587]

for T in 300:1:450
    for Ση in 100:10:15000
        RT = T * CubicEoS.GAS_CONSTANT_SI
        N = (Ση * V) .* χ
        converged = false
        singlephase = false
        try
            state = vt_flash(C₁C₅, N, V, RT)
            converged = state.converged
            singlephase = state.singlephase
        catch e
            @warn "VTFlash not converged" T Ση e
        end
        print(join([T, Ση], '\t'))
        print('\t', join([converged, singlephase], '\t'))
        println()
    end
end

max_stars_repo_name: vvpisarev/CubicEoS

max_issues_repo_name: vvpisarev/CubicEoS

max_forks_repo_name: vvpisarev/CubicEoS


text:
using Plots
using Statistics
using ProgressBars
L=100
function BoundaryCondition(i)
    if i==L+1
        return 1
    elseif i==0
        return L
    else
        return i
    end
end
function height_cal(i,arr) #calculating the height of the layer in a column
    i=BoundaryCondition(i)
    row=L #starting from the top layer
    height=L
    for j in 1:L
        if arr[row,i]==0.0 #when the entry is empty
            height-=1 #decreasing height by one
            row-=1
        else
            break #exit the loop for the next column
        end
    end
    return height
end
function neighbor_checking(arr, col)
    h=[]
    for i in (col-1):1:(col+1)
        height= height_cal(i,arr)
        if i==col
            push!(h, (height+1))
        else
            push!(h, height)
        end
    end
    return maximum(h)
end
function mean_calculater(arr,L)
    h=[]
    for i in 1:L
        i=BoundaryCondition(i)
        height= height_cal(i,arr)
        push!(h, height)
    end
    return mean(h)
end
function deposing(L,n)
    arr=zeros((L, L))
    meanList=[]
    color=1
    finalP=n
    anim = Animation()
    for particle in ProgressBar(1:n)
        col=rand(1:L)
        height=neighbor_checking(arr, col)
        if height>=L
            finalP=particle
            break
        end
        arr[height, col]=color
        if particle%(10*L*color)==0
            color+=1
        end
        #meanNum=mean_calculater(arr,L)
        #push!(meanList, meanNum)
        p = heatmap(hcat(arr), c=cgrad(:roma, 10, categorical = true, scale = :exp), xlabel="L")
        frame(anim, p)
    end
    gif(anim,"C:\\Users\\Narges\\Documents\\GitHub\\computational_physics\\chapter3\\Fig\\3.3.gif",fps=70)
    return arr, meanList, finalP
end
arr,meanList, finalP=deposing(L,15000)
heatmap(hcat(arr), c=cgrad(:roma, 10, categorical = true, scale = :exp), xlabel="L")
savefig("C:\\Users\\Narges\\Documents\\GitHub\\computational_physics\\chapter3\\Fig\\3.3_1.png")
scatter(1:finalP-1,meanList, xlabel="time", ylabel="mean height of the layer")
savefig("C:\\Users\\Narges\\Documents\\GitHub\\computational_physics\\chapter3\\Fig\\3.3_2.png")


text:
##--- Learnable Transforms ---##

type CenteringTransformation{T<:Real}
    shift::ShiftTransform{T}
end
CenteringTransformation() = CenteringTransformation(ShiftTransform())
function CenteringTransformation(x::AbstractArray)
    tfm = CenteringTransformation()
    learn!(tfm, x)
end
learn!(tfm::CenteringTransformation, x) = (tfm.shift = mean(x); tfm)
transform!(tfm::CenteringTransformation, x) = (transform!(tfm.shift, x); x)
invert!(tfm::CenteringTransformation, x) = (invert!(tfm.shift, x); x)
## What should this return ??
# Base.inv(tfm::CenteringTransformation) = ???

type StandardizeTransformation{T<:Real}
    shift::ShiftTransformation{T}
    scale::ScaleTransformation{T}
end
function learn!{T}(tfm::StandardizeTransformation{T}, x)
    tfm.shift = mean(x)
    tfm.scale = one(T)/std(x)
end
function transform!(tfm::StandardizeTransformation, x)
    transform!(tfm.shift, x)
    transform!(tfm.scale, x)
end
function invert!(tfm::StandardizeTransformation, x)
    invert!(tfm.scale, x)
    invert!(tfm.shift, x)
end
## What should this return ??
# Base.inv(tfm::StandardizeTransformation) = ???


text:
using Test, SimplePadics, SimplePolynomials

x = Padic{5}(-1)

@test x * x == 1
@test x == -1

y = sqrt(x)
@test y * y == -1

@test x + x == 2x
@test x + 1 == 0
@test x - 1 == -Padic{5}(2)
@test valuation(x) == 0
@test valuation(x / 5) == -1
@test valuation(5x) == 1
@test base(x) == 5
@test x // 5 == x / 5
@test x / x == 1
@test x // x == Padic{5}(1)

y = exp(5x)
@test log(y) == -5

@test Padic{5}(2 // 3) * 3 == 2

s = string(x)
@test s[1] == '…'
@test s[end-5:end] == ".0_{5}"

a = Padic{17}(1000)
d = digits(a)
@test d[1] + 17 * d[2] + 17^2 * d[3] == 1000

set_precision(20)
@test get_precision() == 20

a = Padic{5}(100)
@test abs(a) == 1 / 25
@test abs(a) == abs(-a)
@test abs(0a) == 0

x = getx()
F = x^2 + 1
@test has_p_root(F, 5)[1]
t = p_root(F, 5)
@test F(t) == 0


text:
using SmartMeterAnalytics
using Base.Test

@test tou(DateTime(2012, 10, 10, 12)) == :OnPeak # Cooling weekday on-peak
@test tou(DateTime(2012, 10, 10, 8)) == :MidPeak # Cooling weekday mid-peak
@test tou(DateTime(2012, 10, 10, 2)) == :OffPeak # Cooling weekday off-peak

@test tou(DateTime(2012, 1, 10, 9)) == :OnPeak # Heating weekday on-peak
@test tou(DateTime(2012, 1, 10, 12)) == :MidPeak# Heating weekday mid-peak
@test tou(DateTime(2012, 1, 10, 22)) == :OffPeak # Heating weekday off-peak

@test tou(DateTime(2012, 10, 13, 12)) == :OffPeak # Cooling weekend on-peak
@test tou(DateTime(2012, 10, 13, 8)) == :OffPeak # Cooling weekend mid-peak
@test tou(DateTime(2012, 10, 13, 2)) == :OffPeak # Cooling weekend off-peak

@test tou(DateTime(2012, 1, 14, 9)) == :OffPeak # Heating weekend on-peak
@test tou(DateTime(2012, 1, 14, 12)) == :OffPeak# Heating weekend mid-peak
@test tou(DateTime(2012, 1, 14, 22)) == :OffPeak # Heating weekend off-peak

@test tou(DateTime(2012, 7, 1, 12)) == :OffPeak # Cooling holiday on-peak
@test tou(DateTime(2012, 7, 1, 8)) == :OffPeak # Cooling holiday mid-peak
@test tou(DateTime(2012, 7, 1, 2)) == :OffPeak # Cooling holiday off-peak

@test tou(DateTime(2012, 2, 20, 9)) == :OffPeak # Heating holiday on-peak
@test tou(DateTime(2012, 2, 20, 12)) == :OffPeak# Heating holiday mid-peak
@test tou(DateTime(2012, 2, 20, 22)) == :OffPeak # Heating holiday off-peak

max_stars_repo_name: GordStephen/SmartMeterAnalytics

max_issues_repo_name: GordStephen/SmartMeterAnalytics

max_forks_repo_name: GordStephen/SmartMeterAnalytics


text:
"""
BlockOptBackend

Performs the iteration of Algorithm 7.1 after being initiated by
a `Simulation` instance.
"""
mutable struct BlockOptBackend
    model::Model
    driver::Driver
    n::Int
    w::Int
    xₖ::Vector{Float64}

    # following fields storage is initilized in construction
    ∇fₖ::Vector
    hₖ::Vector
    pₖ::Vector
    xₜ::Vector{Float64}
    fₖ::Float64
    fₜ::Float64
    ρ::Float64
    Δₖ::Float64
    ∇fₖ_norm::Float64
    pₖ_norm::Float64
    Sₖ::Matrix{Float64}
    Yₖ::Matrix{Float64}
    Uₖ::Matrix{Float64}
    Vₖ::Matrix{Float64}
    Hₖ::Symmetric{Float64,Matrix{Float64}}
    Qₖ::Matrix{Float64}
    P::Symmetric{Float64,Matrix{Float64}}
    b::Vector{Float64}
    C::Symmetric{Float64,Matrix{Float64}}
    aₖ::Vector{Float64}
    qₖ::Vector{Float64}

    function BlockOptBackend(m, d)
        # TODO: the driver & model must be locked from modifications
        # during simulation, or copy data. Look into a julia Mutex
        n, w = dimension(m), Int(samples(d) / 2)

        xₖ = copy(initial_iterate(m))

        ∇fₖ = hₖ = pₖ = xₜ = qₖ = similar(xₖ)

        fₖ = fₜ = ρ = Δₖ = ∇fₖ_norm = pₖ_norm = 0.0

        Sₖ = orth(randn(n, 2w - 1))

        Yₖ = similar(Sₖ)

        Uₖ = zeros(n, 2w)

        Vₖ = similar(Uₖ)

        Hₖ = Symmetric(zeros(n, n) + I)

        Qₖ = zeros(n, 2w + 1)

        P = Symmetric(zeros(2w + 1, 2w + 1))

        b = zeros(2w + 1)

        C = similar(P)

        aₖ = similar(b)

        new(
            m,
            d,
            n,
            w,
            xₖ,
            ∇fₖ,
            hₖ,
            pₖ,
            xₜ,
            fₖ,
            fₜ,
            ρ,
            Δₖ,
            ∇fₖ_norm,
            pₖ_norm,
            Sₖ,
            Yₖ,
            Uₖ,
            Vₖ,
            Hₖ,
            Qₖ,
            P,
            b,
            C,
            aₖ,
            qₖ,
        )
    end
end

# Accessors used in simulation.jl to weave! (TODO: Restrict this class as well)
model(b::BlockOptBackend) = getfield(b, :model)

driver(b::BlockOptBackend) = getfield(b, :driver)

fₖ(b::BlockOptBackend) = getfield(b, :fₖ)

∇fₖ_norm(b::BlockOptBackend) = getfield(b, :∇fₖ_norm)

pₖ_norm(b::BlockOptBackend) = getfield(b, :pₖ_norm)

Δₖ(b::BlockOptBackend) = getfield(b, :Δₖ)

ρ(b::BlockOptBackend) = getfield(b, :ρ)

# Forwarded Model Methods
obj(b::BlockOptBackend, x) = obj(model(b), x)

grad(b::BlockOptBackend, x) = grad(model(b), x)

grad!(b::BlockOptBackend, out, x) = grad!(model(b), out, x)

dimension(b::BlockOptBackend) = dimension(model(b))

# Forwarded Driver/DriverOptions Methods
QN_update(b::BlockOptBackend) = QN_update(driver(b))

S_update(b::BlockOptBackend) = S_update(driver(b))

pflag(b::BlockOptBackend) = pflag(driver(b))

samples(b::BlockOptBackend) = samples(driver(b))

Δ_max(b::BlockOptBackend) = Δ_max(driver(b))

δ_tol(b::BlockOptBackend) = δ_tol(driver(b))

ϵ_tol(b::BlockOptBackend) = ϵ_tol(driver(b))

max_iterations(b::BlockOptBackend) = max_iterations(driver(b))

"""
terminal(b::BlockOptBackend, k::Int)

The number of trust-region subproblem solves is given by `k`.

## Development:
Remove BlockOptBackend terminal method for one manged in
simulation.jl. Alowing terminal access to trs_counter and
can observe other traced aspects (giving more state options).
"""
function terminal(b::BlockOptBackend, k::Int)

    if ∇fₖ_norm(b) < ϵ_tol(b) || max_iterations(b) ≤ k

        return true
    end

    return false
end

function secantQN(b::BlockOptBackend)

    b.Hₖ = QN_update(b)(b.Hₖ, b.pₖ, grad(b, b.xₜ) - b.∇fₖ, δ_tol(b))

    return nothing
end

function blockQN(b::BlockOptBackend)

    b.Hₖ = QN_update(b)(b.Hₖ, b.Uₖ, b.Vₖ, δ_tol(b))

    return nothing
end

function update_Sₖ(b::BlockOptBackend)

    b.Sₖ = S_update(b)(b.Sₖ, b.Yₖ, b.pₖ)

    return nothing
end

function build_trs(b::BlockOptBackend)

    b.Qₖ = orth([b.hₖ / b.∇fₖ_norm b.∇fₖ / b.∇fₖ_norm b.Yₖ])

    b.P = Symmetric(b.Qₖ' * b.Hₖ * b.Qₖ)

    b.b = b.Qₖ' * b.Hₖ * b.∇fₖ

    b.C = Symmetric(b.Qₖ' * b.Hₖ * b.Hₖ * b.Qₖ)

    return nothing
end

function solve_trs(b::BlockOptBackend)

    aₖ, info = trs_small(b.P, b.b, b.Δₖ, b.C, compute_local = false)

    b.aₖ = aₖ[:, 1]

    b.qₖ = b.Qₖ * b.aₖ

    b.pₖ = b.Hₖ * b.qₖ

    b.pₖ_norm = norm(b.pₖ)

    return nothing
end

function build_trial(b::BlockOptBackend)

    b.xₜ = b.xₖ + b.pₖ

    b.fₜ = obj(b, b.xₜ)

    b.ρ = (b.fₖ - b.fₜ) / (0 - (0.5 * dot(b.qₖ, b.Hₖ, b.qₖ) + (b.Hₖ * b.∇fₖ)' * b.qₖ))

    return nothing
end

function update_Δₖ(b::BlockOptBackend)

    if b.ρ < 0.25

        b.Δₖ = 0.25 * b.Δₖ

    elseif b.ρ > 0.75 && b.pₖ_norm ≈ b.Δₖ

        b.Δₖ = min(2 * b.Δₖ, Δ_max(b))

    end

    return nothing
end

function accept_trial(b::BlockOptBackend)
    if b.ρ > 0

        b.xₖ = b.xₜ

        b.fₖ = b.fₜ

        return true
    end

    return false
end

function gAD(b::BlockOptBackend, S)

    Sdual = Dual{Float64}.(b.xₖ, eachcol(S)...)

    Ydual = Dual{Float64}.(similar(b.xₖ), eachcol(S)...)

    Ydual = grad!(b, Ydual, Sdual)

    g = similar(b.xₖ)

    Y = similar(S)

    @views for i = 1:length(g)
        g[i] = Ydual[i].value

        Y[i, :] = Ydual[i].partials[:]
    end

    return g, Y
end

function build_UV(b::BlockOptBackend)

    b.Uₖ = [b.Sₖ b.∇fₖ / b.∇fₖ_norm]

    b.Vₖ = [b.Yₖ b.hₖ / b.∇fₖ_norm]

    return nothing
end

function gHS(b::BlockOptBackend)

    w = b.w

    g, Y₁ = gAD(b, b.Sₖ[:, 1:w])

    _, Y₂ = gAD(b, [b.Sₖ[:, w+1:end] g])

    b.∇fₖ_norm, b.∇fₖ, b.hₖ, b.Yₖ = norm(g), g, Y₂[:, end], [Y₁ Y₂[:, 1:end-1]]

    build_UV(b)

    return nothing
end

function initialize(b::BlockOptBackend)

    b.fₖ = obj(b, b.xₖ)

    gHS(b)

    α = mean(eigvals(b.Sₖ' * b.Yₖ))

    b.Δₖ = min(1.1 * b.∇fₖ_norm / (2 * α), Δ_max(b))

    b.Hₖ = 1 / α * b.Hₖ

    blockQN(b)

    return nothing
end


text:
using RCall

function load_R_library(libname)
    reval("library($libname)")
end


text:
module BenchmarkingInterpolationKernels

using BenchmarkTools, Statistics, Printf
using InterpolationKernels
using InterpolationKernels:
    compute_weights, compute_offset_and_weights,
    brief, support, infimum, supremum

function compute_weights!(y::Array{T,2},
                          ker::Kernel{T,1},
                          x::Array{T,1}) where {T}
    @assert size(y) == (1,length(x))
    @inbounds @simd for i in eachindex(x)
        wgts = compute_weights(ker, x[i])
        y[1,i] = wgts[1]
    end
    nothing
end

function compute_weights!(y::Array{T,2},
                          ker::Kernel{T,2},
                          x::Array{T,1}) where {T}
    @assert size(y) == (2,length(x))
    @inbounds @simd for i in eachindex(x)
        wgts = compute_weights(ker, x[i])
        y[1,i] = wgts[1]
        y[2,i] = wgts[2]
    end
    nothing
end

function compute_weights!(y::Array{T,2},
                          ker::Kernel{T,3},
                          x::Array{T,1}) where {T}
    @assert size(y) == (3,length(x))
    @inbounds @simd for i in eachindex(x)
        wgts = compute_weights(ker, x[i])
        y[1,i] = wgts[1]
        y[2,i] = wgts[2]
        y[3,i] = wgts[3]
    end
    nothing
end

function compute_weights!(y::Array{T,2},
                          ker::Kernel{T,4},
                          x::Array{T,1}) where {T}
    @assert size(y) == (4,length(x))
    @inbounds @simd for i in eachindex(x)
        wgts = compute_weights(ker, x[i])
        y[1,i] = wgts[1]
        y[2,i] = wgts[2]
        y[3,i] = wgts[3]
        y[4,i] = wgts[4]
    end
    nothing
end

function compute_offset_and_weights!(y::Array{T,2},
                                     ker::Kernel{T,1},
                                     x::Array{T,1}) where {T}
    @assert size(y) == (2,length(x))
    @inbounds @simd for i in eachindex(x)
        off, wgts = compute_offset_and_weights(ker, x[i])
        y[1,i] = off
        y[2,i] = wgts[1]
    end
    nothing
end

function compute_offset_and_weights!(y::Array{T,2},
                                     ker::Kernel{T,2},
                                     x::Array{T,1}) where {T}
    @assert size(y) == (3,length(x))
    @inbounds @simd for i in eachindex(x)
        off, wgts = compute_offset_and_weights(ker, x[i])
        y[1,i] = off
        y[2,i] = wgts[1]
        y[3,i] = wgts[2]
    end
    nothing
end

function compute_offset_and_weights!(y::Array{T,2},
                                     ker::Kernel{T,3},
                                     x::Array{T,1}) where {T}
    @assert size(y) == (4,length(x))
    @inbounds @simd for i in eachindex(x)
        off, wgts = compute_offset_and_weights(ker, x[i])
        y[1,i] = off
        y[2,i] = wgts[1]
        y[3,i] = wgts[2]
        y[4,i] = wgts[3]
    end
    nothing
end

function compute_offset_and_weights!(y::Array{T,2},
                                     ker::Kernel{T,4},
                                     x::Array{T,1}) where {T}
    @assert size(y) == (5,length(x))
    @inbounds @simd for i in eachindex(x)
        off, wgts = compute_offset_and_weights(ker, x[i])
        y[1,i] = off
        y[2,i] = wgts[1]
        y[3,i] = wgts[2]
        y[4,i] = wgts[3]
        y[5,i] = wgts[4]
    end
    nothing
end

function simdmap!(dst::AbstractArray{<:Any,N},
                  f,
                  src::AbstractArray{<:Any,N}) where {N}
    @inbounds @simd for i in eachindex(dst, src)
        dst[i] = f(src[i])
    end
    return dst
end

number_of_operations(f, ker::Kernel) = number_of_operations(f, typeof(ker))
number_of_operations(ker::Kernel) = number_of_operations(ker, ker)

# For calling the function on a single value, the number of operations is the
# maximum possible.  In case a constant is returned, the number of operations is
# assumed to be one (not zero).  Function compute_offset_and_weights usually adds 3
# operations compared to compute_weights.
for (K, n1, n2, n3) in ((BSpline{1},               2,  2,  2),
                        (BSplinePrime{1},          1,  1,  2), # n1 should be 0
                        (BSpline{2},               2,  1,  4),
                        (BSplinePrime{2},          2,  1,  4),
                        (BSpline{3},               3,  6,  9),
                        (BSplinePrime{3},          2,  3,  6),
                        (BSpline{4},               6, 15, 18),
                        (BSplinePrime{4},          5, 13, 16),
                        (CubicSpline,              6, 15, 18),
                        (CubicSplinePrime,         6, 11, 14),
                        (CardinalCubicSpline,      6, 11, 14),
                        (CardinalCubicSplinePrime, 6, 24, 27), # weights not optimized
                        (CatmullRomSpline,         6, 10, 13),
                        (CatmullRomSplinePrime,    5, 11, 14),)
    @eval begin
        number_of_operations(::$K, ::Type{<:$K}) = $n1
        number_of_operations(::typeof(compute_weights),::Type{<:$K}) = $n2
        number_of_operations(::typeof(compute_offset_and_weights),::Type{<:$K}) = $n3
    end
end

gflops(nops::Real, ns::Real) = (ns, nops/ns)

function runtests(n::Int = 1000; table::Bool=false, verbose::Bool=false, io::IO=stdout)
    # Benchmark calling kernels.
    if table
        print(io, "Benchmark calling kernels ",
              "(minimum times in nanoseconds for $n calls):\n")
        print(io, "| Julia version                | $VERSION |\n")
        print(io, "|:-----------------------------|------:|\n")
    end
    for T in (Float32, Float64)
        for ker in (BSpline{1,T}(),
                    BSpline{2,T}(),
                    BSpline{3,T}(),
                    BSpline{4,T}(),
                    CubicSpline{T}(-1/2,1/6),
                    CatmullRomSpline{T}(),
                    CardinalCubicSpline{T}(-1/2),
                    MitchellNetravaliSpline{T}(1/3,1/3),)
            xmin = infimum(support(ker))
            xmax = supremum(support(ker))
            x = (xmax - xmin)*rand(T, n) .+ xmin;
            @assert eltype(x) === T
            y = Array{T}(undef, n)
            nops = number_of_operations(ker)*n
            if !table
                print(io, "Testing calling ", summary(ker), " ", n, " times: ")
                flush(io)
            end
            b = @benchmark $(simdmap!)($y, $ker, $x)
            t = b.times
            if table
                @printf(io, "| %-28s | %5.0f |\n", summary(ker), minimum(t))
            else
                @printf(io, "%.3f Gflops\n", nops/minimum(t))
                if verbose
                    @printf(io, "  - memory: %d allocations, %d bytes\n",
                            b.allocs, b.memory)
                    @printf(io, "  - minimum time: %7.0f ns, %.3f Gflops\n",
                            gflops(nops, minimum(t))...)
                    @printf(io, "  - median time:  %7.0f ns, %.3f Gflops\n",
                            gflops(nops, median(t))...)
                    @printf(io, "  - mean time:    %7.0f ns, %.3f Gflops\n",
                            gflops(nops, mean(t))...)
                    @printf(io, "  - maximum time: %7.0f ns, %.3f Gflops\n",
                            gflops(nops, maximum(t))...)
                    println(io)
                end
            end
        end
    end

    # Benchmark compute_weights.
    if table
        print(io, "\nBenchmarking `compute_weights` ",
              "(minimum times in nanoseconds for $n calls):\n")
        print(io, "| Julia version                | $VERSION |\n")
        print(io, "|:-----------------------------|------:|\n")
    end
    for T in (Float32, Float64)
        for ker in (BSpline{1,T}(),
                    BSpline{2,T}(),
                    BSpline{3,T}(),
                    BSpline{4,T}(),
                    CubicSpline{T}(-1/2,1/6),
                    CatmullRomSpline{T}(),
                    CardinalCubicSpline{T}(-1/2),
                    MitchellNetravaliSpline{T}(1/3,1/3),)
            if isodd(length(ker))
                x = rand(T, n) .- one(T)/2
            else
                x = rand(T, n)
            end
            @assert eltype(x) === T
            z = Array{T}(undef, length(ker), n)
            nops = number_of_operations(compute_weights, ker)*n
            if !table
                print(io, "Testing `compute_weights` for ",
                      summary(ker), " ", n, " times: ")
                flush(io)
            end
            b = @benchmark $(compute_weights!)($z, $ker, $x)
            t = b.times
            if table
                @printf(io, "| %-28s | %5.0f |\n", summary(ker), minimum(t))
            else
                @printf(io, "%.3f Gflops\n", nops/minimum(t))
                if verbose
                    @printf(io, "  - memory: %d allocations, %d bytes\n",
                            b.allocs, b.memory)
                    @printf(io, "  - minimum time: %7.0f ns, %.3f Gflops\n",
                            gflops(nops, minimum(t))...)
                    @printf(io, "  - median time:  %7.0f ns, %.3f Gflops\n",
                            gflops(nops, median(t))...)
                    @printf(io, "  - mean time:    %7.0f ns, %.3f Gflops\n",
                            gflops(nops, mean(t))...)
                    @printf(io, "  - maximum time: %7.0f ns, %.3f Gflops\n",
                            gflops(nops, maximum(t))...)
                    println(io)
                end
            end
        end
    end

    # Benchmark compute_offset_and_weights.
    if table
        print(io, "\nBenchmarking `compute_offset_and_weights` ",
              "(minimum times in nanoseconds for $n calls):\n")
        print(io, "| Julia version                | $VERSION |\n")
        print(io, "|:-----------------------------|------:|\n")
    end
    for T in (Float32, Float64)
        x = 200*rand(T, n) .- 100;
        for ker in (BSpline{1,T}(),
                    BSpline{2,T}(),
                    BSpline{3,T}(),
                    BSpline{4,T}(),
                    CubicSpline{T}(-1/2,1/6),
                    CatmullRomSpline{T}(),
                    CardinalCubicSpline{T}(-1/2),
                    MitchellNetravaliSpline{T}(1/3,1/3),)
            nops = number_of_operations(compute_offset_and_weights, ker)*n
            if !table
                print(io, "Testing `compute_offset_and_weights` for ",
                      summary(ker), " ", n, " times: ")
                flush(io)
            end
            z = Array{T}(undef, 1 + length(ker), n)
            b = @benchmark $(compute_offset_and_weights!)($z, $ker, $x)
            t = b.times
            if table
                @printf(io, "| %-28s | %5.0f |\n", summary(ker), minimum(t))
            else
                @printf(io, "%.3f Gflops\n", nops/minimum(t))
                if verbose
                    @printf(io, "  - memory: %d allocations, %d bytes\n",
                            b.allocs, b.memory)
                    @printf(io, "  - minimum time: %7.0f ns, %.3f Gflops\n",
                            gflops(nops, minimum(t))...)
                    @printf(io, "  - median time:  %7.0f ns, %.3f Gflops\n",
                            gflops(nops, median(t))...)
                    @printf(io, "  - mean time:    %7.0f ns, %.3f Gflops\n",
                            gflops(nops, mean(t))...)
                    @printf(io, "  - maximum time: %7.0f ns, %.3f Gflops\n",
                            gflops(nops, maximum(t))...)
                    println(io)
                end
            end
        end
    end
end

end # module

if !isinteractive()
    let verbose = false,
        table=false,
        n = 1000

        option = true
        for i in 1:length(ARGS)
            arg = ARGS[i]
            if option && (arg == "-h" || arg == "--help")
                println("Usage:  julia -O3 $(@__FILE__) [-h|--help] [-v] [--] [NUMBER]")
                println("Options:  julia -O3 SCRIPT [-h|--help] [-v] [--] [NUMBER]")
                println("    -h|--help      Print this help.")
                println("    -v|--verbose   Print detailed results.")
                println("    -t|--table     Print results as a table.")
                println("    NUMBER         Vector size (default: $n")
                exit(0)
            elseif option && (arg == "-v" || arg == "--verbose")
                verbose = true
                continue
            elseif option && (arg == "-t" || arg == "--table")
                table = true
                continue
            elseif option && (arg == "--")
                option = false
                continue
            end
            if i == length(ARGS)
                n = tryparse(Int, ARGS[i], base=10)
                if n === nothing
                    println(stderr, "invalid number of elements")
                    exit(1)
                end
            else
                println(stderr, "too many arguments")
                exit(1)
            end
        end
        BenchmarkingInterpolationKernels.runtests(n, verbose=verbose, table=table)
    end
end
nothing


text:
using EquationEditor

ex = :(x/(y+x)^2)
lst = latexify(ex)
Base.find_package("WebSockets")
using Test

@testset "handlefiles" begin
    include("handlefiles_test.jl")
end
@testset "serve" begin
    include("serve_test.jl")
end

max_stars_repo_name: hustf/EquationEditor

max_issues_repo_name: hustf/EquationEditor

max_forks_repo_name: hustf/EquationEditor