text:
using SparseArrays
using StatsBase
using LinearAlgebra
using Random
using ProgressMeter
using ForwardDiff

mutable struct Param{T}
	A::SparseMatrixCSC{T,Int64}
	At::SparseMatrixCSC{T,Int64}  # for adjoint solver
	dAdx::Vector{SparseMatrixCSC{T,Int64}}  # necessary for the gradient
	x1::Vector{T}
	x2::Vector{T}
	nx::Int64
	nz::Int64
	dx::Float64
	dz::Float64
end

include("operatorA_Neumann.jl")
#include("operatorA_Dirichlet.jl")

function Param(nx::Int,nz::Int,mpars=nothing)

	dz=inv(nz-1);
	dx=inv(nx-1);
	z=(0:nz-1)*dz;
	x=(0:nx-1)*dx;

	return Param([z,x], mpars)
end

function Param(mgrid, mpars=nothing)
	dz=Float64(step(mgrid[1]))
	dx=Float64(step(mgrid[2]))
	nz,nx=length.(mgrid)
		nznx=nz*nx

	if(!(mpars===nothing))
		T=eltype(mpars)
	else
		T=Float64
	end
	A=spzeros(T, nznx+1,nznx)
	At=spzeros(T, nznx,nznx+1)

	dAdx=dAdmpars(dx,dz,nx,nz,T)

	pa=Param(A,At,dAdx,zeros(T, nznx),zeros(T, nznx+1),nx,nz,dx,dz)

	if(!(mpars===nothing))
		updateA!(pa, mpars)
	end

	return pa
end

"""
"""
function applyinvA!(fout, fin, pa; A=pa.A, mute_boundary_source=true)
	nx=pa.nx
	nz=pa.nz
	nznx=nz*nx

	x1=pa.x1
	x2=pa.x2
	fill!(pa.x1, 0.0)
	fill!(pa.x2, 0.0)

	for i in 1:nznx
		x2[i]=fin[i]
	end

	if(mute_boundary_source)
		for j = 2 : nx - 1 # interior x (column) loop:
		    k = 1 + nz*(j - 1);# initial boundary z for this x (k = i + m*(j - 1))
		    for i = [1 nz]  # % north & south boundary z (row) loop:
			x2[k] = 0; # 0 to equate boundary values in source term ... (?)
			k = k + nz - 1; #% from north to south boundary
		    end
		end

		for j = [1 nx] # west & east boundary x (column) loop:
		    k = 2 + nz*(j - 1); #% initial interior z for this x (k = i + m*(j - 1))
		    for i = 2 : nz - 1 # interior z (row) loop:
			x2[k]= 0;  #% 0 to equate boundary values
			k = k + 1; #% to next z
		    end
		 end

		for j = [1 nx] #% west & east boundary x (column) loop:
		    k = 1 + nz*(j - 1); #% start on north boundary
		    for i = [1 nz] #% north & south boundary y (row) loop:
			x2[k] = 0; #% 0 to equate boundary values
			k = k + nz - 1; #% from north to south boundary
		    end
		 end
	 end

#	LinearAlgebra.ldiv!(x1,A,x2) # direct inverse solve using backslash operator
	x1 .= A\x2 # memory allocated here?
	#x11 = A\x2 # memory allocated here?

	for i in 1:nznx
		fout[i]=x1[i]
	end

	return fout
end

function applyinvAt!(fout, fin, pa; At=pa.At)
	nx=pa.nx
	nz=pa.nz
	nznx=nz*nx

	x1=pa.x1
	x2=pa.x2
	fill!(pa.x1, 0.0)
	fill!(pa.x2, 0.0)

	for i in 1:nznx
		x1[i]=fin[i]
	end

#	LinearAlgebra.ldiv!(x1,A,x2) # direct inverse solve using backslash operator
	x2 .= At\x1 # memory allocated here?

	for i in 1:nznx+1
		fout[i]=x2[i]
	end

	return fout
end

"""
Apply `A` to `fin` without allocating memory.
"""
function applyA!(fout, fin, pa; A=pa.A)
	fill!(fout, 0.0) # clear output
	nx=pa.nx
	nz=pa.nz
	nznx=nz*nx

	x1=pa.x1
	x2=pa.x2
	fill!(pa.x1, 0.0)
	fill!(pa.x2, 0.0)
	for i in 1:nznx
		x1[i]=fin[i]
	end

	mul!(x2, A, x1)
	for i in 1:nznx
		fout[i]=x2[i]
	end
	rmul!(fout, -1.0)
	return fout
end

function applyinvA(fin,pa;A=pa.A)
	fout=zeros(size(A,2))
	applyinvA!(fout,fin,pa,A=A)
	return fout
end
function applyinvAt(fin,pa;At=pa.At)
	fout=zeros(size(At,2))
	applyinvAt!(fout,fin,pa,At=At)
	return fout
end

function applyA(fin,pa;A=pa.A)
	fout=zeros(size(fin))
	applyA!(fout,fin,pa,A=A)
	return fout
end


text:
#### plot with range

using GPlot

#############################################################
######################### should be in the ~3 seconds range.
t = @elapsed begin
    tempdir = mktempdir()
    f = Figure("___")
    x, y = [1, 2], [3, 4]
    x2 = range(1, stop=2, length=10)
    y2 = exp.(x2)
    # -- plot
    plot(x, y); plot(x, y, ls="--", lw=0.1, color="blue", label="label1")
    plot!(x2, y2, lwidth=0.02)
    xlabel("xlabel"); ylabel("ylabel"); title("title")
    legend(pos="tl", fontsize=7)
    xticks([1, 2], ["a", "b"]); yticks([1, 2], ["a", "b"])
    xscale("lin"); yscale("lin")
    xlim(min=0,max=5); ylim(min=0,max=5)
    GPlot.assemble_figure(f)
    savefig(f, res=100, format="png", path=tempdir)
    # -- hist
    GPlot.destroy(f)
    f = Figure()
    x = rand(10)
    hist(x, fill="cornflowerblue", color="white", scaling="pdf", nbins=50)
    bar(x, x, 2x; stacked=true, fills=["blue","red"], colors=["red","blue"])
    GPlot.assemble_figure(f)
    # if there's a working GLE backend, warm that up too
    savefig(f, format="pdf", path=tempdir)
    # clear everything
    GPlot.destroy(f)
end; println(".............done in $(round(t, digits=1))s ✅")
#############################################################
#############################################################

text!("hello", (3, -2))

text!("goodbye", (5, 0), font="tt")

text!("another-one", (5, -1), fontsize=0.8)

text!("blah", (3, 0))

clo!()

# so it's preserved, not good

debug_gle(gcf())

################################
################################

continuous_preview(true)

clf()

θ = range(-4pi, 4pi, length=200);
z = range(-2, 2, length=200);
r = z.^2 .+ 1;
x = r .* sin.(θ);
y = r .* cos.(θ);

xx = 4 * ones(10)
yy = range(-4,4,length=10)
zz = 0*yy .- 2

xx2 = range(-4,4,length=10)
yy2 = 4 * ones(10)
zz2 = 0*yy2 .- 2

#cla()

plot3(x,y,z, marker=".", msize=.2, mcol="red")
# plot3!(xx,yy,zz)
# plot3!(xx2,yy2,zz2)

plot3(x,y,z, marker=".", msize=0.01, mcol="salmon")

plot(y,z)

preview()

plot(x, y)

ext: jl


text:
"""
    random_mps(ElT::Type{<:Number},sites::Vector{<: Index},χ::Int64,σ::Float64)

Make a random MPS with bond dimension `χ` using the Hilbert spacee `sites`.
Each bulk tensor has one site index, and two link indices. The components of each
tensor, with type `ElT`, are randomly drawn from a uniform distribution centered
around zero, with width `σ`.
"""
function random_mps(ElT::Type{<:Number}, sites::Vector{<:Index}, χ::Int64, σ::Float64)
  d = dim(sites[1]) # Dimension of the local Hilbert space
  N = length(sites)
  links = [Index(χ; tags="Link, l=$l") for l in 1:(N - 1)]
  M = ITensor[]
  if N == 1
    rand_mat = σ * (ones(d) - 2 * rand(d))
    if ElT <: Complex
      rand_mat += im * σ * (ones(d) - 2 * rand(d))
    end
    push!(M, ITensor(rand_mat, sites[1]))
    return MPS(M)
  end
  # Site 1
  rand_mat = σ * (ones(d, χ) - 2 * rand(d, χ))
  if ElT <: Complex
    rand_mat += im * σ * (ones(d, χ) - 2 * rand(d, χ))
  end
  push!(M, ITensor(rand_mat, sites[1], links[1]))
  for j in 2:(N - 1)
    rand_mat = σ * (ones(χ, d, χ) - 2 * rand(χ, d, χ))
    if ElT <: Complex
      rand_mat += im * σ * (ones(χ, d, χ) - 2 * rand(χ, d, χ))
    end
    push!(M, ITensor(rand_mat, links[j - 1], sites[j], links[j]))
  end
  # Site N
  rand_mat = σ * (ones(χ, d) - 2 * rand(χ, d))
  if ElT <: Complex
    rand_mat += im * σ * (ones(χ, d) - 2 * rand(χ, d))
  end
  push!(M, ITensor(rand_mat, links[N - 1], sites[N]))
  return MPS(M)
end

"""
    random_mpo(ElT::Type{<:Number},sites::Vector{<: Index},χ::Int64,σ::Float64;
               processtags::Bool=false)

Make a random MPO with bond dimension `χ` using the Hilbert spacee `sites`.
Each bulk tensor has two site indices, and two link indices. The components of each
tensor, with type `ElT`, are randomly drawn from a uniform distribution centered
around zero, with width `σ`.
If `processtags=true`, add the tag `input` to the bra, and the tag `output`
to the ket.
"""
function random_mpo(
  ElT::Type{<:Number}, sites::Vector{<:Index}, χ::Int64, σ::Float64; processtags::Bool=false
)
  d = dim(sites[1]) # Dimension of the local Hilbert space
  #@show sites
  N = length(sites)
  links = [Index(χ; tags="Link, l=$l") for l in 1:(N - 1)]

  M = ITensor[]
  if N == 1
    rand_mat = σ * (ones(d, d) - 2 * rand(d, d))
    if ElT <: Complex
      rand_mat += im * σ * (ones(d, d) - 2 * rand(d, d))
    end
    push!(M, ITensor(rand_mat, sites[1]', sites[1]))
    return MPO(M)
  end

  # Site 1
  rand_mat = σ * (ones(d, χ, d) - 2 * rand(d, χ, d))
  if ElT <: Complex
    rand_mat += im * σ * (ones(d, χ, d) - 2 * rand(d, χ, d))
  end
  push!(M, ITensor(rand_mat, sites[1]', links[1], sites[1]))
  for j in 2:(N - 1)
    rand_mat = σ * (ones(d, χ, d, χ) - 2 * rand(d, χ, d, χ))
    if ElT <: Complex
      rand_mat += im * σ * (ones(d, χ, d, χ) - 2 * rand(d, χ, d, χ))
    end
    push!(M, ITensor(rand_mat, sites[j]', links[j - 1], sites[j], links[j]))
  end
  # Site N
  rand_mat = σ * (ones(d, χ, d) - 2 * rand(d, χ, d))
  if ElT <: Complex
    rand_mat += im * σ * (ones(d, χ, d) - 2 * rand(d, χ, d))
  end
  push!(M, ITensor(rand_mat, sites[N]', links[N - 1], sites[N]))

  U = MPO(M)

  proc_tagcheck = (any(x -> hastags(x, "Input"), sites))
  if proc_tagcheck
    removetags!(U, "Input")
  end
  return U
end

"""
    random_lpdo(ElT::Type{<:Number},sites::Vector{<: Index},χ::Int64,ξ::Int64,σ::Float64)

Make a random LPDO with bond dimension `χ`, kraus dimension `ξ` ,and using
the Hilbert spacee `sites`. Each bulk tensor has one site index, one kraus index,
and two link indices. The components of each tensor, with type `ElT`, are
randomly drawn from a uniform distribution centered around zero, with width `σ`.
"""
function random_lpdo(
  ElT::Type{<:Number},
  sites::Vector{<:Index},
  χ::Int64,
  ξ::Int64,
  σ::Float64;
  purifier_tags=default_purifier_tags,
)
  d = dim(sites[1]) # Dimension of the local Hilbert space
  N = length(sites)
  links = [Index(χ; tags="Link, l=$l") for l in 1:(N - 1)]
  kraus = [Index(ξ; tags=addtags(purifier_tags, "k=$s")) for s in 1:N]

  M = ITensor[]
  if N == 1
    # Site 1
    rand_mat = σ * (ones(d, ξ) - 2 * rand(rng, d, ξ))
    if ElT <: Complex
      rand_mat += im * σ * (ones(d, ξ) - 2 * rand(rng, d, ξ))
    end
    push!(M, ITensor(rand_mat, sites[1], kraus[1]))
    return LPDO(MPO(M))
  end

  # Site 1
  rand_mat = σ * (ones(d, χ, ξ) - 2 * rand(d, χ, ξ))
  if ElT <: Complex
    rand_mat += im * σ * (ones(d, χ, ξ) - 2 * rand(d, χ, ξ))
  end
  push!(M, ITensor(rand_mat, sites[1], links[1], kraus[1]))
  # Site 2..N-1
  for j in 2:(N - 1)
    rand_mat = σ * (ones(d, χ, ξ, χ) - 2 * rand(d, χ, ξ, χ))
    if ElT <: Complex
      rand_mat += im * σ * (ones(d, χ, ξ, χ) - 2 * rand(d, χ, ξ, χ))
    end
    push!(M, ITensor(rand_mat, sites[j], links[j - 1], kraus[j], links[j]))
  end
  # Site N
  rand_mat = σ * (ones(d, χ, ξ) - 2 * rand(d, χ, ξ))
  if ElT <: Complex
    rand_mat += im * σ * (ones(d, χ, ξ) - 2 * rand(d, χ, ξ))
  end
  push!(M, ITensor(rand_mat, sites[N], links[N - 1], kraus[N]))

  return LPDO(MPO(M), purifier_tags)
end

"""
    random_choi(ElT::Type{<:Number},sites::Vector{<: Index},χ::Int64,ξ::Int64,σ::Float64;
                purifier_tags = default_purifier_tags)

Make a random Choi matrix with bond dimension `χ`, kraus dimension `ξ` ,and using
the Hilbert spacee `sites`. Each bulk tensor has two site indices (corresponding
to input and output indices, one kraus index,and two link indices. The components
of each tensor, with type `ElT`, are randomly drawn from a uniform distribution
centered around zero, with width `σ`.
"""
function random_choi(
  ElT::Type{<:Number},
  sites::Vector{<:Index},
  χ::Int64,
  ξ::Int64,
  σ::Float64;
  purifier_tags=default_purifier_tags,
)
  d = dim(sites[1]) # Dimension of the local Hilbert space
  N = length(sites)
  links = [Index(χ; tags="Link, l=$l") for l in 1:(N - 1)]
  kraus = [Index(ξ; tags=addtags(purifier_tags, "k=$s")) for s in 1:N]

  M = ITensor[]
  if N == 1
    # Site 1
    rand_mat = σ * (ones(d, d, ξ) - 2 * rand(d, d, ξ))
    if ElT <: Complex
      rand_mat += im * σ * (ones(d, d, ξ) - 2 * rand(d, d, ξ))
    end
    push!(M, ITensor(rand_mat, sites[1], sites[1]', kraus[1]))
    return LPDO(choitags(MPO(M)), purifier_tags)
  end

  # Site 1
  rand_mat = σ * (ones(d, d, χ, ξ) - 2 * rand(d, d, χ, ξ))
  if ElT <: Complex
    rand_mat += im * σ * (ones(d, d, χ, ξ) - 2 * rand(d, d, χ, ξ))
  end
  push!(M, ITensor(rand_mat, sites[1], sites[1]', links[1], kraus[1]))
  # Site 2..N-1
  for j in 2:(N - 1)
    rand_mat = σ * (ones(d, d, χ, ξ, χ) - 2 * rand(d, d, χ, ξ, χ))
    if ElT <: Complex
      rand_mat += im * σ * (ones(d, d, χ, ξ, χ) - 2 * rand(d, d, χ, ξ, χ))
    end
    push!(M, ITensor(rand_mat, sites[j], sites[j]', links[j - 1], kraus[j], links[j]))
  end
  # Site N
  rand_mat = σ * (ones(d, d, χ, ξ) - 2 * rand(d, d, χ, ξ))
  if ElT <: Complex
    rand_mat += im * σ * (ones(d, d, χ, ξ) - 2 * rand(d, d, χ, ξ))
  end
  push!(M, ITensor(rand_mat, sites[N], sites[N]', links[N - 1], kraus[N]))

  Λ = MPO(M)
  has_inputtags = (any(x -> hastags(x, "Input"), sites))
  has_outputtags = (any(x -> hastags(x, "Output"), sites))
  if !has_inputtags && !has_outputtags
    addtags!(Λ, "Input"; plev=0, tags="Qubit")
    addtags!(Λ, "Output"; plev=1, tags="Qubit")
  elseif !has_outputtags
    for j in 1:N
      replacetags!(Λ, "Input", "Output"; plev=1)
    end
  end
  noprime!(Λ)
  return LPDO(Λ, purifier_tags)
end

"""
    randomstate(N::Int64; kwargs...)

    randomstate(ElT::Type{<: Number}, N::Int64; kwargs...)

Generates a random quantum state of N qubits.

Optionally, specify an element type, such as `randomstate(Float64, 10)` for a random real state (by default it is complex).

# Arguments
  - `N`: number of qubits
  - `mixed`: if false (default), generate a random MPS; if true, generates a random LPDO
  - `alg`: algorithm used for initialization: `"rand"` initializes random tensor elements;
    `"circuit"` initializes with a random quantum circuit (MPS only).
  - `σ`: size of the 0-centered uniform distribution in `alg="rand"`.
  - `χ`: bond dimension of the MPS/LPDO
  - 'ξ`: kraus dimension (LPDO)
  - `normalize`: if true, return normalized state
"""
function randomstate(ElT::Type{<:Number}, N::Int64; kwargs...)
  sites = siteinds("Qubit", N)
  return randomstate(ElT, sites; kwargs...)
end

randomstate(N::Int64; kwargs...) = randomstate(ComplexF64, N; kwargs...)

function randomstate(ElT::Type{<:Number}, sites::Vector{<:Index}; kwargs...)
  mixed::Bool = get(kwargs, :mixed, false)
  lpdo::Bool = get(kwargs, :lpdo, true)
  if !mixed
    return randomstate(ElT, MPS, sites; kwargs...)
  else
    if lpdo
      return randomstate(ElT, LPDO, sites; kwargs...)
    else
      return randomstate(ElT, MPO, sites; kwargs...)
    end
  end
end

randomstate(sites::Vector{<:Index}; kwargs...) = randomstate(ComplexF64, sites; kwargs...)

function randomstate(ElT::Type{<:Number}, T::Type, sites::Vector{<:Index}; kwargs...)
  χ::Int64 = get(kwargs, :χ, 1)
  ξ::Int64 = get(kwargs, :ξ, 1)
  σ::Float64 = get(kwargs, :σ, 0.1)
  purifier_tags = get(kwargs, :purifier_tags, default_purifier_tags)
  alg::String = get(kwargs, :alg, "rand")
  normalize::Bool = get(kwargs, :normalize, false)

  if ξ > 1
    M = random_lpdo(ElT, sites, χ, ξ, σ; purifier_tags=purifier_tags)
  else
    if T == MPS
      # Build MPS by random parameter initialization
      alg == "rand" && (M = random_mps(ElT, sites, χ, σ))
      alg == "circuit" && (M = randomMPS(ElT, sites, χ))
    elseif T == MPO
      error("initialization of random MPO density matrix not yet implemented.")
    elseif T == LPDO
      M = random_lpdo(ElT, sites, χ, ξ, σ; purifier_tags=purifier_tags)
    else
      error("ansatz type not recognized")
    end
  end
  normalize && normalize!(M)
  return M
end

"""
    randomstate(M::Union{MPS,MPO,LPDO}; kwargs...)

Generate a random state with same Hilbert space (i.e. site indices)
of a reference state `M`.
"""
function randomstate(ElT::Type{<:Number}, M::Union{MPS,MPO,LPDO}; kwargs...)
  hM = hilbertspace(M)
  return randomstate(ElT, hM; kwargs...)
end

randomstate(M::Union{MPS,MPO,LPDO}; kwargs...) = randomstate(ComplexF64, M; kwargs...)

"""
    randomprocess(N::Int64; kwargs...)

    randomprocess(ElT::Type{<: Number}, N::Int64; kwargs...)

Generates a random quantum procecss of N qubits.

Optionally choose the element type with calls like `randomprocess(Float64, 10)` (by default it is complex).

# Arguments
  - `N`: number of qubits
  - `mixed`: if false (default), generates a random MPO; if true, generates a random LPDO.
  - `alg`: initialization criteria, set to `"randompars"` (see `randomstate`).
  - `σ`: size of the 0-centered uniform distribution in `alg="rand"`.
  - `χ`: bond dimension of the MPO/LPDO.
  - 'ξ`: kraus dimension (LPDO).
"""
function randomprocess(ElT::Type{<:Number}, N::Int64; kwargs...)
  sites = siteinds("Qubit", N)
  return randomprocess(ElT, sites; kwargs...)
end

randomprocess(N::Int64; kwargs...) = randomprocess(ComplexF64, N; kwargs...)

function randomprocess(ElT::Type{<:Number}, sites::Vector{<:Index}; kwargs...)
  mixed::Bool = get(kwargs, :mixed, false)
  if mixed
    return randomprocess(ElT, LPDO, sites; kwargs...)
  end
  return randomprocess(ElT, MPO, sites; kwargs...)
end

function randomprocess(sites::Vector{<:Index}; kwargs...)
  return randomprocess(ComplexF64, sites; kwargs...)
end

function randomprocess(ElT::Type{<:Number}, T::Type, sites::Vector{<:Index}; kwargs...)
  χ::Int64 = get(kwargs, :χ, 1)
  ξ::Int64 = get(kwargs, :ξ, 1)
  σ::Float64 = get(kwargs, :σ, 0.1)
  purifier_tags = get(kwargs, :purifier_tags, default_purifier_tags)
  alg::String = get(kwargs, :alg, "rand")
  normalize::Bool = get(kwargs, :normalize, false)
  processtags = !(any(x -> hastags(x, "Input"), sites))

  if ξ > 1
    M = random_choi(ElT, sites, χ, ξ, σ; purifier_tags=purifier_tags)
  else
    if T == MPO
      if alg == "rand"
        M = random_mpo(ElT, sites, χ, σ; processtags=processtags)
      else
        error("randomMPO with circuit initialization not implemented yet")
      end
    elseif T == LPDO
      M = random_choi(ElT, sites, χ, ξ, σ; purifier_tags=purifier_tags)
    else
      error("ansatz type not recognized")
    end
  end
  if normalize
    if M isa MPO
      Φ = normalize!(unitary_mpo_to_choi_mps(M))
      Φ = Φ * √2^length(M)
      M = choi_mps_to_unitary_mpo(Φ)
    else
      normalize!(M; localnorm = 2)
    end
  end
  return M
end

function randomprocess(T::Type, sites::Vector{<:Index}; kwargs...)
  return randomprocess(ComplexF64, T, sites; kwargs...)
end

"""
    randomprocess(M::Union{MPS,MPO}; kwargs...)

Generate a random process with same Hilbert space (i.e. input
and output indices)of a reference process `M`.
"""
function randomprocess(ElT::Type{<:Number}, M::Union{MPS,MPO}; kwargs...)
  mixed = get(kwargs, :mixed, false)
  N = length(M)
  s = Index[]
  for j in 1:N
    push!(s, firstind(M[j]; tags="Site", plev=0))
  end
  proc = randomprocess(ElT, s; mixed=mixed, kwargs...)
  return proc
end

randomprocess(M::Union{MPS,MPO}; kwargs...) = randomprocess(ComplexF64, M; kwargs...)

randomprocess(ElT::Type{<:Number}, L::LPDO; kwargs...) = randomprocess(ElT, L.X; kwargs...)

randomprocess(L::LPDO; kwargs...) = randomprocess(ComplexF64, L; kwargs...)


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:2] <= 1)

                     @variable(m, -4.847398418942631 <= q <= 3.5461678827819525)

                     add_NL_constraint(m, :(sigmoid(-0.19745953189550702 + -0.698583212677272*sigmoid(-0.2024495957518173 + -0.1501654471281335*sigmoid(-0.06497322351233636 + -0.4876409540129627*$(x[1]) + -0.026006048535915927*$(x[2])) + 0.7965619538518913*sigmoid(-0.4201092923480956 + 0.6085672502473995*$(x[1]) + -0.465188513706837*$(x[2])) + 0.8656395810315112*sigmoid(-0.9600965095983032 + 0.022609326108050443*$(x[1]) + -0.12508836215483532*$(x[2]))) + 0.07925279857906231*sigmoid(-0.4774715210476401 + -0.18966514690787006*sigmoid(-0.06497322351233636 + -0.4876409540129627*$(x[1]) + -0.026006048535915927*$(x[2])) + 0.1687755600264067*sigmoid(-0.4201092923480956 + 0.6085672502473995*$(x[1]) + -0.465188513706837*$(x[2])) + 0.9882336248352352*sigmoid(-0.9600965095983032 + 0.022609326108050443*$(x[1]) + -0.12508836215483532*$(x[2]))) + -0.8658353105986021*sigmoid(-0.5967554189497073 + 0.7632066186026765*sigmoid(-0.06497322351233636 + -0.4876409540129627*$(x[1]) + -0.026006048535915927*$(x[2])) + -0.767873838532454*sigmoid(-0.4201092923480956 + 0.6085672502473995*$(x[1]) + -0.465188513706837*$(x[2])) + 0.02410968604346264*sigmoid(-0.9600965095983032 + 0.022609326108050443*$(x[1]) + -0.12508836215483532*$(x[2])))) + sigmoid(0.9047173820312868 + 0.09375520961098927*sigmoid(-0.2024495957518173 + -0.1501654471281335*sigmoid(-0.06497322351233636 + -0.4876409540129627*$(x[1]) + -0.026006048535915927*$(x[2])) + 0.7965619538518913*sigmoid(-0.4201092923480956 + 0.6085672502473995*$(x[1]) + -0.465188513706837*$(x[2])) + 0.8656395810315112*sigmoid(-0.9600965095983032 + 0.022609326108050443*$(x[1]) + -0.12508836215483532*$(x[2]))) + -0.9282166334997988*sigmoid(-0.4774715210476401 + -0.18966514690787006*sigmoid(-0.06497322351233636 + -0.4876409540129627*$(x[1]) + -0.026006048535915927*$(x[2])) + 0.1687755600264067*sigmoid(-0.4201092923480956 + 0.6085672502473995*$(x[1]) + -0.465188513706837*$(x[2])) + 0.9882336248352352*sigmoid(-0.9600965095983032 + 0.022609326108050443*$(x[1]) + -0.12508836215483532*$(x[2]))) + 0.3928411969784107*sigmoid(-0.5967554189497073 + 0.7632066186026765*sigmoid(-0.06497322351233636 + -0.4876409540129627*$(x[1]) + -0.026006048535915927*$(x[2])) + -0.767873838532454*sigmoid(-0.4201092923480956 + 0.6085672502473995*$(x[1]) + -0.465188513706837*$(x[2])) + 0.02410968604346264*sigmoid(-0.9600965095983032 + 0.022609326108050443*$(x[1]) + -0.12508836215483532*$(x[2])))) + sigmoid(0.16193188153327842 + 0.4964335766578394*sigmoid(-0.2024495957518173 + -0.1501654471281335*sigmoid(-0.06497322351233636 + -0.4876409540129627*$(x[1]) + -0.026006048535915927*$(x[2])) + 0.7965619538518913*sigmoid(-0.4201092923480956 + 0.6085672502473995*$(x[1]) + -0.465188513706837*$(x[2])) + 0.8656395810315112*sigmoid(-0.9600965095983032 + 0.022609326108050443*$(x[1]) + -0.12508836215483532*$(x[2]))) + -0.6165233713312808*sigmoid(-0.4774715210476401 + -0.18966514690787006*sigmoid(-0.06497322351233636 + -0.4876409540129627*$(x[1]) + -0.026006048535915927*$(x[2])) + 0.1687755600264067*sigmoid(-0.4201092923480956 + 0.6085672502473995*$(x[1]) + -0.465188513706837*$(x[2])) + 0.9882336248352352*sigmoid(-0.9600965095983032 + 0.022609326108050443*$(x[1]) + -0.12508836215483532*$(x[2]))) + -0.4911670927879337*sigmoid(-0.5967554189497073 + 0.7632066186026765*sigmoid(-0.06497322351233636 + -0.4876409540129627*$(x[1]) + -0.026006048535915927*$(x[2])) + -0.767873838532454*sigmoid(-0.4201092923480956 + 0.6085672502473995*$(x[1]) + -0.465188513706837*$(x[2])) + 0.02410968604346264*sigmoid(-0.9600965095983032 + 0.022609326108050443*$(x[1]) + -0.12508836215483532*$(x[2])))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
function logistic_test()
    @testset "Logistic" begin
        n = 15;
        xV = range(-2.0, 3.0, length = n);

        lb = -2.5;
        ub = 1.5;
        lrShifter = 3.0;
        yV = logistic(xV; lb = lb, ub = ub, lrShifter = lrShifter);
        @test isa(yV, Vector{Float64})
        @test size(yV) == size(xV)
        @test all(yV .> lb)  &&  all(yV .< ub)

        glf = GeneralizedLogistic(lb = lb, ub = ub, lrShifter = lrShifter);
        y2V = logistic(glf, xV);
        @test yV ≈ y2V
    end
end

@testset "Logistic" begin
    logistic_test();
end

# ----------


text:
##################################################
## Classical derivative-based, iterative, root-finding algorithms
##
## If ri = f^(i-1)/f^(i), then these have an update step `xn - delta` where:
##
## * Newton: delta = r1  # order 2 if simple root (multiplicity 1)
## * Halley: delta = 2*r2/(2r2 - r1) * r1 # order 3 if simple root
## * Schroder: delta = r2  / (r2 - r1) * r1  # order 2
## * Thukral(3): delta =  (-2*r3)*(r2 - r1)/(r1^2 - 3*r1*r3 + 2*r2*r3) * r1 # order 3
## * Thukral(4): delta =  3*r1*r2*r4*(r1^2 - 3*r1*r3 + 2*r2*r3)/(-r1^3*r2 + 4*r1^2*r2*r4 + 3*r1^2*r3*r4 - 12*r1*r2*r3*r4 + 6*r2^2*r3*r4) # order 4
##
## The latter two come from
## [Thukral](http://article.sapub.org/10.5923.j.ajcam.20170702.05.html). They are not implemented.

## Newton
abstract type AbstractNewtonLikeMethod <: AbstractUnivariateZeroMethod end
fn_argout(::AbstractNewtonLikeMethod) = 2
struct Newton <: AbstractNewtonLikeMethod end
"""

    Roots.Newton()

Implements Newton's [method](http://tinyurl.com/b4d7vls):
`xᵢ₊₁ =  xᵢ - f(xᵢ)/f'(xᵢ)`.  This is a quadratically convergent method requiring
one derivative. Two function calls per step.

Example
```
find_zero((sin,cos), 3.0, Roots.Newton())
```

If function evaluations are expensive one can pass in a function which returns (f, f/f') as follows

```
find_zero(x -> (sin(x), sin(x)/cos(x)), 3.0, Roots.Newton())
```

This can be advantageous if the derivative is easily computed from the
value of f, but otherwise would be expensive to compute.

The error, `eᵢ = xᵢ - α`, can be expressed as `eᵢ₊₁ = f[xᵢ,xᵢ,α]/(2f[xᵢ,xᵢ])eᵢ²` (Sidi, Unified treatment of regula falsi, Newton-Raphson, secant, and Steffensen methods for nonlinear equations).

"""
Newton

# we store x0,x1,fx0,fx1 **and** Δ = fx1/f'(x1)
struct NewtonState{T,S} <: AbstractUnivariateZeroState{T,S}
    xn1::T
    xn0::T
    Δ::T
    fxn1::S
    fxn0::S
end

function init_state(M::Newton, F::Callable_Function, x)
    x₀ = float(first(x))
    fx₀, Δ = F(x₀)
    x₁ = x₀ - Δ
    state = init_state(M, F, x₀, x₁, fx₀, fx₀)
end

# compute fx₁, Δ
function init_state(::Newton, F, x₀, x₁, fx₀, fx₁)
    fx₁, Δ = F(x₁)
    NewtonState(x₁, x₀, Δ, fx₁, fx₀)
end

initial_fncalls(M::Newton) = 2

function update_state(method::Newton, F, o::NewtonState{T,S}, options, l=NullTracks()) where {T, S}

    xn0, xn1 = o.xn0, o.xn1
    fxn0, fxn1 = o.fxn0, o.fxn1
    Δ = o.Δ

    if isissue(Δ)
        log_message(l, "Issue with `f/f′'")
        return o, true
    end

    xn0, xn1 = xn1, xn1-Δ
    fxn0 = fxn1
    fxn1, Δ = F(xn1)
    incfn(l,2)

    @set! o.xn0 = xn0
    @set! o.xn1 = xn1
    @set! o.Δ = Δ
    @set! o.fxn0 = fxn0
    @set! o.fxn1 = fxn1

    return o, false

end

"""
    Roots.newton(f, fp, x0; kwargs...)

Implementation of Newton's method: `xᵢ₊₁ =  xᵢ - f(xᵢ)/f'(xᵢ)`.

Arguments:

* `f::Function` -- function to find zero of

* `fp::Function` -- the derivative of `f`.

* `x0::Number` -- initial guess. For Newton's method this may be complex.

With the `FowardDiff` package derivatives may be computed automatically. For example,  defining
`D(f) = x -> ForwardDiff.derivative(f, float(x))` allows `D(f)` to be used for the first derivative.

Keyword arguments are passed to `find_zero` using the `Roots.Newton()` method.

See also `Roots.newton((f,fp), x0)` and `Roots.newton(fΔf, x0)` for simpler implementations.

"""
newton(f, fp, x0; kwargs...) = find_zero((f, fp), x0, Newton(); kwargs...)

## Halley
abstract type AbstractHalleyLikeMethod <: AbstractUnivariateZeroMethod end
fn_argout(::AbstractHalleyLikeMethod) = 3

"""
    Roots.Halley()

Implements Halley's [method](http://tinyurl.com/yd83eytb),
`xᵢ₊₁ = xᵢ - (f/f')(xᵢ) * (1 - (f/f')(xᵢ) * (f''/f')(xᵢ) * 1/2)^(-1)`
This method is cubically converging, but requires more function calls per step (3) than
other methods.

Example
```
find_zero((sin, cos, x->-sin(x)), 3.0, Roots.Halley())
```

If function evaluations are expensive one can pass in a function which
returns `(f, f/f',f'/f'')` as follows

```
find_zero(x -> (sin(x), sin(x)/cos(x), -cos(x)/sin(x)), 3.0, Roots.Halley())
```

This can be advantageous if the derivatives are easily computed from
the computation for f, but otherwise would be expensive to compute separately.

The error, `eᵢ = xᵢ - α`, satisfies
`eᵢ₊₁ ≈ -(2f'⋅f''' -3⋅(f'')²)/(12⋅(f'')²) ⋅ eᵢ³` (all evaluated at `α`).

"""
struct Halley <: AbstractHalleyLikeMethod
end

struct HalleyState{T,S} <: AbstractUnivariateZeroState{T,S}
    xn1::T
    xn0::T
    Δ::T
    ΔΔ::T
    fxn1::S
    fxn0::S
end

# we compute one step here to get x₁
function init_state(M::AbstractHalleyLikeMethod, F::Callable_Function, x)
    x₀ = float(first(x))
    fx₀, (Δ, ΔΔ) = F(x₀)
    x₁ = x₀ - 2ΔΔ/(2ΔΔ - Δ) * Δ
    state = init_state(M, F, x₀, x₁, fx₀, fx₀)
end

function init_state(::AbstractHalleyLikeMethod, F, x₀, x₁, fx₀, fx₁)
    fx₁, (Δ, ΔΔ) = F(x₁)
    HalleyState(x₁, x₀, Δ, ΔΔ, fx₁, fx₀)
end

initial_fncalls(M::AbstractHalleyLikeMethod) = 2*3

function update_state(method::Halley, F, o::HalleyState{T,S}, options::UnivariateZeroOptions, l=NullTracks()) where {T,S}
    xn = o.xn1
    fxn = o.fxn1
    r1, r2 = o.Δ, o.ΔΔ

    Δ =  2*r2/(2r2 - r1) * r1
    if isissue(Δ)
        log_message(l, "Issue with computing `Δ`")
        return (o, true)
    end

    xn1::T = xn - Δ
    fxn1::S, (r1::T, r2::T) = F(xn1)
    incfn(l,3)

    @set! o.xn0 = xn
    @set! o.xn1 = xn1
    @set! o.fxn0 = fxn
    @set! o.fxn1 = fxn1
    @set! o.Δ = r1
    @set! o.ΔΔ = r2

    return o, false

end

"""
    Roots.halley(f, fp, fpp, x0; kwargs...)

Implementation of Halley's method (cf `?Roots.Halley()`).

Arguments:

* `f::Function` -- function to find zero of

* `fp::Function` -- derivative of `f`.

* `fpp:Function` -- second derivative of `f`.

* `x0::Number` -- initial guess

With the `FowardDiff` package derivatives may be computed automatically. For example,  defining
`D(f) = x -> ForwardDiff.derivative(f, float(x))` allows `D(f)` and `D(D(f))` to be used for the first and second
derivatives, respectively.

Keyword arguments are passed to `find_zero` using the `Roots.Halley()` method.

"""
halley(f, fp, fpp, x0; kwargs...) = find_zero((f, fp, fpp), x0, Halley(); kwargs...)


text:
using InteractionNetworkModels, Distributions, BenchmarkTools, Distances, StructuredDistances
using Plots

# The Model(s)
model_mode = Hollywood(-3.0, Poisson(7), 10)
S = sample(model_mode, 10)
S = [[1,1,1,1], [2,2,2,2], [3,3,3,3]]
V = 1:20
d = FastEditDistance(FastLCS(51),51)
model = SIS(S, 4.9, d, V, 50, 50)
# model = SIS(S, 2.9, d, V)

mcmc_sampler = SisMcmcInsertDelete(
    ν_ed=2, β=0.6, ν_td=2,
    lag=1,
    K=200)

mcmc_sampler.prop_pointers

@time out = mcmc_sampler(
    model,
    lag=20,
    init=model.mode,
    desired_samples=2000,
    burn_in=0
)

plot(out)
summaryplot(out)
S
out.sample


text:
import JuLIP.Potentials: @pot, @D, evaluate!, evaluate_d!, PairPotential, @analytic
using LinearAlgebra: norm
using BenchmarkTools
using JuLIP, Test

h3("generate hand-coded morse potential")

mutable struct Morseold <: PairPotential
   e0::Float64
   A::Float64
   r0::Float64
end

@pot Morseold

evaluate!(tmp, p::Morseold, r::Number) =
      p.e0 * (exp(-2*p.A*(r/p.r0-1.0)) - 2.0*exp(-p.A*(r/p.r0-1.0)))
evaluate_d!(tmp, p::Morseold, r::Number) =
      -2.0*p.e0*p.A/p.r0*(exp(-2*p.A*(r/p.r0-1.0)) - exp(-p.A*(r/p.r0-1.0)))

A = 4.1
e0 = 0.99
r0 = 1.05

morseold = Morseold(e0, A, r0)

h3("generate AD morse potential")

morse1 = let A = A, e0 = e0, r0 = r0
   @analytic( r -> e0 * (exp(-2*A*(r/r0-1.0)) - 2.0*exp(-A*(r/r0-1.0))) )
end
@show typeof(morse1)

h3("Check consistency of hand-coded and analytic Morse potentials...")
rr = collect(range(0.9, stop=2.1, length=100))
morse_r = morse1.(rr)
morseold_r = morseold.(rr)
dmorse_r = [(@D morse1(r)) for r in rr]
dmorseold_r = [ (@D morseold(r)) for r in rr ]

println(@test norm(morse_r - morseold_r, Inf) < 1e-12)
println(@test norm(dmorse_r - dmorseold_r, Inf) < 1e-12)

function runn(f, x, N)
   s = 0.0
   for n = 1:N
      x += 0.0001
      s += f(x)
   end
   return s
end

function runn_d(f, x, N)
   s = 0.0
   for n = 1:N
      x += 0.0001
      s += @D f(x)
   end
   return s
end

h2("Performance tests: @analytic vs hand-coded")
x = 1.0+rand()
println("Evaluations of @analytic Potential")
@btime runn($morse1, $x, 1_000)
println("Evaluations hand-coded Potential")
@btime runn($morseold, $x, 1_000)
println("Grad of @analytic Potential")
@btime runn_d($morse1, $x, 1_000)
println("Grad of hand-coded Potential")
@btime runn_d($morseold, $x, 1_000)


text:
# VERTICES TEST
# Plotting scenarios and its convex hull
p= scatter(DemandData15,WindData15, xlabel="Demand [p.u.]", ylabel = "Wind [p.u.]", label = "", markersize = 16)

    Points = [DemandData15 WindData15]

    points = N -> [Points[i,:] for i in 1:N]
    v = points(15)
    hull = convex_hull(v)
    K = length(hull)

    plot!(p, VPolygon(hull), alpha=0.2, line =4)
    Hull = hull2array(hull)
    scatter!(p, Hull[:,1], Hull[:,2], color = :red, markersize= 16 , label = "")
    display(p)
    savefig(figsDir*"/15Hull.pdf")

# Plotting RES scenarios and its convex hull
Points = WindData[1:2,1:50]
fig = plot(Points, xlabel="Hour", ylabel = "RES [p.u.]", label = "", line = 4)

v = vrep(Points')

    # Constructs a polyhedon from this V-representation with the QHull library
p = polyhedron(v, QHull.Library())
    # Removing redundant points, i.e. points which are in the interior of the convex hull
removevredundancy!(p)
    # Show remaining points, i.e. the non-redundant ones
@show vrep(p)
    # Show the H-representation, the facets describing the polytope
@show hrep(p)
hull = vrep(p)
    # plot!(p, VPolygon(hull), alpha=0.2, line =4)
    # Hull = hull2array(hull)
plot(p, color = :red , label = "Convex Hull")
    # savefig(figsDir*"/15Hull.pdf")

max_stars_repo_name: gonzalez-alvaro/MultidimensionalGAPM

max_issues_repo_name: gonzalez-alvaro/MultidimensionalGAPM

max_forks_repo_name: gonzalez-alvaro/MultidimensionalGAPM


text:
using AutoHashEquals

@auto_hash_equals struct LPValue
    m::Number
    n::Number
end

LPValue(n) = LPValue(0,n)
Base.one(::Type{LPValue}) = LPValue(1)
Base.zero(::Type{LPValue}) = LPValue(0)
Base.zero(::LPValue) = LPValue(0)
Base.convert(::Type{Number},x::LPValue) = x.m*100000+x.n
Base.transpose(x::String) = x

function Base.isless(x::LPValue,y::LPValue)
    if x.m != y.m
        return x.m < y.m
    else
        return x.n < y.n
    end
    return true
end

function Base.isless(x::LPValue,y::Number)
    return isless(x,LPValue(y))
end

function Base.min(x::LPValue,y::LPValue)
    if x < y return x end
    if y < x return y end
    if x == y return x end
end

function Base.abs(x::LPValue)
    return LPValue(abs(x.m),abs(x.n))
end

function Base.inv(x::LPValue)
    @assert x.m == 0
    return LPValue(0,inv(x.n))
end

function Base.adjoint(x::LPValue)
    @assert x.m == 0
    return LPValue(0,inv(x.n))
end

function Base.:(*)(x::LPValue,y::LPValue)
    @assert y.m == 0
    return LPValue(x.m*y.n,x.n*y.n)
end

function Base.:(*)(x::LPValue,y::Number)
    return LPValue(x.m*y,x.n*y)
end

function Base.:(/)(x::LPValue,y::LPValue)
    @assert x.m == 0
    return LPValue(x.m/y.n,x.n/y.n)
end

function Base.:(+)(x::LPValue,y::LPValue)
    return LPValue(x.m+y.m,x.n+y.n)
end

function Base.:(+)(x::LPValue,y::Number)
    return LPValue(x.m,x.n+y)
end

function Base.:(+)(x::Number,y::LPValue)
    return LPValue(y.m,y.n+x)
end

function Base.:(-)(x::LPValue,y::LPValue)
    return LPValue(x.m-y.m,x.n-y.n)
end

function Base.:(-)(x::LPValue,y::Number)
    return LPValue(x.m,x.n-y)
end

function Base.:(-)(x::Number,y::LPValue)
    return LPValue(-y.m,y.n-x)
end

function Base.show(io::IO, x::LPValue)
    hasm = false
    if x.m != 0
        print(io,string(x.m))
        print(io,"M")
        hasm = true
    end
    if x.n != 0
        if x.n > 0 && hasm
            print(io,"+")
        end
        print(io,string(x.n))
    end
    if x.n == 0 && x.m == 0
        print(io,"0")
    end
end

max_stars_repo_name: mkualquiera/egypt

max_issues_repo_name: mkualquiera/egypt

max_forks_repo_name: mkualquiera/egypt


text:
export run_acdcpf

""
function run_acdcpf(file::String, model_type::Type, solver; kwargs...)
    data = _PM.parse_file(file)
    PowerModelsACDC.process_additional_data!(data)
    return run_acdcpf(data::Dict{String,Any}, model_type, solver; ref_extensions = [add_ref_dcgrid!], kwargs...)
end

""
function run_acdcpf(data::Dict{String,Any}, model_type::Type, solver; kwargs...)
    return _PM.run_model(data, model_type, solver, post_acdcpf; ref_extensions = [add_ref_dcgrid!], kwargs...)
end

""
function post_acdcpf(pm::_PM.AbstractPowerModel)
    _PM.variable_bus_voltage(pm, bounded = false)
    _PM.variable_gen_power(pm, bounded = false)
    _PM.variable_branch_power(pm, bounded = false)

    # dirty, should be improved in the future TODO
    if typeof(pm) <: _PM.SOCBFPowerModel
        _PM.variable_branch_current(pm, bounded = false)
    end

    variable_active_dcbranch_flow(pm, bounded = false)
    variable_dcbranch_current(pm, bounded = false)
    variable_dc_converter(pm, bounded = false)
    variable_dcgrid_voltage_magnitude(pm, bounded = false)

    _PM.constraint_model_voltage(pm)
    constraint_voltage_dc(pm)

    for (i,bus) in _PM.ref(pm, :ref_buses)
        @assert bus["bus_type"] == 3
        _PM.constraint_theta_ref(pm, i)
        _PM.constraint_voltage_magnitude_setpoint(pm, i)
    end

    for (i, bus) in _PM.ref(pm, :bus)# _PM.ids(pm, :bus)
        constraint_power_balance_ac(pm, i)
        # PV Bus Constraints
        if length(_PM.ref(pm, :bus_gens, i)) > 0 && !(i in _PM.ids(pm,:ref_buses))
            # this assumes inactive generators are filtered out of bus_gens
            @assert bus["bus_type"] == 2
            _PM.constraint_voltage_magnitude_setpoint(pm, i)
            for j in _PM.ref(pm, :bus_gens, i)
                _PM.constraint_gen_setpoint_active(pm, j)
            end
        end
    end

    for i in _PM.ids(pm, :branch)
        # dirty, should be improved in the future TODO
        if typeof(pm) <: _PM.SOCBFPowerModel
            _PM.constraint_power_losses(pm, i)
            _PM.constraint_voltage_magnitude_difference(pm, i)
            _PM.constraint_branch_current(pm, i)
        else
            _PM.constraint_ohms_yt_from(pm, i)
            _PM.constraint_ohms_yt_to(pm, i)
        end
    end
    for i in _PM.ids(pm, :busdc)
        constraint_power_balance_dc(pm, i)
    end
    for i in _PM.ids(pm, :branchdc)
        constraint_ohms_dc_branch(pm, i)
    end
    for (c, conv) in _PM.ref(pm, :convdc)
        constraint_conv_transformer(pm, c)
        constraint_conv_reactor(pm, c)
        constraint_conv_filter(pm, c)
        if conv["type_dc"] == 2
            constraint_dc_voltage_magnitude_setpoint(pm, c)
            constraint_reactive_conv_setpoint(pm, c)
        else
            if conv["type_ac"] == 2
                constraint_active_conv_setpoint(pm, c)
            else
                constraint_active_conv_setpoint(pm, c)
                constraint_reactive_conv_setpoint(pm, c)
            end
        end
        constraint_converter_losses(pm, c)
        constraint_converter_current(pm, c)
    end
end


text:
@testset "1035.uncrossed-lines.jl" begin
    @test max_uncrossed_lines([1, 4, 2], [1, 2, 4]) == 2
    @test max_uncrossed_lines([2, 5, 1, 2, 5], [10, 5, 2, 1, 5, 2]) == 3
    @test max_uncrossed_lines([1, 3, 7, 1, 7, 5], [1, 9, 2, 5, 1]) == 2
end


text:
# load requried packages
using IncrementalInference
using Test

## parameters

lm_prior_noise = 0.01
meas_noise = 0.25
odom_noise = 0.1
n_samples = 200

# initialize mean landmark locations
l0 = 0.0
l1 = 10.0
l2 = 40.0

# "Ground-truth" robot poses
x0 = 0.0
x1 = 10.0
x2 = 20.0
x3 = 40.0

@testset "2door basic binary multihypothesis tests..." begin

## Initialize empty factor graph
fg = initfg()
fg.solverParams.N = n_samples
# Place strong prior on locations of three "doors"
addVariable!(fg, :l0, ContinuousScalar, N=n_samples)
addFactor!(fg, [:l0], Prior(Normal(l0, lm_prior_noise)))

addVariable!(fg, :l1, ContinuousScalar, N=n_samples)
addFactor!(fg, [:l1], Prior(Normal(l1, lm_prior_noise)))

# Add first pose
addVariable!(fg, :x0, ContinuousScalar, N=n_samples)

# Make first "door" measurement
# addFactor!(fg, [:x0; :l0], LinearRelative(Normal(0, meas_noise)))
addFactor!(fg, [:x0; :l0; :l1], LinearRelative(Normal(0, meas_noise)), multihypo=[1.0; 1.0/2.0; 1.0/2.0])

# Add second pose
addVariable!(fg, :x1, ContinuousScalar, N=n_samples)

# Gaussian transition model
addFactor!(fg, [:x0; :x1], LinearRelative(Normal(x1-x0, odom_noise)))

# Make second "door" measurement
# addFactor!(fg, [:x1; :l1], LinearRelative(Normal(0, meas_noise)) )
addFactor!(fg, [:x1; :l0; :l1], LinearRelative(Normal(0, meas_noise)), multihypo=[1.0; 1.0/2.0; 1.0/2.0])

## Add one more pose/odometry to invoke issue #236

# Add third pose
addVariable!(fg, :x2, ContinuousScalar, N=n_samples)
addFactor!(fg, [:x1; :x2], LinearRelative(Normal(x2-x1, odom_noise)))

## Solve graph

tree, smt, hist = solveTree!(fg)
# drawGraph(fg)
# drawTree(tree, show=true)

@test abs(getKDEMean(getKDE(fg, :x0))[1]-x0) < 2.0
@test abs(getKDEMean(getKDE(fg, :x1))[1]-x1) < 2.0
@test abs(getKDEMean(getKDE(fg, :x2))[1]-x2) < 2.0

@test abs(getKDEMean(getKDE(fg, :l0))[1]-l0) < 3.0
@test abs(getKDEMean(getKDE(fg, :l1))[1]-l1) < 3.0

0
end

#
# using RoMEPlotting
# Gadfly.set_default_plot_size(35cm,25cm)
#
# plotKDE(fg, [:l0;:l1])
# plotKDE(fg, [:x0;:x1;:x2])

@testset "3door basic trinary multihypo test..." begin

# Initialize empty factor graph
fg = initfg()

# Place strong prior on locations of three "doors"
addVariable!(fg, :l0, ContinuousScalar, N=n_samples)
addFactor!(fg, [:l0], IIF.Prior(Normal(l0, lm_prior_noise)))

addVariable!(fg, :l1, ContinuousScalar, N=n_samples)
addFactor!(fg, [:l1], IIF.Prior(Normal(l1, lm_prior_noise)))

addVariable!(fg, :l2, ContinuousScalar, N=n_samples)
addFactor!(fg, [:l2], IIF.Prior(Normal(l2, lm_prior_noise)))

# Add first pose
addVariable!(fg, :x0, ContinuousScalar, N=n_samples)

# Make first "door" measurement
addFactor!(fg, [:x0; :l0; :l1; :l2], LinearRelative(Normal(0, meas_noise)), multihypo=[1.0; 1.0/3.0; 1.0/3.0; 1.0/3.0])

# Add second pose
addVariable!(fg, :x1, ContinuousScalar, N=n_samples)

# # Gaussian transition model
addFactor!(fg, [:x0; :x1], LinearRelative(Normal(x1-x0, odom_noise)))

# # Make second "door" measurement
addFactor!(fg, [:x1; :l0; :l1; :l2], LinearRelative(Normal(0, meas_noise)), multihypo=[1.0; 1.0/3.0; 1.0/3.0; 1.0/3.0])

# # Add third pose
addVariable!(fg, :x2, ContinuousScalar, N=n_samples)

# # Gaussian transition model
addFactor!(fg, [:x1; :x2], LinearRelative(Normal(x2-x1, odom_noise)))

# Add fourth pose
addVariable!(fg, :x3, ContinuousScalar, N=n_samples)

# Add odometry transition and new landmark sighting
addFactor!(fg, [:x2, :x3], LinearRelative(Normal(x3-x2, odom_noise)))
addFactor!(fg, [:x3; :l0; :l1; :l2], LinearRelative(Normal(0, meas_noise)), multihypo=[1.0; 1.0/3.0; 1.0/3.0; 1.0/3.0])

ensureAllInitialized!(fg)

# drawGraph(fg)
solveTree!(fg)

@test abs(getKDEMean(getKDE(fg, :x0))[1]-x0) < 2.0
@test abs(getKDEMean(getKDE(fg, :x1))[1]-x1) < 2.0
@test abs(getKDEMean(getKDE(fg, :x2))[1]-x2) < 2.0

@test abs(getKDEMean(getKDE(fg, :l0))[1]-l0) < 3.0
@test abs(getKDEMean(getKDE(fg, :l1))[1]-l1) < 3.0
@test abs(getKDEMean(getKDE(fg, :l2))[1]-l2) < 3.0

end

# using RoMEPlotting
# Gadfly.set_default_plot_size(35cm,25cm)
#
# plotKDE(fg, [:l0;:l1;:l2])
# plotKDE(fg, [:x0;:x1;:x2])

#


text:
module BenchTuple

using Benchmarkplotting
using Statistics
using Gadfly
using MLStyle
using DataFrames
import Match

import Base.getindex
getindex(asoc_lst :: Vector{Pair{Symbol, T}}, key ::Symbol) where T =
    for (search_key, value) in asoc_lst
        if search_key === key
            return value
        end
    end

"""
3 cases could succeed in validation:

(1, (_, "2", _), ("3", 4, 5))
(_, "1", 2, _, (3, "4", _), _)
(_, 1, _, 2, _, 3, _, 4, _, 5)
((1, 2, 3, _), (4, 5, 6, _, (7, 8, 9, _, (11, 12, 13))))
"""

data = [
    Symbol("ok_case1") => (1, (42, "2", 42), ("3", 4, 5)),
    Symbol("ok_case2") => (42, "1", 2, 42, (3, "4", 42), 42),
    Symbol("ok_case3") => (42, 1, 42, 2, 42, 3, 42, 4, 42, 5),
    Symbol("ok_case4") => ((1, 2, 3, 42), (4, 5, 6, 42, (7, 8, 9, 42, (11, 12, 13)))),
    Symbol("fail1") => (1, (42, 42, 42), ("3", 4, 5)),
    Symbol("fail2") => (42, "1", 2, 42, (3, 42, 42), 42),
    Symbol("fail3") => (42, 1, 42, 2, 42, 3, 42, 4, 42, 42)
]

implementations = [
    :MLStyle => (@λ begin
        (1, (_, "2", _), ("3", 4, 5))   -> 1
        (_, "1", 2, _, (3, "4", _), _)  -> 2
        (_, 1, _, 2, _, 3, _, 4, _, 5)  -> 3
        ((1, 2, 3, _), (4, 5, 6, _, (7, 8, 9, _, (11, 12, 13)))) -> 4
        _ -> 5
    end),
    Symbol("Match.jl") => function (x)
        Match.@match x begin
            (1, (_, "2", _), ("3", 4, 5))   => 1
            (_, "1", 2, _, (3, "4", _), _)  => 2
            (_, 1, _, 2, _, 3, _, 4, _, 5)  => 3
            ((1, 2, 3, _), (4, 5, 6, _, (7, 8, 9, _, (11, 12, 13)))) => 4
            _ => 5
        end
    end,
    :HandWritten => function(tp)
        !(tp isa Tuple)  ? 5 :
        let n = length(tp)
            if n === 3 && tp[1] === 1
                tp_ = tp[2]
                if !(tp isa Tuple) || tp_[2] != "2"
                    return 5
                end
                tp[3] == ("3", 4, 5) ? 1 : 5
            elseif n === 6 && tp[2] == "1" && tp[3] === 2
                tp = tp[5]
                !(tp isa Tuple) ? 5 :
                tp[1] === 3 && tp[2] == "4" ? 2 : 5
            elseif n ===  10 &&
                   tp[2] === 1 &&
                   tp[4] === 2 &&
                   tp[6] === 3 &&
                   tp[8] === 4 &&
                   tp[10] === 5
                3
            elseif n === 2 && tp[1] isa Tuple && tp[2] isa Tuple
                @inline eqtp(a, slice, v) =
                    all(slice) do i
                        a[i] === v[i]
                    end
                (a, b) = tp
                if eqtp(a, 1:3, (1, 2, 3)) &&
                   eqtp(b, 1:3, (4, 5, 6)) &&
                    let a = b[5]
                        eqtp(a, 1:3, (7, 8, 9)) &&
                        a[5] === (11, 12, 13)
                    end
                    4
                end
            else
                5
            end
        end
    end
]

criterion(x) = (meantime = mean(x.times), allocs = 1 + x.allocs)
df = bcompare(criterion, data, implementations)
@info df
theme = Theme(
    guide_title_position = :left,
    colorkey_swatch_shape = :circle,
    minor_label_font = "Consolas",
    major_label_font = "Consolas",
    point_size=5px
)
report_meantime, df_time = report(:meantime, df, Scale.y_log10, theme)
report_allocs, df_allocs = report(:allocs, df, theme)

open("stats/vs-hw(tuple).txt", "w") do f
    write(f, string(df))
end

draw(SVG("stats/vs-hw(tuple)-on-time.svg", 10inch, 4inch), report_meantime);
draw(SVG("stats/vs-hw(tuple)-on-allocs.svg", 10inch, 4inch), report_allocs);

end


text:
using DataFrames
using CSV
using Plots
#%%
data_v = open("pendulum.dat")
df_v=CSV.read(data_v, header=false, delim=" ")
t_v=df_v[1]
theta_v=df_v[2]
omega_v=df_v[3]

data_rk = open("pendulum_rk.dat")
df_rk=CSV.read(data_rk, header=false, delim=" ")
t_rk=df_rk[1]
theta_rk=df_rk[2]
omega_rk=df_rk[3]
#%%
plot(df_v[1], df_v[2], label="theta", dpi=600)
plot(df_rk[1], df_rk[2], label="theta", dpi=600)
# xaxis!(:log10)
# yaxis!(:log10)
#%%
E_v = (omega_v.^2)/2 + (ones(length(omega_v))-cos.(theta_v))
E_rk = (omega_rk.^2)/2 + (ones(length(omega_rk))-cos.(theta_rk))

plot(t_v, E_v, label="Verlet", dpi=600)
plot!(t_rk, E_rk, label="RK2", dpi=600)
plot!(legend=:topleft)
# ylims!(0,0.4)
ylabel!("E/(m L^2)")
xlabel!("step")
savefig("pendulum.png")

max_stars_repo_name: aurelio-amerio/corso-cpp


text:
# This file is a part of ValueShapes.jl, licensed under the MIT License (MIT).

"""
    ConstValueDist <: Distributions.Distribution

Represents a delta distribution for a constant value of arbritrary type.

Calling `varshape` on a `ConstValueDist` will yield a
[`ConstValueShape`](@ref).
"""
struct ConstValueDist{VF<:VariateForm,T} <: Distribution{VF,Continuous}
    value::T
end

export ConstValueDist

ConstValueDist(x::T) where {T<:Real} = ConstValueDist{Univariate,T}(x)
ConstValueDist(x::T) where {T<:AbstractVector{<:Real}} = ConstValueDist{Multivariate,T}(x)
ConstValueDist(x::T) where {T<:AbstractMatrix{<:Real}} = ConstValueDist{Matrixvariate,T}(x)

Distributions.pdf(d::ConstValueDist{Univariate}, x::Real) = d.value == x ? float(eltype(d))(Inf) : float(eltype(d))(0)
Distributions._logpdf(d::ConstValueDist, x::AbstractArray{<:Real}) = d.value == x ? float(eltype(d))(Inf) : float(eltype(d))(0)

Distributions.cdf(d::ConstValueDist{Univariate}, x::Real) = d.value <= x ? Float32(1) : Float32(0)
Distributions.quantile(d::ConstValueDist{Univariate}, q::Real) = d.value # Sensible?
Distributions.minimum(d::ConstValueDist{Univariate}) = d.value
Distributions.maximum(d::ConstValueDist{Univariate}) = d.value
Distributions.insupport(d::ConstValueDist{Univariate}, x::Real) = x == d.value

StatsBase.mode(d::ConstValueDist) = d.value

Base.size(d::ConstValueDist) = size(d.value)
Base.length(d::ConstValueDist) = prod(size(d))
Base.eltype(d::ConstValueDist) = eltype(d.value)

Random.rand(rng::AbstractRNG, d::ConstValueDist) = d.value

function Random._rand!(rng::AbstractRNG, d::ConstValueDist, x::AbstractArray{<:Real})
    copyto!(x, d.value)
end

ValueShapes.varshape(d::ConstValueDist) = ConstValueShape(d.value)


text:
nt = 1000
m = ABCDQR_simula(nt)

y = m["y"]
u = m["u"]

############################### starting point
nx = 4
i = nx+1
Ai,Ci,Qi,Ri,Si = EmACQR.ACQRS_sub(y,nx,i)
Bi = zeros(nx,2)
Di = zeros(2,2)
m1i = Ci\y[:,1:1]
P1i = 0*Ai

#######################
max_iter=10
tol=1e-6
txo=false
print("Testing ABCDQR_em: ")
A,B,C,D,Q,R,m1,P1,loglik,aic = ABCDQR_em(y,u,Ai,Bi,Ci,Di,Qi,Ri,m1i,P1i,max_iter,tol,txo)
@test sum(A.^2) > 1e-6
println("OK")

#######################
max_iter=10
tol=1e-10
txo=false
print("Testing ABCDQR_em_s: ")
A,B,C,D,Q,R,m1,P1,loglik,aic = ABCDQR_em_s(y,u,Ai,Bi,Ci,Di,Qi,Ri,m1i,max_iter,tol,txo)
@test sum(A.^2) > 1e-6
println("OK")


text:
const ActionSource = Union{PkgId,Module}
const sym2loader = Dict{Symbol,Vector{ActionSource}}()
const sym2saver  = Dict{Symbol,Vector{ActionSource}}()

for (applicable_, add_, dict_) in (
        (:applicable_loaders, :add_loader, :sym2loader),
        (:applicable_savers,  :add_saver,  :sym2saver))
    @eval begin
        function $applicable_(sym::Symbol)
            ret = get($dict_, sym, nothing)
            ret === nothing && error(string("No ", $applicable_, " found for ", sym))
            return ret
        end
        $add_(@nospecialize(fmt::Type), id::Union{ActionSource,Pair}) = $add_(formatname(fmt)::Symbol, id)
        function $add_(sym::Symbol, id::ActionSource)
            list = get!(Vector{ActionSource}, $dict_, sym)
            push!(list, id)
        end
        $add_(sym::Symbol, pkg::Pair{<:Union{String,Symbol}, UUID}) = $add_(sym, Base.PkgId(pkg.second, String(pkg.first)))
    end
end

"""
    add_loader(fmt, :Package=>uuid)
    add_loader(fmt, [:Package=>uuid, specifiers...])

Declare that format `fmt` can be loaded with package `:Package`.
Specifiers include `OSX`, `Unix`, `Windows` and `Linux` to restrict usage to particular operating systems.

See also [`add_format`](@ref) which can combine package support with the format declaration.
"""
add_loader

"""
    add_saver(fmt, :Package=>uuid)
    add_saver(fmt, [:Package=>uuid, specifiers...])

Declare that format `fmt` can be saved with package `:Package`.
Specifiers include `OSX`, `Unix`, `Windows` and `Linux` to restrict usage to particular operating systems.

See also [`add_format`](@ref) which can combine package support with the format declaration.
"""
add_saver

"""
- `load(filename)` loads the contents of a formatted file, trying to infer
  the format from `filename` and/or magic bytes in the file (see [`query`](@ref)).
- `load(strm)` loads from an `IOStream` or similar object. In this case,
  there is no filename extension, so we rely on the magic bytes for format
  identification.
- `load(File{format"PNG"}(filename))` specifies the format directly, and bypasses the format [`query`](@ref).
- `load(Stream{format"PNG"}(io))` specifies the format directly, and bypasses the format [`query`](@ref).
- `load(f; options...)` passes keyword arguments on to the loader.
"""
load

"""
Some packages may implement a streaming API, where the contents of the file can
be read in chunks and processed, rather than all at once. Reading from these
higher-level streams should return a formatted object, like an image or chunk of
video or audio.

- `loadstreaming(filename)` loads the contents of a formatted file, trying to infer
  the format from `filename` and/or magic bytes in the file. It returns a streaming
  type that can be read from in chunks, rather than loading the whole contents all
  at once.
- `loadstreaming(strm)` loads the stream from an `IOStream` or similar object.
  In this case, there is no filename extension, so we rely on the magic bytes
  for format identification.
- `loadstreaming(File{format"WAV"}(filename))` specifies the format directly, and
  bypasses the format [`query`](@ref).
- `loadstreaming(Stream{format"WAV"}(io))` specifies the format directly, and
  bypasses the format [`query`](@ref).
- `loadstreaming(f; options...)` passes keyword arguments on to the loader.
"""
loadstreaming

"""
- `save(filename, data...)` saves the contents of a formatted file,
  trying to infer the format from `filename`.
- `save(Stream{format"PNG"}(io), data...)` specifies the format directly, and
  bypasses the format [`query`](@ref).
- `save(File{format"PNG"}(filename), data...)` specifies the format directly, and
  bypasses the format [`query`](@ref).
- `save(f, data...; options...)` passes keyword arguments on to the saver.
"""
save

"""
Some packages may implement a streaming API, where the contents of the file can
be written in chunks, rather than all at once. These higher-level streams should
accept formatted objects, like an image or chunk of video or audio.

- `savestreaming(filename, data...)` saves the contents of a formatted file,
  trying to infer the format from `filename`.
- `savestreaming(File{format"WAV"}(filename))` specifies the format directly, and
  bypasses the format [`query`](@ref).
- `savestreaming(Stream{format"WAV"}(io))` specifies the format directly, and
  bypasses the format [`query`](@ref).
- `savestreaming(f, data...; options...)` passes keyword arguments on to the saver.
"""
savestreaming

# if a bare filename or IO stream are given, query for the format and dispatch
# to the formatted handlers below
for fn in (:load, :loadstreaming, :metadata)
    fnq = QuoteNode(fn)
    @eval function $fn(file, args...; options...)
        checkpath_load(file)
        sym = querysym(file)
        libraries = applicable_loaders(sym)
        return action($fnq, libraries, sym, file, args...; options...)
    end
    # Version that bypasses format-inference
    @eval function $fn(@nospecialize(file::Formatted), args...; options...)
        checkpath_load(filename(file))
        sym = formatname(file)::Symbol
        libraries = applicable_loaders(sym)
        return action($fnq, libraries, file, args...; options...)
    end
end
for fn in (:save, :savestreaming)
    fnq = QuoteNode(fn)
    @eval function $fn(file, args...; options...)
        checkpath_save(file)
        sym = querysym(file; checkfile=false)
        libraries = applicable_savers(sym)
        return action($fnq, libraries, sym, file, args...; options...)
    end
    @eval function $fn(@nospecialize(file::Formatted), args...; options...)
        checkpath_save(filename(file))
        sym = formatname(file)::Symbol
        libraries = applicable_savers(sym)
        return action($fnq, libraries, file, args...; options...)
    end
    @eval function $fn(@nospecialize(fmt::Type), file, args...; options...)
        checkpath_save(file)
        sym = formatname(fmt)::Symbol
        libraries = applicable_savers(sym)
        return action($fnq, libraries, sym, file, args...; options...)
    end
end

# return a save function, so you can do `thing_to_save |> save("filename.ext")`
function save(file; options...)
    sym = querysym(file; checkfile=false)
    libraries = applicable_savers(sym)
    return data -> action(:save, libraries, sym, file, data; options...)
end

# do-syntax for streaming IO
for fn in (:loadstreaming, :savestreaming)
    @eval function $fn(@nospecialize(f::Function), @nospecialize(args...); @nospecialize(kwargs...))
        str = $fn(args...; kwargs...)
        try
            f(str)
        finally
            close(str)
        end
    end
end

function checkpath_load(file)
    file === nothing && return nothing   # likely stream io
    isa(file, IO) && return nothing
    _isurl(file) || isfile(file) || throw(ArgumentError("No file exists at given path: $file"))
    return nothing
end
function checkpath_save(file)
    file === nothing && return nothing
    isa(file, IO) && return nothing
    isdir(file) && throw(ArgumentError("Given file path is a directory: $file"))
    dn = dirname(file)
    !isdir(dn) && mkpath(dn)
    return nothing
end

# TODO: we may use URI when cases become complex
_isurl(path::AbstractString) = startswith(path, "http://") || startswith(path, "https://")
_isurl(path) = false

action(call::Symbol, libraries::Vector{ActionSource}, sym::Symbol, io::IO, args...; options...) =
    action(call, libraries, Stream{DataFormat{sym}}(io), args...; options...)
action(call::Symbol, libraries::Vector{ActionSource}, sym::Symbol, file, args...; options...) =
    action(call, libraries, File{DataFormat{sym}}(file), args...; options...)

# To test for broken packages which extend FileIO functions
const fileiofuncs = Dict{Symbol,Function}(:load => load,
                                          :loadstreaming => loadstreaming,
                                          :metadata => metadata,
                                          :save => save,
                                          :savestreaming => savestreaming)

const require_lock = ReentrantLock()
function action(call::Symbol, libraries::Vector{ActionSource}, @nospecialize(file::Formatted), args...; options...)
    issave = call ∈ (:save, :savestreaming)
    failures = Tuple{Any,ActionSource,Vector}[]
    pkgfuncname = Symbol("fileio_", call)
    local mod
    for library in libraries
        try
            mod = isa(library, Module) ? library : lock(() -> Base.require(library), require_lock)
            f = if isdefined(mod, pkgfuncname)
                getfield(mod, pkgfuncname)
            elseif isdefined(mod, call)
                getfield(mod, call)
            else
                throw(SpecError(mod, call))
            end
            if f === fileiofuncs[call]
                argtyps = map(Core.Typeof, args)
                m = which(f, (typeof(file), argtyps...))
                if m == which(f, (Formatted, argtyps...))  # is it equal to a fallback?
                    throw(SpecError(mod, call))            # prevent a StackOverflowError
                end
                @warn "$mod incorrectly extends FileIO functions (see FileIO documentation)"
            end
            return Base.invokelatest(f, file, args...; options...)
        catch e
            if isa(e, SpecError)
                str = "neither $call nor $pkgfuncname is defined"
                e = issave ? WriterError(string(mod), str, e) : LoaderError(string(mod), str, e)
            end
            push!(failures, (e, library, catch_backtrace()))
        end
    end
    handle_exceptions(failures, "$call $(repr(file))")
end


text:
@testset "get_output_str" begin
    @test GraphQLClient.get_output_str("field1") == "field1,"
    @test GraphQLClient.get_output_str(Alias("alias_name", "field_name")) == "alias_name:field_name,"
    @test GraphQLClient.get_output_str(["field1"]) == "field1,"
    @test GraphQLClient.get_output_str(["field1", "field2"]) == "field1,field2,"
    @test GraphQLClient.get_output_str(Dict("field1" => "field2")) == "field1{field2,},"
    @test GraphQLClient.get_output_str(Dict("field1" => ["field2"])) == "field1{field2,},"
    @test GraphQLClient.get_output_str(Dict("field1" => ["field2", "field3"])) == "field1{field2,field3,},"
    @test GraphQLClient.get_output_str(["field1", Dict("field2" => ["field3"])]) == "field1,field2{field3,},"
end

@testset "_get_fields_str" begin
    client = Client("url", "ws"; introspect=false)

    # Not an object, list of objects or nonnull
    field = build_type(nothing, "field_name", type=build_type("SCALAR", "String"))
    @test GraphQLClient._get_field_str(client, field, [], []) == "field_name,"

    # Empty Object
    field = build_type(nothing, "object_field_name", type=build_type("OBJECT", "MyObject"))
    client.type_to_fields_map["MyObject"] = Dict()
    @test GraphQLClient._get_field_str(client, field, [], []) == ""

    # Object with two fields
    field1 = build_type(nothing, "field1", type=build_type("SCALAR", "String"))
    field2 = build_type(nothing, "field2", type=build_type("SCALAR", "Float"))
    client.type_to_fields_map["MyObject"] = Dict(
        "field1" => field1,
        "field2" => field2,
    )
    @test GraphQLClient._get_field_str(client, field, [], []) == "object_field_name{field1,field2,}"

    # List of scalars
    field = build_type(nothing, "scalar_field_name",
        type=build_type("LIST", nothing,
            ofType=build_type("SCALAR", "String")))
    @test GraphQLClient._get_field_str(client, field, [], []) == "scalar_field_name,"

    # List of objects with two fields
    field = build_type(nothing, "object_field_name",
        type=build_type("LIST", nothing,
            ofType=build_type("OBJECT", "MyObject")))
    @test GraphQLClient._get_field_str(client, field, [], []) == "object_field_name{field1,field2,}"

    # Nonnull scalar, enum and list of enums
    field = build_type(nothing, "field_name",
        type=build_type("NON_NULL", nothing,
            ofType=build_type("SCALAR", "String")))
    @test GraphQLClient._get_field_str(client, field, [], []) == "field_name,"
    field = build_type(nothing, "enum_field_name",
        type=build_type("NON_NULL", nothing,
            ofType = build_type("ENUM", "MyEnum")))
    @test GraphQLClient._get_field_str(client, field, [], []) == "enum_field_name,"
    field = build_type(nothing, "enum_field_name",
        type=build_type("NON_NULL", nothing,
            ofType=build_type("LIST", nothing,
                ofType=build_type("ENUM", "MyEnum"))))
    @test GraphQLClient._get_field_str(client, field, [], []) == "enum_field_name,"
    field = build_type(nothing, "enum_field_name",
        type=build_type("NON_NULL", nothing,
            ofType=build_type("LIST", nothing,
                ofType=build_type("LIST", nothing,
                    ofType=build_type("ENUM", "MyEnum")))))
    @test GraphQLClient._get_field_str(client, field, [], []) == "enum_field_name,"

    # Nonnull object
    field = build_type(nothing, "object_field_name",
    type=build_type("NON_NULL", nothing,
            ofType=build_type("OBJECT", "MyObject")))
    @test GraphQLClient._get_field_str(client, field, [], []) == "object_field_name{field1,field2,}"

    # Nonnull list of objects
    field = build_type(nothing, "object_field_name",
        type=build_type("NON_NULL", nothing,
            ofType=build_type("LIST", nothing,
                ofType=build_type("OBJECT", "MyObject"))))
    @test GraphQLClient._get_field_str(client, field, [], []) == "object_field_name{field1,field2,}"

    # Recursion in objects should warn
    field1 = build_type(nothing, "field1", type=build_type("SCALAR", "String"))
    recursive_object = build_type(nothing, "object_field_name", type=build_type("OBJECT", "MyRecursiveObject"))
    client.type_to_fields_map["MyRecursiveObject"] = Dict(
        "field1" => field1,
        "field2" => recursive_object,
    )
    objects_that_recurse = String[]
    @test GraphQLClient._get_field_str(client, recursive_object, [], objects_that_recurse) == "object_field_name{field1,}"
    @test in("MyRecursiveObject", objects_that_recurse)

    # Same when object is nonnull
    recursive_object = build_type(nothing, "object_field_name",
        type=build_type("NON_NULL", nothing,
            ofType=build_type("OBJECT", "MyRecursiveObject")))
    objects_that_recurse = String[]
    @test GraphQLClient._get_field_str(client, recursive_object, [], objects_that_recurse) == "object_field_name{field1,}"
    @test in("MyRecursiveObject", objects_that_recurse)

    # Same when object is list
    recursive_object = build_type(nothing, "object_field_name",
        type=build_type("LIST", nothing,
            ofType=build_type("OBJECT", "MyRecursiveObject")))
    objects_that_recurse = String[]
    @test GraphQLClient._get_field_str(client, recursive_object, [], objects_that_recurse) == "object_field_name{field1,}"
    @test in("MyRecursiveObject", objects_that_recurse)

    # Unknown type
    field = build_type(nothing, "field_name", type=build_type("UNION", "String"))
    @test_throws GraphQLClient.GraphQLClientException GraphQLClient._get_field_str(client, field, [], [])
    # Unknown NON_NULLtype
    field = build_type(nothing, "field_name",
        type=build_type("NON_NULL", nothing,
            ofType=build_type("UNION", "String")))
    @test_throws GraphQLClient.GraphQLClientException GraphQLClient._get_field_str(client, field, [], [])
end

@testset "get_all_output_fields_str" begin
    # type_to_fields_map
    client = Client("url", "ws"; introspect=false)
    type = "MyType"
    client.type_to_fields_map[type] = Dict(
        "field1" => build_type(nothing, "field1", type=build_type("SCALAR", "Boolean")),
        "field2" => build_type(nothing, "field2", type=build_type("OBJECT", "MyObject")),
    )
    client.type_to_fields_map["MyObject"] = Dict(
        "field2" => build_type(nothing, "field2", type=build_type("SCALAR", "Float")),
    )
    @test GraphQLClient.get_field_names_string(client, type, [], []) == "field1,field2{field2,}"

    # Type with no fields
    type = "MyEmptyType"
    client.type_to_fields_map[type] = Dict()
    @test GraphQLClient.get_field_names_string(client, type, [], []) == ""

    # get_all_output_fields_str
    query_name = "myQuery"
    client.query_to_type_map[query_name] = "MyType"
    @test GraphQLClient.get_all_output_fields_str(client, query_name, []) == "field1,field2{field2,}"
    query_name = "myQuery"
    client.query_to_type_map[query_name] = "MyEmptyType"
    @test GraphQLClient.get_all_output_fields_str(client, query_name, []) == ""

    # test objects_to_ignore is working
    type = "MyRecursiveType"
    client.type_to_fields_map[type] = Dict(
        "field1" => build_type(nothing, "field1", type=build_type("SCALAR", "Boolean")),
        "field2" => build_type(nothing, "field2", type=build_type("OBJECT", "MyRecursiveType")),
    )
    client.query_to_type_map[query_name] = "MyRecursiveType"
    @test GraphQLClient.get_all_output_fields_str(client, query_name, []) == "field1,"
    @test_logs (:warn, Regex("Can't query all output fields due to recursion of these object\\(s\\):")) GraphQLClient.get_all_output_fields_str(client, query_name, [])
end

@testset "issue #7" begin
    output_fields=[
        Dict(
            "level1" =>Dict(
                "level2" => Dict(
                    "level3" => Dict(
                        "level4_1" => [
                            "level4_1_1",
                        ],
                        "level4_2" => [
                            "level4_2_1",
                        ],
                    )
                )
            )
        )
    ]
    output_str = GraphQLClient.get_output_str(output_fields)
    @test occursin("level4_1{level4_1_1,}" , output_str)
    @test occursin("level4_2{level4_2_1,}" , output_str)
end


text:
function part1()

    line1, line2 = readlines("input13.txt")
    timestamp = parse(Int, line1)
    bus_ids = parse.(Int, filter( x-> !occursin("x",x), split(line2,",")))

    flag = true
    wait = 0

    while flag
        for bus_id in bus_ids
            if (timestamp+wait) % bus_id == 0
                return wait * bus_id
            end
        end
        wait += 1
    end

end

function check( timestamp, departs )

   for (bus_id,t) in departs
      if (timestamp+t) % bus_id !== 0
         return false
      end
   end
   return true
end

# Three days to find the solution
# This page helps me https://www.geeksforgeeks.org/chinese-remainder-theorem-set-1-introduction/

function part2(departs)

    bus_ids = split(departs,",")

    departs = Tuple{Int,Int}[]
    for (minute,bus_id) in enumerate(bus_ids)
        if !occursin( "x", bus_id)
            push!(departs, (parse(Int,bus_id),minute-1))
        end
    end

    timestamp = 1
    increment = timestamp

    for (bus, offset) in departs
        while (timestamp + offset) % bus != 0
            timestamp += increment
        end
        increment = lcm(increment, bus)
    end

    return timestamp

end

println(part1())
@assert part2("17,x,13,19") == 3417
@assert part2("67,7,59,61") == 754018
@assert part2("67,x,7,59,61") == 779210
@assert part2("7,13,x,x,59,x,31,19") == 1068781
@assert part2("67,7,x,59,61") == 1261476
@time @assert part2(1, "1789,37,47,1889") == 1202161486
_, departs = readlines("input13.txt")
@time println(part2(departs))

max_stars_repo_name: pnavaro/AdventOfCode

max_issues_repo_name: pnavaro/AdventOfCode

max_forks_repo_name: pnavaro/AdventOfCode


text:
"""
    toSVector(v::AbstractVector, len::Int)

Convert a Vector `v` to a StaticVector of length `D`. If the length of `v` is larger than `len`, the vector is truncated, otherwise the vector is padded to the right with zeros.
"""
function toSVector(v::AbstractVector, len::Int)

    u = zeros(eltype(v), len)

    for i=1:min(len, length(v))
        u[i] = v[i]
    end

    return SVector{len}(u)
end

"""
    toSVector(v::AbstractVector)

Convert a Vector `v` to a StaticVector.
"""
toSVector(v::AbstractVector) = SVector{length(v)}(v)

"""
    cdot(a::StaticVector{d1, T1}, b::StaticVector{d2, T2}) where {d1, d2, T1, T2}

Take the inner product of two vectors of unequal size. The longer vector is truncated to the size of the smaller one.
"""
@generated function cdot(a::StaticVector{d1, T1}, b::StaticVector{d2, T2}) where {d1, d2, T1, T2}
    s = min(d1, d2)

    expr = :(a[1]*b[1])

    for j = 2:s
        expr = :($expr + a[$j] * b[$j])
    end

    return :($expr)
end

"""
    cdot(a::AbstractVector, b::AbsctractVector)

Take the inner product of two vectors of unequal size. The longer vector is truncated to the size of the smaller one.
"""
function cdot(a::AbstractVector, b::AbstractVector)
    s = min(length(a), length(b))

    accu = a[1]*b[1]

    for j = 2:s
        accu += accu + a[j]*b[j]
    end

    return accu
end


text:
#-------------------------------------------------------------------
#* EMSO Model Library (EML) Copyright (C) 2004 - 2007 ALSOC.
#*
#* This LIBRARY is free software; you can distribute it and/or modify
#* it under the therms of the ALSOC FREE LICENSE as available at
#* http://www.enq.ufrgs.br/alsoc.
#*
#* EMSO Copyright (C) 2004 - 2007 ALSOC, original code
#* from http://www.rps.eng.br Copyright (C) 2002-2004.
#* All rights reserved.
#*
#* EMSO is distributed under the therms of the ALSOC LICENSE as
#* available at http://www.enq.ufrgs.br/alsoc.
#*--------------------------------------------------------------------
#* Author: Gerson Balbueno Bicca
#* $Id: PHE.mso 250 2007-04-27 16:32:02Z bicca $
#*------------------------------------------------------------------
type Main_PHE
	Main_PHE()=begin
		new(
			PHE_HeatTransfer(Dict{Symbol,Any}(
				:Brief=>"PHE Heat Transfer",
				:Symbol=>" "
			)),
			PHE_PressureDrop(Dict{Symbol,Any}(
				:Brief=>"PHE Pressure Drop",
				:Symbol=>" "
			)),
			Physical_Properties(Dict{Symbol,Any}(
				:Brief=>"PHE Properties",
				:Symbol=>" "
			)),
			[:HeatTransfer,:PressureDrop,:Properties,]
		)
	end
	HeatTransfer::PHE_HeatTransfer
	PressureDrop::PHE_PressureDrop
	Properties::Physical_Properties
	variables::Array{Symbol,1}
	attributes::Dict{Symbol,Any}
end
export Main_PHE
function atributes(in::Main_PHE,_::Dict{Symbol,Any})
	fields::Dict{Symbol,Any}=Dict{Symbol,Any}()
	fields[:Pallete]=false
	fields[:Brief]="to be documented"
	fields[:Info]="to be documented"
	drive!(fields,_)
	return fields
end
Main_PHE(_::Dict{Symbol,Any})=begin
	newModel=Main_PHE()
	newModel.attributes=atributes(newModel,_)
	newModel
end


text:
function initialize_raw_loop_async(connection::HTTPConnection, buffer; skip_preface=false)
    channel_act_raw = connection.channel_act_raw
    channel_evt_raw = connection.channel_evt_raw

    ## Initialize the connection
    CLIENT_PREFACE = bytearr("PRI * HTTP/2.0\r\n\r\nSM\r\n\r\n")

    if !skip_preface
        if connection.isclient
            write(buffer, CLIENT_PREFACE)
        else
            @assert read(buffer, length(CLIENT_PREFACE)) == CLIENT_PREFACE
        end
    end

    connection.evt_raw_processor = @async begin
        try
            while true
                if connection.closed
                    break
                end

                if eof(buffer)
                    connection.closed = true
                    continue
                end

                frame = try
                    Frame.decode(buffer)
                catch
                    goaway!(connection, ProtocolError("Decode error."))
                    break
                end

                # Abstract atom headers frame away
                if typeof(frame) == HeadersFrame || typeof(frame) == PushPromiseFrame
                    continuations = Array{ContinuationFrame, 1}()
                    while !(frame.is_end_headers ||
                            (!frame.is_end_headers && length(continuations) == 0) ||
                            continuations[length(continuations)].is_end_headers)

                        continuation = try
                            Frame.decode(buffer)
                        catch
                            goaway!(connection, ProtocolError("Decode error."))
                            break
                        end

                        if !(typeof(continuation) == ContinuationFrame &&
                             continuation.stream_identifier == frame.stream_identifier)
                            goaway!(connection, ProtocolError("Headers must be followed by continuations if it is not the end."))
                        end
                        push!(continuations, continuation)
                    end

                    fragment = copy(frame.fragment)
                    if length(continuations) > 0
                        for i = 1:length(continuations)
                            append!(fragment, continuations[i].fragment)
                        end
                    end

                    if typeof(frame) == HeadersFrame
                        put!(channel_evt_raw, HeadersFrame(frame.is_end_stream,
                                                           true,
                                                           frame.is_priority,
                                                           frame.stream_identifier,
                                                           frame.exclusive,
                                                           frame.dependent_stream_identifier,
                                                           frame.weight,
                                                           fragment))
                    else
                        put!(channel_evt_raw, PushPromiseFrame(true,
                                                               frame.stream_identifier,
                                                               frame.promised_stream_identifier,
                                                               fragment))
                    end
                else
                    put!(channel_evt_raw, frame)

                    if typeof(frame) == DataFrame && !frame.is_end_stream
                        put!(channel_act_raw, WindowUpdateFrame(0, length(frame.data)))
                        put!(channel_act_raw, WindowUpdateFrame(frame.stream_identifier, length(frame.data)))
                    end

                    if typeof(frame) == GoawayFrame
                        connection.closed = true
                    end
                end
            end
        catch ex
            @info("Exception in connection.evt_raw_processor $ex")
            rethrow(ex)
        end
    end

    connection.act_raw_processor = @async begin
        try
            while true
                if connection.closed
                    break
                end

                frame = take!(channel_act_raw)

                if typeof(frame) == HeadersFrame || typeof(frame) == PushPromiseFrame
                    write(buffer, Frame.encode(frame))
                    while !frame.is_end_headers
                        continuation = take!(channel_act_raw)
                        if !(typeof(continuation) == ContinuationFrame &&
                             continuation.stream_identifier == frame.stream_identifier)
                            goaway!(connection, InternalError("Headers must be followed by continuations if it is not the end."))
                        end
                        write(buffer, Frame.encode(continuation))
                    end
                else
                    encoded = Frame.encode(frame)

                    if typeof(frame) == DataFrame
                        stream = get_stream(connection, frame.stream_identifier)
                        if stream.window_size < length(encoded)
                            #@show("putting back as window_size=$(stream.window_size), len=$(length(encoded))")
                            put!(channel_act_raw, frame)
                            continue
                        end

                        stream.window_size -= length(encoded)
                    end

                    write(buffer, encoded)

                    if typeof(frame) == GoawayFrame
                        connection.closed = true
                    end
                end
            end
        catch ex
            @info("Exception in connection.act_raw_processor $ex")
            rethrow(ex)
        end
    end

    put!(channel_act_raw, SettingsFrame(false, Array{Tuple{Frame.SETTING_IDENTIFIER, UInt32}, 1}()))
end

function process_channel_act(connection::HTTPConnection)
    channel_act_raw = connection.channel_act_raw
    channel_act = connection.channel_act
    channel_evt_raw = connection.channel_evt_raw
    channel_evt = connection.channel_evt

    act = take!(channel_act)

    if connection.last_stream_identifier <= act.stream_identifier
        connection.last_stream_identifier = act.stream_identifier
    end

    stream = get_stream(connection, act.stream_identifier)
    if stream.state == IDLE && (connection.settings.max_concurrent_streams !== nothing) &&
        concurrent_streams_count(connection) > connection.settings.max_concurrent_streams
        put!(channel_act, act)
        return
    end

    if typeof(act) == ActSendHeaders
        if (connection.settings.max_header_list_size !== nothing)
            sum = 0

            for k in keys(act.headers)
                sum += length(k) + length(act.headers[k]) + 32
            end

            if sum > connection.settings.max_header_list_size
                goaway!(connection, InternalError("Header list size exceeded."))
                return
            end
        end

        frame = send_stream_headers(connection, act)
        handle_stream_state!(connection, frame, true)
    elseif typeof(act) == ActSendData
        frame = send_stream_data(connection, act)
        handle_stream_state!(connection, frame, true)
    elseif typeof(act) == ActPromise
        if !connection.settings.push_enabled
            goaway!(connection, InternalError("Push is disabled."))
            return
        end

        frame = send_stream_push_promise(connection, act)
        handle_stream_state!(connection, frame, true)
    else
        goaway!(connection, InternalError("Unknown action for channel."))
    end
end

function process_channel_evt(connection::HTTPConnection)
    channel_act_raw = connection.channel_act_raw
    channel_act = connection.channel_act
    channel_evt_raw = connection.channel_evt_raw
    channel_evt = connection.channel_evt

    frame = take!(channel_evt_raw)

    ## Frames where stream identifier is 0x0

    if typeof(frame) == SettingsFrame
        if !frame.is_ack
            parameters = frame.parameters
            if length(parameters) > 0
                for i = 1:length(parameters)
                    handle_setting!(connection, parameters[i][1], parameters[i][2])
                end
            end
            put!(channel_act_raw, SettingsFrame(true, Array{Tuple{Frame.SETTING_IDENTIFIER, UInt32}, 1}()))
        end
        return
    end
    if typeof(frame) == PingFrame
        put!(channel_act_raw, PingFrame(true, frame.data))
        return
    end
    if typeof(frame) == GoawayFrame
        @assert frame.error_code == 0x0
        put!(channel_evt, EvtGoaway())
        return
    end
    if typeof(frame) == WindowUpdateFrame && frame.stream_identifier == 0x0
        connection.window_size += frame.window_size_increment
        return
    end

    ## Frames where stream identifier is not 0x0
    @assert frame.stream_identifier != 0x0
    handle_stream_state!(connection, frame, false)
    stream = get_stream(connection, frame.stream_identifier)

    if typeof(frame) == DataFrame
        recv_stream_data(connection, frame)
        return
    end
    if typeof(frame) == HeadersFrame
        recv_stream_headers(connection, frame)
        return
    end
    if typeof(frame) == PriorityFrame
        recv_stream_priority(connection, frame)
        return
    end
    if typeof(frame) == RstStreamFrame
        recv_stream_rst_stream(connection, frame)
        return
    end
    if typeof(frame) == PushPromiseFrame
        recv_stream_push_promise(connection, frame)
        return
    end
    if typeof(frame) == WindowUpdateFrame
        recv_stream_window_update(connection, frame)
        return
    end

    @assert false
end

function initialize_in_loop_async(connection::HTTPConnection)
    connection.act_processor = @async begin
        try
            while !connection.closed
                process_channel_act(connection)
            end
            put!(channel_evt, EvtGoaway())
        catch ex
            @info("Exception in connection.act_processor $ex")
            rethrow(ex)
        end
    end

    connection.evt_raw_processor = @async begin
        try
            while !connection.closed # || (connection.closed && isready(connection.channel_evt_raw))
                process_channel_evt(connection)
            end
        catch ex
            @info("Exception in connection.evt_raw_processor $ex")
            rethrow(ex)
        end
    end
    nothing
end

function initialize_loop_async(connection::HTTPConnection, buffer; skip_preface=false)
    initialize_raw_loop_async(connection, buffer; skip_preface=skip_preface)
    initialize_in_loop_async(connection)
end


text:
#SUMPOOL

@polly function sumpool_forward(nn, nd, ih, iw, ow, oh, dh, dw, sh, sw, inp_F, out_F)

    for n=1:nn, d=1:nd, r=1:oh, c=1:ow
        val=Float32(0.0)
        for h=10*r+1:min(10*r+dh,ih), w=10*c:min(10*c+dw,iw)
            val += inp_F[n,d,h,w]
        end
        out_F[n,d,r,c] = val
    end
end

@polly function sumpool_backward(nn, nd, ih, iw, ow, oh, dh, dw, err_in, err_out)

    for n=1:nn, d=1:nd, r=1:oh, c=1:ow
        val=Float32(0.0)
        for h=10*r+1:min(10*r+dh,ih), w=10*c:min(10*c+dw,iw)
            err_in[n,d,h,w] += err_out[n,d,r,c]
        end
    end
end

let
    nn = 150
    nd = 120
    ih = 100
    iw = 100
    dh = 15
    dw = 15
    sh = 10
    sw = 10
    oh = div(ih - dh,sh) + 1
    ow = div(iw - dw,sw) + 1

    inp_F=zeros(Float32,nn,nd,ih,iw)
    out_F=zeros(Float32,nn,nd,oh,ow)
    err_in=zeros(Float32,nn,nd,ih,iw)
    err_out=zeros(Float32,nn,nd,oh,ow)

    #init_array

	for a=1:nn,b=1:nd,c=1:oh,d=1:ow
		out_F[a,b,c,d] = Float32((a*b + c*d) % nn)
		err_out[a,b,c,d] = Float32((a+b+c+d) % nn)
	end

	for a=1:nn,b=1:nd,c=1:iw,d=1:ih
		inp_F[a,b,c,d] = Float32((a*b + c*d) % nd)
		err_in[a,b,c,d] = Float32((a+b+c+d) % nd)
	end

	if length(ARGS)==0 || ARGS[1]=="forward"
		SUITE["sumpool"] = @benchmarkable sumpool_forward($nn,$nd,$ih,$iw,$ow,$oh,$dh,$dw,$sh,$sw,inp_F,out_F) setup = (inp_F=copy($inp_F); out_F=copy($out_F); err_in=copy($err_in); err_out=copy($err_out))
	else
		SUITE["sumpool"] = @benchmarkable sumpool_backward($nn,$nd,$ih,$iw,$ow,$oh,$dh,$dw, err_in, err_out)  setup = (inp_F=copy($inp_F); out_F=copy($out_F); err_in=copy($err_in); err_out=copy($err_out))
	end
end

max_stars_repo_name: akileshbadrinaaraayanan/IITH

max_issues_repo_name: akileshbadrinaaraayanan/IITH

max_forks_repo_name: akileshbadrinaaraayanan/IITH


text:
# using OMEinsum
using BackwardsLinalg

@doc raw"
    trg(a, χ, niter)

return the partition-function of a two-dimensional system of size `2^niter`
described by the tensor `a` calculated via the tensor renormalization group
algorithm.
`a` is a rank-4 tensor with the following indices:

        |1
    4--[a]--2
       3|
"
function trg(a::AbstractArray{T,4}, χ, niter; tol::Float64 = 1e-16) where T
    lnZ = zero(T)
    for n in 1:niter
        maxval = maximum(abs.(a))
        a /= maxval
        lnZ += 2.0^(1-n)*log(maxval)

        dr_ul = ein"urdl -> drul"(a)
        ld_ru = ein"urdl -> ldru"(a)
        dr, ul = trg_svd(dr_ul, χ, tol)
        ld, ru = trg_svd(ld_ru, χ, tol)

        a = ein"npu,por,dom,lmn -> urdl"(dr,ld,ul,ru)
    end
    trace = ein"ijij -> "(a)[]
    lnZ += log(trace)/2.0^niter
    return lnZ
end

function trg_svd(t, dmax, tol)
    d1, d2, d3, d4 = size(t)
    tmat = reshape(t, d1*d2, d3*d4)
    u, s, v = svd(tmat)
    dmax = min(searchsortedfirst(s, tol, rev=true), dmax, length(s))
    FS = s[1:dmax]
    sqrtFSp = sqrt.(FS)
    u = reshape(ein"ij,j -> ij"(u[:,1:dmax],  sqrtFSp), (d1, d2, dmax))
    v = reshape(ein"ij,i -> ij"(copy(v')[1:dmax,:], sqrtFSp), (dmax, d3, d4))

    return u, v
end


text:
include("basic.jl")
include("tracer_provider.jl")


text:
"""
    Axis(k[, v=OneTo(length(k))])

Subtypes of `AbstractAxis` that maps keys to values. The first argument specifies
the keys and the second specifies the values. If only one argument is specified
then the values span from 1 to the length of `k`.

## Examples

The value for all of these is the same.
```jldoctest axis_examples
julia> using AxisIndices

julia> x = Axis(2.0:11.0)  # when only one argument is specified assume it's the keys
Axis(2.0:1.0:11.0 => SimpleAxis(1:10))

julia> y = Axis(1:10)
Axis(1:10 => SimpleAxis(1:10))
```

Standard indexing returns the same values
```jldoctest axis_examples
julia> x[2]
2

julia> x[2] == y[2]
true

julia> x[1:2]
Axis(2.0:1.0:3.0 => SimpleAxis(1:2))

julia> y[1:2]
Axis(1:2 => SimpleAxis(1:2))

julia> x[1:2] == y[1:2]
true
```

Functions that return `true` or `false` may be used to search the keys for their
corresponding index. The following is equivalent to the previous example.
```julia
julia> x[==(3.0)]
2

julia> x[3.0] ==  # 3.0 is the 2nd key of x
       y[==(2)]   # 2 is the 2nd key of z
true

julia> x[<(4.0)]  # all keys less than 4.0 are 2.0:3.0 which correspond to values 1:2
Axis(2.0:1.0:3.0 => SimpleAxis(1:2))

julia> y[<=(3.0)]  # all keys less than or equal to 3.0 are 2.0:3.0 which correspond to values 1:2
Axis(2.0:1.0:3.0 => SimpleAxis(1:2))

julia> z[<(3)]  # all keys less than or equal to 3 are 1:2 which correspond to values 1:2
Axis(1:2 => SimpleAxis(1:2))

julia> x[<(4.0)] == y[<=(3.0)] == z[<(3)]
true
```
Notice that `==` returns a single value instead of a collection of all elements
where the key was found to be true. This is because all keys must be unique so
there can only ever be one element returned.
"""
struct Axis{K,I,Ks,Inds<:AbstractRange{I}} <: AbstractAxis{I,Inds}
    keys::Ks
    parent::Inds

    function Axis{K,I,Ks,Inds}(ks::Ks, inds::Inds) where {K,I,Ks<:AbstractVector{K},Inds}
        check_axis_length(ks, inds)
        check_unique_keys(k)
        return new{K,I,Ks,Inds}(ks, inds)
    end

    function Axis{K,V,Ks,Vs}(x::AbstractUnitRange{<:Integer}) where {K,V,Ks,Vs}
        if x isa Ks
            if x isa Vs
                return new{K,V,Ks,Vs}(x, x)
            else
                return new{K,V,Ks,Vs}(x, Vs(x))
            end
        else
            if x isa Vs
                return new{K,V,Ks,Vs}(Ks(x), x)
            else
                return  new{K,V,Ks,Vs}(Ks(x), Vs(x))
            end
        end
    end

    function Axis{K,I,Ks,Inds}(axis::SimpleAxis) where {K,I,Ks,Inds}
        return new{K,I,Ks,Inds}(axis, axis)
    end

    function Axis{K,I,Ks,Inds}(axis::AbstractAxis) where {K,I,Ks,Inds}
        return new{K,I,Ks,Inds}(Ks(keys(axis)), Inds(parent(axis)))
    end

    function Axis{K,I,Ks,Inds}(axis::Axis{K,I,Ks,Inds}) where {K,I,Ks,Inds}
        if can_change_size(axis)
            return copy(axis)
        else
            return axis
        end
    end

    # Axis{K,I}
    function Axis{K,I}() where {K,I}
        ks = Vector{K}()
        inds = SimpleAxis()
        return new{K,I,typeof(ks),typeof(inds)}()
    end

    function Axis{K,I}(ks::AbstractVector) where {K,I}
        check_unique_keys(ks)
        if can_change_size(ks)
            inds = SimpleAxis(OneToMRange{I}(length(ks)))
        else
            inds = SimpleAxis(indices(ks))
        end
        return new{K,I,typeof(ks),typeof(inds)}()
    end
    function Axis{K,I}(ks::AbstractVector, inds::AbstractAxis) where {K,I}
        if eltype(ks) <: K
            if eltype(inds) <: I
                return Axis{K,I,typeof(ks),typeof(inds)}(ks, inds)
            else
                return Axis{K,I}(ks, AbstractUnitRange{I}(inds))
            end
        else
            return Axis{K,I}(AbstractVector{K}(ks), inds)
        end
    end
    function Axis{K,I}(ks::AbstractVector, inds::AbstractUnitRange) where {K,I}
        return Axis{K,I}(ks, compose_axis(inds))
    end

    # Axis
    function Axis(axis::AbstractAxis)
        if can_change_size(axis)
            return axis
        else
            return copy(axis)
        end
    end

    Axis{K}() where {K} = Axis{K,Int}()

    Axis() = Axis{Any}()

    Axis(x::Pair) = Axis(x.first, x.second)

    function Axis(ks::AbstractVector, inds::AbstractAxis)
        return new{eltype(ks),eltype(inds),typeof(ks),typeof(inds)}(ks, inds)
    end

    Axis(ks::AbstractVector, inds::AbstractUnitRange) = Axis(ks, compose_axis(inds))

    function Axis(ks::AbstractVector)
        check_unique_keys(ks)
        if can_change_size(ks)
            inds = SimpleAxis(OneToMRange(length(ks)))
        else
            inds = compose_axis(static_first(eachindex(ks)):static_length(ks))
        end
        new{eltype(ks),eltype(inds),typeof(ks),typeof(inds)}(ks, inds)
    end
end

initialize_axis(ks, inds) = initialize_axis(ks, compose_axis(inds))
function initialize_axis(ks, inds::AbstractAxis)
    return unsafe_initialize(
        Axis{eltype(ks),eltype(inds),typeof(ks),typeof(inds)},
        (ks, inds)
    )
end

## interface
Base.keys(axis::Axis) = getfield(axis, :keys)
@inline Base.getproperty(axis::Axis, k::Symbol) = getproperty(parent(axis), k)

function ArrayInterface.unsafe_reconstruct(axis::Axis{K,I,Ks,Inds}, inds; keys=nothing) where {K,I,Ks,Inds}
    if keys === nothing
        ks = Base.keys(axis)
        p = parent(axis)
        kindex = firstindex(ks)
        pindex = first(p)
        if kindex === pindex
            return initialize_axis(@inbounds(ks[inds]), inds)
        else
            return initialize_axis(@inbounds(ks[inds .+ (pindex - kindex)]), inds)
        end
    else
        return initialize_axis(keys, inds)
    end
end

@inline function ArrayInterface.to_axis(::IndexStyle, axis::Axis, inds)
    if allunique(inds)
        ks = Base.keys(axis)
        p = parent(axis)
        kindex = firstindex(ks)
        pindex = first(p)
        if kindex === pindex
            return initialize_axis(@inbounds(ks[inds]), to_axis(parent(axis), inds))
        else
            return initialize_axis(
                @inbounds(ks[inds .+ (pindex - kindex)]),
                to_axis(parent(axis), inds)
            )
        end
    else
        return unsafe_reconstruct(axis, to_axis(parent(axis), inds))
    end
end

## other stuff
function StaticRanges.set_last!(axis::Axis, val)
    can_set_last(axis) || throw(MethodError(set_last!, (axis, val)))
    set_last!(parent(axis), val)
    resize_last!(keys(axis), length(parent(axis)))
    return axis
end

function StaticRanges.set_last(axis::Axis, val)
    vs = set_last(parent(axis), val)
    return unsafe_reconstruct(axis, resize_last(keys(axis), length(vs)); keys=vs)
end

function StaticRanges.set_first(axis::Axis, val)
    vs = set_first(parent(axis), val)
    return unsafe_reconstruct(axis, vs; keys=resize_first(keys(axis), length(vs)))
end

function StaticRanges.set_first!(axis::Axis, val)
    can_set_first(axis) || throw(MethodError(set_first!, (axis, val)))
    set_first!(parent(axis), val)
    resize_first!(keys(axis), length(parent(axis)))
    return axis
end

function maybe_unsafe_reconstruct(axis::Axis, inds::AbstractUnitRange{I}; keys=nothing) where {I<:Integer}
    if keys === nothing
        return unsafe_reconstruct(axis, SimpleAxis(inds); keys=@inbounds(Base.keys(axis)[inds]))
    else
        return unsafe_reconstruct(axis, SimpleAxis(inds); keys=keys)
    end
end
function maybe_unsafe_reconstruct(axis::Axis, inds::AbstractArray)
    if keys === nothing
        axs = (unsafe_reconstruct(axis, SimpleAxis(eachindex(inds))),)
    elseif allunique(inds)
        axs = (unsafe_reconstruct(axis, SimpleAxis(eachindex(inds)); keys=@inbounds(keys(axis)[inds])),)
    else  # not all indices are unique so will result in non-unique keys
        axs = (SimpleAxis(eachindex(inds)),)
    end
    return AxisArray{eltype(axis),ndims(inds),typeof(inds),typeof(axs)}(inds, axs)
end

@inline function Base.length(axis::Axis{K,I,Ks,Inds}) where {K,I,Ks,Inds}
    if known_length(Ks) === nothing
        return length(parent(axis))
    else
        return known_length(Ks)
    end
end

function StaticRanges.set_length!(axis::Axis, len)
    can_set_length(axis) || error("Cannot use set_length! for instances of typeof $(typeof(axis)).")
    set_length!(parent(axis), len)
    set_length!(keys(axis), len)
    return axis
end

function StaticRanges.set_length(axis::Axis, len)
    return unsafe_reconstruct(
        axis,
        set_length(parent(axis), len);
        keys=set_length(keys(axis), len)
    )
end

for f in (:grow_last!, :grow_first!, :shrink_last!, :shrink_first!)
    @eval begin
        function StaticRanges.$f(axis::Axis, n::Integer)
            can_set_length(axis) ||  throw(MethodError($f, (axis, n)))
            StaticRanges.$f(keys(axis), n)
            StaticRanges.$f(parent(axis), n)
            return axis
        end
    end
end

for f in (:grow_last, :grow_first, :shrink_last, :shrink_first, :resize_first, :resize_last)
    @eval begin
        @inline function StaticRanges.$f(axis::Axis, n::Integer)
            return unsafe_reconstruct(
                axis,
                StaticRanges.$f(parent(axis), n);
                keys = StaticRanges.$f(keys(axis), n)
            )
        end

    end
end

@inline function StaticRanges.shrink_last(axis::Axis, n::AbstractUnitRange{<:Integer})
    return unsafe_reconstruct(axis, n; keys=shrink_last(keys(axis), length(axis) - length(n)))
end
@inline function StaticRanges.shrink_first(axis::Axis, n::AbstractUnitRange{<:Integer})
    return unsafe_reconstruct(axis, n; keys=shrink_first(keys(axis), length(axis) - length(n)))
end
function StaticRanges.grow_last(axis::Axis, n::AbstractUnitRange{<:Integer})
    return unsafe_reconstruct(axis, n; keys=grow_last(keys(axis), length(n) - length(axis)),)
end
function StaticRanges.grow_first(axis::Axis, n::AbstractUnitRange{<:Integer})
    return unsafe_reconstruct(axis, n;keys=grow_first(keys(axis), length(n) - length(axis)))
end
function StaticRanges.resize_last(axis::Axis, n::AbstractUnitRange{<:Integer})
    return unsafe_reconstruct(axis, n; keys=resize_last(keys(axis), length(n)))
end
function StaticRanges.resize_first(axis::Axis, n::AbstractUnitRange{<:Integer})
    return unsafe_reconstruct(axis, n; keys=resize_first(keys(axis), length(n)))
end

function Base.pop!(axis::Axis)
    can_set_last(axis) || throw(MethodError(pop!, axis))
    pop!(keys(axis))
    return pop!(parent(axis))
end

function Base.popfirst!(axis::Axis)
    can_set_last(axis) || throw(MethodError(pop!, axis))
    popfirst!(keys(axis))
    return popfirst!(parent(axis))
end

function push_key!(axis::Axis, key)
    push!(keys(axis), key)
    grow_last!(parent(axis), 1)
    return nothing
end

function popfirst_axis!(axis::Axis)
    if StaticRanges.can_set_first(axis)
        StaticRanges.shrink_first!(keys(axis), 1)
    else
        shrink_last!(keys(axis), 1)
    end
    shrink_last!(parent(axis), 1)
    return nothing
end

_keys_type(::Type{T}) where {Ks,T<:Axis{<:Any,<:Any,Ks,<:Any}} = Ks

Base.keytype(::Type{T}) where {K,I,T<:Axis{K,I}} = K
function StaticRanges.can_set_length(::Type{T}) where {K,I,Ks,Inds,T<:Axis{K,I,Ks,Inds}}
    return can_set_length(Ks) & can_set_length(Inds)
end
function StaticRanges.can_set_last(::Type{T}) where {K,I,Ks,Inds,T<:Axis{K,I,Ks,Inds}}
    return can_set_last(Ks) & can_set_last(Inds)
end
function StaticRanges.can_set_first(::Type{T}) where {K,I,Ks,Inds,T<:Axis{K,I,Ks,Inds}}
    return can_set_first(Ks) & can_set_first(Inds)
end

@propagate_inbounds function Base.getindex(axis::Axis, arg::AbstractUnitRange{I}) where {I<:Integer}
    @boundscheck checkbounds(axis, arg)
    ks = Base.keys(axis)
    p = parent(axis)
    kindex = firstindex(ks)
    pindex = first(p)
    if kindex === pindex
        return initialize_axis(
            @inbounds(ks[arg]),
            @inbounds(getindex(p, arg))
        )
    else
        return initialize_axis(
            @inbounds(ks[arg .+ (kindex - pindex)]),
            @inbounds(getindex(p, arg))
        )
    end
end

print_axis(io::IO, axis::Axis) = print(io, "Axis($(keys(axis)) => $(parent(axis)))")


text:
using Test

@testset "FFI unit tests" begin

    @testset "Dataset" begin
        include(joinpath("ffi", "datasets.jl"))
    end

    @testset "Booster" begin
        include(joinpath("ffi", "booster.jl"))
    end

    @testset "Base Utils" begin
        include(joinpath("ffi", "base_utils.jl"))
    end

end

@testset "MLJ interface tests" begin

    @testset "Parameters" begin
        include(joinpath("mlj", "parameters.jl"))
    end

    @testset "Binary LightGBM" begin
        include(joinpath("mlj", "binary_classifier.jl"))
    end

    @testset "Multiclass LightGBM" begin
        include(joinpath("mlj", "multiclass_classifier.jl"))
    end

    @testset "Regression LightGBM" begin
        include(joinpath("mlj", "regression.jl"))
    end

    @testset "MLJ update interface" begin
        include(joinpath("mlj", "update.jl"))
    end

    @testset "MLJ update interface" begin
        include(joinpath("mlj", "user_report.jl"))
    end

end

@testset "Basic tests" begin

    @testset "Estimator parameters" begin
        include(joinpath("basic", "parameters.jl"))
    end

    @testset "Estimator parameters" begin
        include(joinpath("basic", "evaluation_metrics.jl"))
    end

    @testset "Utils" begin
        include(joinpath("basic", "utils.jl"))
    end

    @testset "Fit" begin
        include(joinpath("basic", "fit.jl"))
    end

end

@testset "Integration tests" begin

    @testset "(OLD) Basic Tests" begin
        include("basic_tests.jl")
    end

    @testset "Weights Test" begin
        include("weightsTest.jl")
    end

    @testset "Init Score Test" begin
        include("initScoreTest.jl")
    end

end