text:
"""
A DB environment supports multiple databases, all residing in the same shared-memory map.
"""
type Environment
    handle::Ptr{Void}
    path::String
    Environment() = new(C_NULL, "")
    Environment(h::Ptr{Void}) = new(h, "")
end

"Return the path that was used in `open`"
path(env::Environment) = env.path

"Check if environment is open"
isopen(env::Environment) = env.handle != C_NULL

"Create an LMDB environment handle"
function create()
    env_ref = Ref{Ptr{Void}}(C_NULL)
    ret = ccall( (:mdb_env_create, liblmdb), Cint, (Ptr{Ptr{Void}},), env_ref)
    (ret != 0) && throw(LMDBError(ret))
    return Environment(env_ref[])
end

"Wrapper of `create` for `do` construct"
function create(f::Function)
    env = create()
    try
        f(env)
    finally
        close(env)
    end
end

"""Open an environment handle

`open` function accepts folowing parameters:
* `env` db environment object
* `path` directory in which the database files reside
* `flags` defines special options for the environment
* `mode` UNIX permissions to set on created files

*Note:* A database directory must exist and be writable.
"""
function open(env::Environment, path::String; flags::Cuint=zero(Cuint), mode::Cmode_t = 0o755)
    env.path = path
    cpath = bytestring(path)
    ret = ccall((:mdb_env_open, liblmdb), Cint,
                 (Ptr{Void}, Cstring, Cuint, Cmode_t),
                  env.handle, cpath, flags, mode)
    (ret != 0) && throw(LMDBError(ret))
    return ret::Cint
end

"""Close the environment and release the memory map"""
function close(env::Environment)
    if env.handle == C_NULL
        throw(LMDBError(-1,"Environment is already closed"))
    end
    ccall( (:mdb_env_close, liblmdb), Void, (Ptr{Void},), env.handle)
    env.handle = C_NULL
    env.path = ""
    return
end

"""Flush the data buffers to disk"""
function sync(env::Environment, force::Bool = false)
    fval = force ? 1 : 0
    ret = ccall( (:mdb_env_sync, liblmdb), Cint, (Ptr{Void}, Cint), env.handle, fval)
    (ret != 0) && throw(LMDBError(ret))
    return ret::Cint
end

"""Set environment flags"""
function set!(env::Environment, flag::Cuint)
    ret = ccall( (:mdb_env_set_flags, liblmdb), Cint, (Ptr{Void}, Cuint, Cint), env.handle, flag, one(Cint))
    (ret != 0) && throw(LMDBError(ret))
    return flag
end
set!(env::Environment, flag::EnvironmentFlags) = set!(env, Cuint(flag))

"""Unset environment flags"""
function unset!(env::Environment, flag::Cuint)
    ret = ccall( (:mdb_env_set_flags, liblmdb), Cint, (Ptr{Void}, Cuint, Cint), env.handle, flag, zero(Cint))
    (ret != 0) && throw(LMDBError(ret))
    return flag
end
unset!(env::Environment, flag::EnvironmentFlags) = unset!(env, Cuint(flag))

"""Set environment flags and parameters

`setindex!` accepts folowing parameters:
* `env` db environment object
* `option` symbol which indicates parameter. Currently supported parameters:
    * Flags
    * Readers
    * MapSize
    * DBs
* `value` parameter value

**Note:** Consult LMDB documentation for particual values of environment parameters and flags.
"""
function setindex!(env::Environment, val::Cuint, option::Symbol)
    ret = zero(Cint)
    if option == :Readers
        ret = ccall( (:mdb_env_set_maxreaders, liblmdb), Cint, (Ptr{Void}, Cuint), env.handle, val)
    elseif option == :MapSize
        ret = ccall( (:mdb_env_set_mapsize, liblmdb), Cint, (Ptr{Void}, Cuint), env.handle, val)
    elseif option == :DBs
        ret = ccall( (:mdb_env_set_maxdbs, liblmdb), Cint, (Ptr{Void}, Cuint), env.handle, val)
    else
        warn("Cannot set $(string(option)) value")
    end
    (ret != 0) && throw(LMDBError(ret))
    return ret
end
setindex!(env::Environment, val::Int, option::Symbol) = setindex!(env, Cuint(val), option)

"""Get environment flags and parameters

`getindex` accepts folowing parameters:
* `env` db environment object
* `option` symbol which indicates parameter. Currently supported parameters:
    * Flags
    * Readers
    * KeySize

**Note:** Consult LMDB documentation for particual values of environment parameters and flags.
"""
function getindex(env::Environment, option::Symbol)
    value = Cuint[0]
    if option == :Flags
        flags = Cuint[0]
        ret = ccall( (:mdb_env_get_flags, liblmdb), Cint, (Ptr{Void}, Ptr{Cuint}), env.handle, value)
        (ret != 0) && throw(LMDBError(ret))
    elseif option == :Readers
        ret = ccall( (:mdb_env_get_maxreaders, liblmdb), Cint, (Ptr{Void}, Ptr{Cuint}), env.handle, value)
        (ret != 0) && throw(LMDBError(ret))
    elseif option == :KeySize
        value[1] = ccall( (:mdb_env_get_maxkeysize, liblmdb), Cint, (Ptr{Void}, ), env.handle)
    else
        warn("Cannot get $(string(option)) value")
    end
    return value[1]
end

"""Information about the environment"""
immutable EnvironmentInfo
    mapaddr::Ptr{Void}
    mapsize::Csize_t     # Size of the data memory map
    last_pgno::Csize_t   # ID of the last used page
    last_txnid::Csize_t  # ID of the last committed transaction
    maxreaders::Cuint    # max reader slots in the environment */
    numreaders::Cuint    # max reader slots used in the environment */
    EnvironmentInfo() = new(C_NULL, zero(Csize_t), zero(Csize_t), zero(Csize_t), zero(Cuint), zero(Cuint))
end

"""Return information about the LMDB environment."""
function info(env::Environment)
    ei_ref = Ref(EnvironmentInfo())
    !isopen(env) && return ei_ref[]
    ret = ccall( (:mdb_env_info, liblmdb), Cint, (Ptr{Void}, Ptr{EnvironmentInfo}), env.handle, ei_ref)
    (ret != 0) && throw(LMDBError(ret))
    return ei_ref[]
end

function show(io::IO, env::Environment)
    print(io,"Environment is ", isopen(env) ? (isempty(env.path) ? "created" : "opened") : "closed")
    if !isempty(env.path)
        print(io,"\nDB path: $(path(env))")
        ei = info(env)
        print(io,"\nSize of the data memory map: $(ei.mapsize)")
        print(io,"\nID of the last used page: $(ei.last_pgno)")
        print(io,"\nID of the last committed transaction: $(ei.last_txnid)")
        print(io,"\nMax reader slots in the environment: $(ei.maxreaders)")
        print(io,"\nMax reader slots used in the environment: $(ei.numreaders)")
    end
end


text:
function start_simulation(dataset_path::String, output_path::String, parameters_path::String, info_path::String)
	# Create the check-in writer
	event_writer = EventWriter(output_path)
	#Get parameters
	parameters = JSON.parsefile(parameters_path)
	info = JSON.parsefile(info_path)

	# Define simulation start time
	simulation_time_start = DateTime(2021,1,1,DAY_START)
	simulation_time_end = simulation_time_start + Dates.Day(31) + Dates.Day(31)
	households = _read_household_dataset(dataset_path)

	# Available because workplaces always starts at id 1
	leisure_places = rand(info[WORKPLACES_START_ID]:info[WORKPLACES_START_ID] + parameters["workplaces"],
									 parameters["leisures"])
	serialize("resources/datasets/Salerno/leisure_places.bin",leisure_places)
	count = 0
	while simulation_time_start < simulation_time_end
		count = 0
		for household in households
			for worker in household.workers
				_simulate_worker(worker,household.id,parameters,info,simulation_time_start,event_writer,leisure_places)
				count+=1
			end
			for student in household.students
				_simulate_student(student,household.id,parameters,info,simulation_time_start,event_writer,leisure_places)
				count+=1
			end
			for other in household.unemployeds
				_simulate_home(other,household.id,parameters,info,simulation_time_start,event_writer,leisure_places)
				count+=1
			end
			write_events(event_writer)

		end

		simulation_time_start += Dates.Hour(1)
		if Dates.hour(simulation_time_start) ==  0
			simulation_time_start += Dates.Hour(DAY_START)
		end
	end

	close_writer(event_writer)
end

function _simulate_worker(worker::Person, household_id::Int,parameters,info,time,event_writer::EventWriter,leisure_places)
	# Handle morning
	if worker.movement_state == 1
		add_event(event_writer,Event(worker.id,household_id,PLACE_HOME,time - Dates.Hour(1))) # They wake up 1 hour before work
	end

	if _worker_next_step(time,worker)
		if worker.movement_state == 1 # Home
			event = Event(worker.id,household_id,PLACE_HOME,time)
		end
		if worker.movement_state == 2 # Transport
			event = Event(worker.id,rand(info[TRANSPORTS_START_ID]:info[WORKPLACES_START_ID] + parameters["transports"]),PLACE_TRANSPORT,time)
		end
		if worker.movement_state == 3 # Work
			# Introduce an event each 4 hours
			local_time = time
			while Dates.hour(local_time) < WORK_END
				add_event(event_writer, Event(worker.id,worker.work_school_place,PLACE_WORK,local_time))
				local_time += Dates.Hour(4) # Time for indirect contact
			end
			event = Event(-1,-1,-1,DateTime(0000)) #FIXME: Null event, this can be easily avoided
		end
		if worker.movement_state == 4 # Transport
			event = Event(worker.id,rand(info[TRANSPORTS_START_ID]:info[TRANSPORTS_START_ID] + parameters["transports"]),PLACE_TRANSPORT,time)
		end
		if worker.movement_state == 5 # Home
			event = Event(worker.id,household_id,PLACE_HOME,time)
		end
		if worker.movement_state == 6 # leisure
			leisure_place = rand(leisure_places)
			event = Event(worker.id,leisure_place,PLACE_LEISURE,time)
			# event = Event(-1,-1,-1,DateTime(0000)) #FIXME: Null event, this can be easily avoided
		end
		add_event(event_writer,event)
	end
end

function _simulate_student(student::Person, household_id::Int,parameters,info,time,event_writer::EventWriter,leisure_places)
	# Handle morning
	if student.movement_state == 1
		add_event(event_writer,Event(student.id,household_id,PLACE_HOME,time - Dates.Hour(1)))
	end

	if _student_next_step(time,student)
		if student.movement_state == 1 # Home
			event = Event(student.id,household_id,PLACE_HOME,time)
		end
		if student.movement_state == 2 # Transport
			event = Event(student.id,rand(info[TRANSPORTS_START_ID]:info[TRANSPORTS_START_ID] + parameters["transports"]),PLACE_TRANSPORT,time)
		end
		if student.movement_state == 3 # School
			event = Event(student.id,student.work_school_place,PLACE_WORK,time)
		end
		if student.movement_state == 4 # Transport
			event = Event(student.id,rand(info[TRANSPORTS_START_ID]:info[TRANSPORTS_START_ID] + parameters["transports"]),PLACE_TRANSPORT,time)
		end
		if student.movement_state == 5 # Home
			event = Event(student.id,household_id,PLACE_HOME,time)
		end
		if student.movement_state == 6 # leisure
			leisure_place = rand(leisure_places)
			event = Event(student.id,leisure_place,PLACE_LEISURE,time)
			# event = Event(-1,-1,-1,DateTime(0000)) #FIXME: Null event, this can be easily avoided
		end
		add_event(event_writer,event)
	end
end

function _simulate_home(home::Person, household_id::Int,parameters,info,time,event_writer::EventWriter,leisure_places)
	if home.movement_state == 1
		add_event(event_writer,Event(home.id,household_id,PLACE_HOME,time - Dates.Hour(1)))
	end

	if _home_next_step(time,home)
		if home.movement_state == 1 # Home
			event = Event(home.id,household_id,PLACE_HOME,time)
		end
		if home.movement_state == 2 # leisure
			leisure_place = rand(leisure_places)
			event = Event(home.id,leisure_place,PLACE_LEISURE,time)
			# event = Event(-1,-1,-1,DateTime(0000)) #FIXME: Null event, this can be easily avoided
		end
		add_event(event_writer,event)
	end
end

# function get_households(df_nodes::DataFrame,parameters)::Array{Household,1}
# 	households = Array{Household,1}()
# 	gdf_households = groupby(df_nodes, :hh_id)

# 	lk = ReentrantLock()
# 	Threads.@threads for i in 1:length(gdf_households)
# 		household_record = gdf_households[i]
# 		id = household_record.part_id[1]
# 		household = Household(id,size(household_record,1))
# 		for person_record in eachrow(household_record)
# 			occupation = person_record.part_occupation
# 			if occupation == WORKING
# 				person = create_worker(parameters,person_record)
# 				push!(household.workers,person)
# 			elseif occupation == SCHOOL
# 				person = create_student(parameters,person_record)
# 				push!(household.students,person)
# 			else
# 				person = Person(person_record.part_id,person_record.part_age,person_record.part_occupation,-1,-1)
# 				push!(household.unemployeds,person)
# 			end
# 		end
# 		lock(lk) do
# 			push!(households,household)
# 		end
# 	end
# 	return households
# end

# function get_households(dataset_path::String,parameters::Contact_simulation_options)::Array{Household,1}
# 	df_nodes = CSV.File(dataset_path) |> DataFrame
# 	return get_households(df_nodes,parameters)
# end

# function create_student(parameters::Contact_simulation_options,student_record)
# 	school_place = rand(1:parameters.school_place_num)
# 	transport = rand(1:parameters.transport_place_num)
# 	return Person(student_record.part_id, student_record.part_age, student_record.part_occupation, school_place, transport)
# end

# function create_worker(parameters::Contact_simulation_options,worker_record)
# 	work_place = rand(1:parameters.work_place_num)
# 	transport = rand(1:parameters.transport_place_num)
# 	return Person(worker_record.part_id, worker_record.part_age, worker_record.part_occupation, work_place, transport)
# end

ext: jl

max_stars_repo_name: Luigi-Crisci/Hypergraphs-Epidemics

max_issues_repo_name: Luigi-Crisci/Hypergraphs-Epidemics

max_forks_repo_name: Luigi-Crisci/Hypergraphs-Epidemics


text:
res1 = get_interval(
    [3.0],
    (x)->x[1]^2,
    f_1p,
    :CICO_ONE_PASS;
    scan_bounds=(1e-9,1e9),
    #theta_bounds = [(1e-2,1e2)],
    loss_crit = 8.,
    local_alg = :LN_NELDERMEAD,
    silent = true
)
@test isapprox(res1.result[1].value, (3-sqrt(3))^2,atol=1e-2)
@test isapprox(res1.result[2].value, (3+sqrt(3))^2,atol=1e-2)

res2 = get_interval(
    [3.0],
    (x)->x[1]^2,
    f_1p,
    :CICO_ONE_PASS;
    scan_bounds=(1e-9,1e9),
    #theta_bounds = [(1e-2,1e2)],
    loss_crit = 8.,
    local_alg = :LD_MMA,
    silent = true,
    loss_grad = :AUTODIFF,
    scan_grad = :AUTODIFF
)
@test isapprox(res2.result[1].value, (3-sqrt(3))^2,atol=1e-2)
@test isapprox(res2.result[2].value, (3+sqrt(3))^2,atol=1e-2)

res3 = get_interval(
    [3.0],
    (x) -> log10(x[1]^2),
    f_1p,
    :CICO_ONE_PASS,
    scan_bounds=(1e-9,1e9),
    #theta_bounds = [(1e-2,1e2)],
    loss_crit = 8.,
    local_alg = :LD_MMA,
    silent = true,
    loss_grad = :AUTODIFF,
    scan_grad = :AUTODIFF
)
@test isapprox(res3.result[1].value, 2*log10(3-sqrt(3)),atol=1e-2)
@test isapprox(res3.result[2].value, 2*log10(3+sqrt(3)),atol=1e-2)


text:
## ----------  ObjectId

# """
#     SingleId

# Id for one object or object vector.
# Does not keep track of location in model.
# Keeps track of `name` (Symbol identifier), `index`, and `description` (for display purposes).

# `index` is used when there is a vector or matrix of objects of the same type
# `index` is typically empty (scalar object) or scalar (vector objects)
# """
# struct SingleId
#     name :: Symbol
#     index :: Array{Int}
#     description :: String
# end

# """
#     ObjectId

# Complete, unique ID of a `ModelObject`

# Contains own id and a vector of parent ids, so one knows exactly where the object
# is placed in the model tree.
# """
# struct ObjectId
#     # Store IDs as vector, not tuple (b/c empty tuples are tricky)
#     # "Youngest" member is positioned last in vector
#     ids :: Vector{SingleId}
# end

# -------------


text:
using FFTW
using LinearAlgebra
using Statistics

export DistributionFunction

"""
    DistributionFunction( grid1, grid2 )

"""
struct DistributionFunction

    xgrid :: OneDGrid
    vgrid :: OneDGrid
    values :: Array{AbstractFloat, 2}

    function DistributionFunction( xgrid, vgrid )

        nx = xgrid.len
        nv = vgrid.len
        values = zeros(nx, nv)

        new( xgrid, vgrid, values )

    end

end

export landau!

"""
    landau!( f, α, kx )

Initialize the distribution function for the Landau damping test case

```math
f(x,v,t=0) = \\frac{1}{\\sqrt{2\\pi}} ( 1 + \\cos (kx x) ) \\exp (-\\frac{v^2}{2})
```

"""
function landau!( f :: DistributionFunction, α, kx)

    nx = f.xgrid.len
    nv = f.vgrid.len
    x  = f.xgrid.points
    v  = f.vgrid.points
    f.values .= (1 .+ α*cos.(kx*x))/sqrt(2π) .* transpose(exp.(-0.5*v.^2))

end

export two_stream_instability!

function two_stream_instability!(f; eps = 0.01, xi = 0.90, v0 = 2.4)

    nx = f.xgrid.len
    nv = f.vgrid.len
    xg = f.xgrid.points
    vg = f.vgrid.points

    for (i,x) in enumerate(xg), (j,v) in enumerate(vg)
        f.values[i,j] = ((1 + eps*((cos(4π*x) + cos(6π*x)) / 1.2
                .+ cos(2π*x))) * (1/sqrt(2π)) * ((2-2xi)/(3-2xi))
                    * (1 + 5 * v^2 / (1-xi)) * exp(-.5 * v^2))
    end
end

"""
    bspline(p, j, x)

Return the value at x in [0,1] of the B-spline with integer nodes of degree p with support starting at j.
Implemented recursively using the [De Boor's Algorithm](https://en.wikipedia.org/wiki/De_Boor%27s_algorithm)

```math
B_{i,0}(x) := \\left\\{
\\begin{matrix}
1 & \\mathrm{if}  \\quad t_i ≤ x < t_{i+1} \\\\
0 & \\mathrm{otherwise}
\\end{matrix}
\\right.
```

```math
B_{i,p}(x) := \\frac{x - t_i}{t_{i+p} - t_i} B_{i,p-1}(x)
+ \\frac{t_{i+p+1} - x}{t_{i+p+1} - t_{i+1}} B_{i+1,p-1}(x).
```
"""
function bspline(p::Int, j::Int, x::Float64)
   if p == 0
      j == 0 ? (return 1.0) : (return 0.0)
   else
      w = (x - j) / p
      w1 = (x - j - 1) / p
   end
   return w * bspline(p - 1, j, x) + (1 - w1) * bspline(p - 1, j + 1, x)
end

"""
    advection!(f, grid, v, dt; p = 5)

Advect the distribution function `f` with velocity `v` along first `f` dimension
with a time step `dt`. Interpolation method uses bspline periodic of order 5 by default. Real type version.
"""
function advection!(f    :: Array{AbstractFloat, 2},
                    grid :: OneDGrid,
                    v,
                    dt;
                    p = 5)

   nx = grid.len
   nv = length(v)
   dx = grid.step
   modes = [2π * i / nx for i in 0:nx-1]

   # compute eigenvalues of degree p b-spline matrix
   eig_bspl  = zeros(Float64, nx)
   eig_bspl .= bspline(p, -div(p+1,2), 0.0)
   for i in 1:div(p+1,2)-1
      eig_bspl .+= bspline(p, i - (p+1)÷2, 0.0) * 2 .* cos.(i * modes)
   end
   eigalpha = zeros(Complex{Float64}, nx)

   ft = fft(f,1)

   for j in 1:nv
      alpha = dt * v[j] / dx
      # compute eigenvalues of cubic splines evaluated
      # at displaced points
      ishift = floor(-alpha)
      beta   = -ishift - alpha
      fill!(eigalpha, 0.0im)
      for i in -div(p-1,2):div(p+1,2)
         eigalpha .+= (bspline(p, i-div(p+1,2), beta)
                        .* exp.((ishift+i) * 1im .* modes))
      end

      # compute interpolating spline using fft and properties
      # of circulant matrices

      ft[:,j] .*= eigalpha ./ eig_bspl

   end

   f .= real(ifft(ft,1))

end

"""
    advection!(f, grid, v, dt; p = 5)

Advect the distribution function `f` with velocity `v` along first `f` dimension
with a time step `dt`. Interpolation method uses bspline periodic of order 5 by default. Complex type version.
"""
function advection!(f    :: Array{ComplexF64, 2},
                    grid :: OneDGrid,
                    v,
                    dt;
                    p = 5)

   nx = grid.len
   nv = length(v)
   dx = grid.step
   modes = [2π * i / nx for i in 0:nx-1]

   # compute eigenvalues of degree p b-spline matrix
   eig_bspl  = zeros(Float64, nx)
   eig_bspl .= bspline(p, -div(p+1,2), 0.0)
   for i in 1:div(p+1,2)-1
      eig_bspl .+= bspline(p, i - (p+1)÷2, 0.0) * 2 .* cos.(i * modes)
   end
   eigalpha = zeros(Complex{Float64}, nx)

   fft!(f,1)

   for j in 1:nv
      alpha = dt * v[j] / dx
      # compute eigenvalues of cubic splines evaluated
      # at displaced points
      ishift = floor(-alpha)
      beta   = -ishift - alpha
      fill!(eigalpha,0.0im)
      for i in -div(p-1,2):div(p+1,2)
         eigalpha .+= (bspline(p, i-div(p+1,2), beta)
                        .* exp.((ishift+i) * 1im .* modes))
      end

      # compute interpolating spline using fft and properties
      # of circulant matrices

      f[:,j] .*= eigalpha ./ eig_bspl

   end

   ifft!(f,1)

end

"""
    compute_rho(f)

Compute charge density
ρ(x,t) = ∫ f(x,v,t) dv
"""
function compute_rho( f::DistributionFunction)

   dv = f.vgrid.step
   rho = dv * sum(real(f.values), dims=2)
   vec(rho .- mean(rho)) # vec squeezes the 2d array returned by sum function

end

# + {"slideshow": {"slide_type": "slide"}}
"""
    compute_e(f)

compute Ex using that -ik*Ex = rho
"""
function compute_e( f::DistributionFunction )

   rho = compute_rho( f )
   nx = f.xgrid.len
   xmin = f.xgrid.start
   xmax = f.xgrid.stop
   k =  2π / (xmax - xmin)
   modes = k * vcat(0:div(nx,2)-1, -div(nx,2):-1)
   modes[1] = 1.0
   rhok = fft(rho) ./ modes .* (-1im)
   real(ifft!(rhok))

end

function advection_x!( f :: DistributionFunction, dt, p)

    advection!(f.values, f.xgrid, f.vgrid.points, dt; p)

end

function advection_v!( f :: DistributionFunction, dt, p)

    dx = f.xgrid.step
    nx = f.xgrid.len
    e  = compute_e(f)
    nrj = 0.5*log(sum(e.*e)*dx)
    fᵗ = collect(transpose(f.values))
    advection!(fᵗ, f.vgrid, e, dt; p)
    f.values .= transpose(fᵗ)
    nrj

end


text:
# This file is part of the IntervalArithmetic.jl package; MIT licensed

doc"""The `@interval` macro is the main method to create an interval.
It converts each expression into a narrow interval that is guaranteed to contain the true value passed by the user in the one or two expressions passed to it.
When passed two expressions, it takes the hull of the resulting intervals
to give a guaranteed containing interval.

Examples:
```
    @interval(0.1)

    @interval(0.1, 0.2)

    @interval(1/3, 1/6)

    @interval(1/3^2)
```
"""
macro interval(expr1, expr2...)
    make_interval(:(parameters.precision_type), expr1, expr2)
end

doc"The `@floatinterval` macro constructs an interval with `Float64` entries."
macro floatinterval(expr1, expr2...)
    make_interval(Float64, expr1, expr2)
end

doc"The `@biginterval` macro constructs an interval with `BigFloat` entries."
macro biginterval(expr1, expr2...)
    make_interval(BigFloat, expr1, expr2)
end

doc"""`transform` transforms a string by applying the function `f` and type
`T` to each argument, i.e. `:(x+y)` is transformed to `:(f(T, x) + f(T, y))`
"""
transform(x::Symbol, f, T) = :($f($T, $(esc(x))))   # use if x is not an expression
transform(x, f, T) = :($f($T, $x))   # use if x is not an

function transform(expr::Expr, f::Symbol, T)

    if expr.head in ( :(.), :ref )   # of form  a.lo  or  a[i]
        return :($f($T, $(esc(expr))))
    end

    new_expr = copy(expr)

    first = 1  # where to start processing arguments

    if expr.head == :call
        if expr.args[1] ∈ (:+, :-, :*, :/, :^)
            first = 2  # skip operator
        else   # escape standard function:
            new_expr.args[1] = :($(esc(expr.args[1])))
            first = 2
        end
    end

    if expr.head == :macrocall  # handles BigInts etc.
        return :($f($T, $(esc(expr))))  # hack: pass straight through
    end

    for (i, arg) in enumerate(expr.args)
        i < first && continue
        #@show i,arg

        new_expr.args[i] = transform(arg, f, T)
    end

    return new_expr
end

# Called by interval and floatinterval macros
doc"""`make_interval` does the hard work of taking expressions
and making each literal (0.1, 1, etc.) into a corresponding interval construction,
by calling `transform`."""

function make_interval(T, expr1, expr2)
    # @show expr1, expr2

    expr1 = transform(expr1, :atomic, :(Interval{$T}))

    if isempty(expr2)  # only one argument
        return :(Interval($expr1))
    end

    expr2 = transform(expr2[1], :atomic, :(Interval{$T}))

    :(interval(($expr1).lo, ($expr2).hi))
end


text:
"""
    MultiParallelNet(layers...)
    MultiParallelNet(layers::Tuple)
    MultiParallelNet(layers::Vector)

Creates a MultiParallelNet mostly used for MultiScale Models. It takes a list of inputs
and passes all of them through each `layer` and returns a tuple of outputs.

## Example

```
Model := MultiParallelNet(L1, L2, L3)

Model(X1, X2) := (Model.L1(X1, X2), Model.L2(X1, X2), Model.L3(X1, X2))
```
"""
struct MultiParallelNet{L}
    layers::L

    function MultiParallelNet(args...)
        layers = tuple(args...)
        return new{typeof(layers)}(layers)
    end

    MultiParallelNet(layers::Tuple) = new{typeof(layers)}(layers)

    MultiParallelNet(layers::Vector) = MultiParallelNet(layers...)
end

Flux.@functor MultiParallelNet

function (mpn::MultiParallelNet)(x::Union{Tuple,Vector})
    buf = Zygote.Buffer(Vector{Any}(undef, length(mpn.layers)))
    for (i, l) in enumerate(mpn.layers)
        buf[i] = l(x...)
    end
    return Tuple(copy(buf))
end

function (mpn::MultiParallelNet)(args...)
    buf = Zygote.Buffer(Vector{Any}(undef, length(mpn.layers)))
    for (i, l) in enumerate(mpn.layers)
        buf[i] = l(args...)
    end
    return Tuple(copy(buf))
end


text:
# Low-level wrappers for Portaudio calls

# General type aliases
const PaTime = Cdouble
const PaError = Cint
const PaSampleFormat = Culong
const PaDeviceIndex = Cint
const PaHostApiIndex = Cint
const PaHostApiTypeId = Cint
# PaStream is always used as an opaque type, so we're always dealing
# with the pointer
const PaStream = Ptr{Cvoid}
const PaStreamCallback = Cvoid
const PaStreamFlags = Culong

const PA_NO_FLAG = PaStreamFlags(0x00)

@enum PortAudioError(
    paNoError = 0,
    paNotInitialized = -10000,
    paUnanticipatedHostError,
    paInvalidChannelCount,
    paInvalidSampleRate,
    paInvalidDevice,
    paInvalidFlag,
    paSampleFormatNotSupported,
    paBadIODeviceCombination,
    paInsufficientMemory,
    paBufferTooBig,
    paBufferTooSmall,
    paNullCallback,
    paBadStreamPtr,
    paTimedOut,
    paInternalError,
    paDeviceUnavailable,
    paIncompatibleHostApiSpecificStreamInfo,
    paStreamIsStopped,
    paStreamIsNotStopped,
    paInputOverflowed,
    paOutputUnderflowed,
    paHostApiNotFound,
    paInvalidHostApi,
    paCanNotReadFromACallbackStream,
    paCanNotWriteToACallbackStream,
    paCanNotReadFromAnOutputOnlyStream,
    paCanNotWriteToAnInputOnlyStream,
    paIncompatibleStreamHostApi,
    paBadBufferPtr
)
const PA_NO_ERROR = 0
const PA_INPUT_OVERFLOWED = -10000 + 19
const PA_OUTPUT_UNDERFLOWED = -10000 + 20

const TYPE_TO_FORMAT = Dict{Type,PaSampleFormat}(
    Float32 => 1,
    Int32 => 2,
    # Int24   => 4,
    Int16 => 8,
    Int8 => 16,
    UInt8 => 32,
    # NonInterleaved => 2^31
)

const PaStreamCallbackResult = Cint

# Callback return values
@enum CallBackResult paContinue = 0 paComplete = 1 paAbort = 2

# because we're calling Pa_ReadStream and PA_WriteStream from separate threads,
# we put a mutex around libportaudio calls
const PA_MUTEX = ReentrantLock()

function Pa_Initialize()
    handle_status(lock(PA_MUTEX) do
        @ccall libportaudio.Pa_Initialize()::PaError
    end)
end

function Pa_Terminate()
    handle_status(lock(PA_MUTEX) do
        @ccall libportaudio.Pa_Terminate()::PaError
    end)
end

Pa_GetVersion() =
    lock(PA_MUTEX) do
        @ccall libportaudio.Pa_GetVersion()::Cint
    end

function Pa_GetVersionText()
    unsafe_string(lock(PA_MUTEX) do
        @ccall libportaudio.Pa_GetVersionText()::Ptr{Cchar}
    end)
end

# Host API Functions

# A Host API is the top-level of the PortAudio hierarchy. Each host API has a
# unique type ID that tells you which native backend it is (JACK, ALSA, ASIO,
# etc.). On a given system you can identify each backend by its index, which
# will range between 0 and Pa_GetHostApiCount() - 1. You can enumerate through
# all the host APIs on the system by iterating through those values.

# PaHostApiTypeId values
const PA_HOST_API_NAMES = Dict{PaHostApiTypeId,String}(
    0 => "In Development", # use while developing support for a new host API
    1 => "Direct Sound",
    2 => "MME",
    3 => "ASIO",
    4 => "Sound Manager",
    5 => "Core Audio",
    7 => "OSS",
    8 => "ALSA",
    9 => "AL",
    10 => "BeOS",
    11 => "WDMKS",
    12 => "Jack",
    13 => "WASAPI",
    14 => "AudioScience HPI",
)

mutable struct PaHostApiInfo
    struct_version::Cint
    api_type::PaHostApiTypeId
    name::Ptr{Cchar}
    device_count::Cint
    default_input_device::PaDeviceIndex
    default_output_device::PaDeviceIndex
end

Pa_GetHostApiInfo(index) = unsafe_load(
    lock(PA_MUTEX) do
        @ccall libportaudio.Pa_GetHostApiInfo(index::PaHostApiIndex)::Ptr{PaHostApiInfo}
    end,
)

# Device Functions

mutable struct PaDeviceInfo
    struct_version::Cint
    name::Ptr{Cchar}
    host_api::PaHostApiIndex
    max_input_channels::Cint
    max_output_channels::Cint
    default_low_input_latency::PaTime
    default_low_output_latency::PaTime
    default_high_input_latency::PaTime
    default_high_output_latency::PaTime
    default_sample_rate::Cdouble
end

Pa_GetDeviceCount() =
    lock(PA_MUTEX) do
        @ccall libportaudio.Pa_GetDeviceCount()::PaDeviceIndex
    end

Pa_GetDeviceInfo(index) = unsafe_load(
    lock(PA_MUTEX) do
        @ccall libportaudio.Pa_GetDeviceInfo(index::PaDeviceIndex)::Ptr{PaDeviceInfo}
    end,
)

Pa_GetDefaultInputDevice() =
    lock(PA_MUTEX) do
        @ccall libportaudio.Pa_GetDefaultInputDevice()::PaDeviceIndex
    end

Pa_GetDefaultOutputDevice() =
    lock(PA_MUTEX) do
        @ccall libportaudio.Pa_GetDefaultOutputDevice()::PaDeviceIndex
    end

# Stream Functions

mutable struct Pa_StreamParameters
    device::PaDeviceIndex
    channel_count::Cint
    sample_format::PaSampleFormat
    suggested_latency::PaTime
    host_API_specific_stream_info::Ptr{Cvoid}
end

mutable struct PaStreamInfo
    struct_version::Cint
    input_latency::PaTime
    output_latency::PaTime
    the_sample_rate::Cdouble
end

# function Pa_OpenDefaultStream(input_channels, output_channels,
#                               sample_format::PaSampleFormat,
#                               the_sample_rate, frames_per_buffer)
#     stream_pointer = Ref{PaStream}(0)
#     handle_status(@ccall libportaudio.Pa_OpenDefaultStream(
#         stream_pointer::Ref{PaStream},
#         input_channels::Cint,
#         output_channels::Cint,
#         sample_format::PaSampleFormat,
#         the_sample_rate::Cdouble,
#         frames_per_buffer::Culong,
#         C_NULL::Ref{Cvoid},
#         C_NULL::Ref{Cvoid}
#     )::PaError)
#
#     stream_pointer[]
# end
#
function Pa_OpenStream(
    input_parameters,
    output_parameters,
    the_sample_rate,
    frames_per_buffer,
    flags::PaStreamFlags,
    callback,
    userdata,
)
    stream_pointer = Ref{PaStream}(0)
    handle_status(
        lock(PA_MUTEX) do
            @ccall libportaudio.Pa_OpenStream(
                stream_pointer::Ref{PaStream},
                input_parameters::Ref{Pa_StreamParameters},
                output_parameters::Ref{Pa_StreamParameters},
                float(the_sample_rate)::Cdouble,
                frames_per_buffer::Culong,
                flags::PaStreamFlags,
                (callback === nothing ? C_NULL : callback)::Ref{Cvoid},
                (userdata === nothing ? C_NULL : pointer_from_objref(userdata))::Ptr{Cvoid},
            )::PaError
        end,
    )
    stream_pointer[]
end

function Pa_StartStream(stream::PaStream)
    handle_status(lock(PA_MUTEX) do
        @ccall libportaudio.Pa_StartStream(stream::PaStream)::PaError
    end)
end

function Pa_StopStream(stream::PaStream)
    handle_status(lock(PA_MUTEX) do
        @ccall libportaudio.Pa_StopStream(stream::PaStream)::PaError
    end)
end

function Pa_CloseStream(stream::PaStream)
    handle_status(lock(PA_MUTEX) do
        @ccall libportaudio.Pa_CloseStream(stream::PaStream)::PaError
    end)
end

function Pa_GetStreamReadAvailable(stream::PaStream)
    available = lock(PA_MUTEX) do
        @ccall libportaudio.Pa_GetStreamReadAvailable(stream::PaStream)::Clong
    end
    available >= 0 || handle_status(available)
    available
end

function Pa_GetStreamWriteAvailable(stream::PaStream)
    available = lock(PA_MUTEX) do
        @ccall libportaudio.Pa_GetStreamWriteAvailable(stream::PaStream)::Clong
    end
    available >= 0 || handle_status(available)
    available
end

function Pa_ReadStream(stream::PaStream, buf::Array, frames::Integer, show_warnings = true)
    # without disable_sigint I get a segfault with the error:
    # "error thrown and no exception handler available."
    # if the user tries to ctrl-C. Note I've still had some crash problems with
    # ctrl-C within `pasuspend`, so for now I think either don't use `pasuspend` or
    # don't use ctrl-C.
    error_code = disable_sigint() do
        fetch(
            @spawn lock(PA_MUTEX) do
                @ccall libportaudio.Pa_ReadStream(
                    stream::PaStream,
                    buf::Ptr{Cvoid},
                    frames::Culong,
                )::PaError
            end
        )
    end
    handle_status(error_code, show_warnings)
    error_code
end

function Pa_WriteStream(stream::PaStream, buf::Array, frames::Integer, show_warnings = true)
    error_code = disable_sigint() do
        fetch(
            @spawn lock(PA_MUTEX) do
                @ccall libportaudio.Pa_WriteStream(
                    stream::PaStream,
                    buf::Ptr{Cvoid},
                    frames::Culong,
                )::PaError
            end
        )
    end
    handle_status(error_code, show_warnings)
    error_code
end

# function Pa_GetStreamInfo(stream::PaStream)
#     unsafe_load(@ccall libportaudio.Pa_GetStreamInfo(stream::PaStream)::Ptr{PaStreamInfo})
# end
#
# General utility function to handle the status from the Pa_* functions
function handle_status(error_code::PaError, show_warnings::Bool = true)
    if error_code == PA_OUTPUT_UNDERFLOWED || error_code == PA_INPUT_OVERFLOWED
        if show_warnings
            @warn(
                "libportaudio: " * unsafe_string(
                    lock(PA_MUTEX) do
                        @ccall libportaudio.Pa_GetErrorText(
                            error_code::PaError,
                        )::Ptr{Cchar}
                    end,
                )
            )
        end
    elseif error_code != PA_NO_ERROR
        throw(
            ErrorException(
                "libportaudio: " * unsafe_string(
                    lock(PA_MUTEX) do
                        @ccall libportaudio.Pa_GetErrorText(
                            error_code::PaError,
                        )::Ptr{Cchar}
                    end,
                ),
            ),
        )
    end
end


text:
import FractionalTimeDG
using OffsetArrays
using PyPlot

include("../FPDE_problem.jl")

Nx = 8
Nt = 4
rx = 4
rt = 3
Mt = 2rt
ppImax = 10
store = FractionalTimeDG.Store(α, rt, Mt, ppImax)

x = OffsetVector(range(0, L, length=Nx+1), 0:Nx)
grid = FractionalTimeDG.SpatialGrid(x, rx, ("essential", "essential"), rx)

U0 = FractionalTimeDG.spatial_L2_projection(u0, grid)
t, U = FractionalTimeDG.FPDEDG(1.0, f, U0, grid, max_t, Nt, rt, Mt, store)

#pts_per_time_interval = 5
#pts_per_space_interval = 5
#xvals, pcwise_t, pcwise_U = FractionalTimeDG.evaluate_pcwise_poly(
#    U, t, grid, pts_per_time_interval, pts_per_space_interval)
τ = [ -1.0, -0.5, 0.0, 0.5, 1.0 ]
pts_per_time_interval = length(τ)
ξ = [ -0.9, -0.5, 0.0, 0.5, 0.9 ]
xvals, pcwise_t, pcwise_U = FractionalTimeDG.evaluate_pcwise_poly(
    U, t, grid, τ, ξ)
τ = [ -1.0, -0.5, 0.0, 0.5, 1.0 ]

figure(1)
Nxvals = length(xvals)
Ntvals = pts_per_time_interval
X = Float64[ xvals[i] for i = 1:Nxvals, j = 1:Ntvals ]
for n = 1:Nt
    T = Float64[ pcwise_t[j,n] for i = 1:Nxvals, j = 1:Ntvals ]
    Z = Float64[ pcwise_U[i,j,n] for i = 1:Nxvals, j = 1:Ntvals ]
    surf(X, T, Z)
end
xlabel(L"$x$")
ylabel(L"$t$")


text:
using JetPack, Jets, Test, LinearAlgebra

n1,n2,n3 = 10,11,12

@testset "JopLeakyIntegration, 1D, correctness, T=$T" for T in (Float32, Float64, Complex{Float32}, Complex{Float64}), α in (1.0, real(rand(T)))
    n1 = 6
    A = JopLeakyIntegration(JetSpace(T,n1), α=α)
    x = ones(domain(A))
    y1 = A * x
    y2 = zeros(domain(A))
    y2[1] = T(1   + 0   + 0   + 0   + 0   +   0)
    y2[2] = T(α^1 + 1   + 0   + 0   + 0   +   0)
    y2[3] = T(α^2 + α^1 + 1   + 0   + 0   +   0)
    y2[4] = T(α^3 + α^2 + α^1 + 1   + 0   +   0)
    y2[5] = T(α^4 + α^3 + α^2 + α^1 + 1   +   0)
    y2[6] = T(α^5 + α^4 + α^3 + α^2 + α^1 +   1)
    @test norm(y1 .- y2) < 10 * eps(real(T))
end

@testset "JopLeakyIntegration, 1D, adjoint, T=$T" for T in (Float32, Float64, Complex{Float32}, Complex{Float64}), α in (1.0, real(rand(T)))
    A = JopLeakyIntegration(JetSpace(T,n1), α=α)
    lhs, rhs = dot_product_test(A, rand(domain(A)), rand(range(A)))
    @test lhs ≈ rhs
end

@testset "JopLeakyIntegration, 2D, adjoint, T=$T" for T in (Float32, Float64, Complex{Float32}, Complex{Float64}), α in (1.0, real(rand(T)))
    A = JopLeakyIntegration(JetSpace(T,n1,n2), α=α)
    lhs, rhs = dot_product_test(A, rand(domain(A)), rand(range(A)))
    @test lhs ≈ rhs
end

@testset "JopLeakyIntegration, 3D, adjoint, T=$T" for T in (Float32, Float64, Complex{Float32}, Complex{Float64}), α in (1.0, real(rand(T)))
    A = JopLeakyIntegration(JetSpace(T,n1,n2,n3), α=α)
    lhs, rhs = dot_product_test(A, rand(domain(A)), rand(range(A)))
    @test lhs ≈ rhs
end

@testset "JopLeakyIntegration, 1D, linearity, T=$T" for T in (Float32, Float64, Complex{Float32}, Complex{Float64}), α in (1.0, real(rand(T)))
    A = JopLeakyIntegration(JetSpace(T,n1), α=α)
    lhs,rhs = linearity_test(A)
    @test lhs ≈ rhs
    lhs,rhs = linearity_test(A')
    @test lhs ≈ rhs
end

@testset "JopLeakyIntegration, 2D, linearity, T=$T" for T in (Float32, Float64, Complex{Float32}, Complex{Float64}), α in (1.0, real(rand(T)))
    A = JopLeakyIntegration(JetSpace(T,n1,n2), α=α)
    lhs,rhs = linearity_test(A)
    @test lhs ≈ rhs
    lhs,rhs = linearity_test(A')
    @test lhs ≈ rhs
end

@testset "JopLeakyIntegration, 2D, linearity, T=$T" for T in (Float32, Float64, Complex{Float32}, Complex{Float64}), α in (1.0, real(rand(T)))
    A = JopLeakyIntegration(JetSpace(T,n1,n2), α=α)
    lhs,rhs = linearity_test(A)
    @test lhs ≈ rhs
    lhs,rhs = linearity_test(A')
    @test lhs ≈ rhs
end

@testset "JopLeakyIntegration, 3D, linearity, T=$T" for T in (Float32, Float64, Complex{Float32}, Complex{Float64}), α in (1.0, real(rand(T)))
    A = JopLeakyIntegration(JetSpace(T,n1,n2,n3), α=α)
    lhs,rhs = linearity_test(A)
    @test lhs ≈ rhs
    lhs,rhs = linearity_test(A')
    @test lhs ≈ rhs
end


text:
using ConvexHull, JuMP, Polyhedra, Test

@testset "Throws error in case of linearity" begin
    simplex = LiftedHRepresentation([0 1 0; 0 0 1; 1 -1 -1], BitSet([3]))
    @test_throws ErrorException double_description(simplex)
    linv = LiftedHRepresentation([1 1], BitSet([1]))
    @test_throws ErrorException double_description(linv)
end

include("polyhedron.jl")


text:
module TestMOIWrapper

using Gurobi
using Random
using Test

const MOI = Gurobi.MOI

function runtests()
    for name in names(@__MODULE__; all = true)
        if !startswith("$(name)", "test_")
            continue
        end
        @testset "$(name)" begin
            if startswith("$(name)", "test_MULTI_ENV")
                try
                    getfield(@__MODULE__, name)()
                catch ex
                    if ex == ErrorException(
                        "Gurobi Error 10009: Failed to obtain a valid license",
                    )
                        @warn(
                            "Skipping a test because there was an issue " *
                            "creating multiple licenses. This is probably " *
                            "because you have a limited license."
                        )
                    else
                        rethrow(ex)
                    end
                end
            else
                getfield(@__MODULE__, name)()
            end
        end
    end
    return
end

const GRB_ENV =
    isdefined(Main, :GRB_ENV) ? Main.GRB_ENV : Gurobi.Env(output_flag = 0)

function test_runtests()
    model =
        MOI.Bridges.full_bridge_optimizer(Gurobi.Optimizer(GRB_ENV), Float64)
    MOI.set(model, MOI.Silent(), true)
    # We set `DualReductions = 0` so that we never return
    # `INFEASIBLE_OR_UNBOUNDED`.
    MOI.set(model, MOI.RawOptimizerAttribute("DualReductions"), 0)
    MOI.set(model, MOI.RawOptimizerAttribute("QCPDual"), 1)
    MOI.set(model, MOI.RawOptimizerAttribute("InfUnbdInfo"), 1)
    MOI.set(model, MOI.RawOptimizerAttribute("NonConvex"), 2)
    MOI.Test.runtests(
        model,
        MOI.Test.Config(atol = 1e-3, rtol = 1e-3),
        exclude = String[
            # TODO(odow): new tests from MOI that need implementing
            "test_infeasible_",
            "test_unbounded_",
            # TODO(odow): investigate errors
            # Gurobi Error 10015: Cannot compute IIS on a feasible model
            # https://www.gurobi.com/documentation/9.5/refman/error_codes.html
            "test_solve_conflict_feasible",
            # SecondOrderCone does not return dual solutions. Tested below.
            "_SecondOrderCone_",
            "test_constraint_PrimalStart_DualStart_SecondOrderCone",
            "_RotatedSecondOrderCone_",
            "_GeometricMeanCone_",
        ],
    )
    MOI.Test.runtests(
        model,
        MOI.Test.Config(
            # Some of these conic ones have very low accuracy.
            atol = 2e-3,
            rtol = 1e-3,
            exclude = Any[MOI.ConstraintDual, MOI.DualObjectiveValue],
        ),
        include = String[
            "_SecondOrderCone_",
            "test_constraint_PrimalStart_DualStart_SecondOrderCone",
            "_RotatedSecondOrderCone_",
            "_GeometricMeanCone_",
        ],
    )
    return
end

function test_User_limit_handling_issue_140()
    # Verify that we return the correct status codes when a mixed-integer
    # problem has been solved to a *feasible* but not necessarily optimal
    # solution. To do that, we will set up an intentionally dumbed-down
    # Gurobi Gurobi.Optimizer (with all heuristics and pre-solve turned off) and
    # ask it to solve a classic knapsack problem. Setting SolutionLimit=1
    # forces the solver to return after its first feasible MIP solution,
    # which tests the right part of the code without relying on potentially
    # flaky or system-dependent time limits.
    m = Gurobi.Optimizer(GRB_ENV)
    MOI.set(m, MOI.Silent(), true)
    MOI.set(m, MOI.RawOptimizerAttribute("SolutionLimit"), 1)
    MOI.set(m, MOI.RawOptimizerAttribute("Heuristics"), 0)
    MOI.set(m, MOI.RawOptimizerAttribute("Presolve"), 0)
    N = 100
    x = MOI.add_variables(m, N)
    for xi in x
        MOI.add_constraint(m, xi, MOI.ZeroOne())
        MOI.set(m, MOI.VariablePrimalStart(), xi, 0.0)
    end
    # Given a collection of items with individual weights and values,
    # maximize the total value carried subject to the constraint that
    # the total weight carried is less than 10.
    Random.seed!(1)
    item_weights = rand(N)
    item_values = rand(N)
    MOI.add_constraint(
        m,
        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(item_weights, x), 0.0),
        MOI.LessThan(10.0),
    )
    MOI.set(
        m,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(.-item_values, x), 0.0),
    )
    MOI.optimize!(m)

    @test MOI.get(m, MOI.TerminationStatus()) == MOI.SOLUTION_LIMIT
    # We should have a primal feasible solution:
    @test MOI.get(m, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT
    # But we have no dual status:
    @test MOI.get(m, MOI.DualStatus()) == MOI.NO_SOLUTION
    return
end

function test_Constant_objective_issue_111()
    m = Gurobi.Optimizer(GRB_ENV)
    x = MOI.add_variable(m)
    MOI.set(
        m,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm{Float64}[], 2.0),
    )
    @test MOI.get(
        m,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
    ).constant == 2.0
    p = Ref{Cdouble}()
    Gurobi.GRBgetdblattr(m, "ObjCon", p)
    @test p[] == 2.0

    MOI.modify(
        m,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
        MOI.ScalarConstantChange(3.0),
    )
    @test MOI.get(
        m,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
    ).constant == 3.0
    Gurobi.GRBgetdblattr(m.inner, "ObjCon", p)
    @test p[] == 3.0
    return
end

function test_user_provided_env()
    model_1 = Gurobi.Optimizer(GRB_ENV)
    @test model_1.env === GRB_ENV
    model_2 = Gurobi.Optimizer(GRB_ENV)
    @test model_2.env === GRB_ENV
    # Check that finalizer doesn't touch GRB_ENV when manually provided.
    finalize(model_1)
    @test GRB_ENV.ptr_env != C_NULL
    return
end

function test_MULTI_ENV()
    # Gurobi tests should pass if the function begins with test_MULTI_ENV and
    # this specific error is thrown.
    return error("Gurobi Error 10009: Failed to obtain a valid license")
end

function test_MULTI_ENV_automatic_env()
    model_1 = Gurobi.Optimizer()
    model_2 = Gurobi.Optimizer()
    @test model_1.env !== model_2.env
    # Check that env is finalized with model when not supplied manually.
    finalize(model_1)
    @test model_1.env.ptr_env == C_NULL
    return
end

function test_user_provided_env_empty()
    model = Gurobi.Optimizer(GRB_ENV)
    @test model.env === GRB_ENV
    @test GRB_ENV.ptr_env != C_NULL
    MOI.empty!(model)
    @test model.env === GRB_ENV
    @test GRB_ENV.ptr_env != C_NULL
    return
end

function test_MULTI_ENV_automatic_env_empty()
    model = Gurobi.Optimizer()
    env = model.env
    MOI.empty!(model)
    @test model.env === env
    @test env.ptr_env != C_NULL
    return
end

function test_MULTI_ENV_manual_finalize()
    env = Gurobi.Env(output_flag = 0)
    model = Gurobi.Optimizer(env)
    finalize(env)
    @test env.finalize_called
    finalize(model)
    @test env.ptr_env == C_NULL
    return
end

function test_QCPDual_1()
    model = Gurobi.Optimizer(GRB_ENV)
    MOI.set(model, MOI.Silent(), true)
    MOI.set(model, MOI.RawOptimizerAttribute("QCPDual"), 1)
    MOI.Utilities.loadfromstring!(
        model,
        """
variables: x, y, z
minobjective: 1.0 * x + 1.0 * y + 1.0 * z
c1: x + y == 2.0
c2: x + y + z >= 0.0
c3: 1.0 * x * x + -1.0 * y * y + -1.0 * z * z >= 0.0
x >= 0.0
y >= 0.0
z >= 0.0
""",
    )
    MOI.optimize!(model)
    @test MOI.get(model, MOI.TerminationStatus()) == MOI.OPTIMAL
    @test MOI.get(model, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT
    @test MOI.get(model, MOI.DualStatus()) == MOI.FEASIBLE_POINT
    c1 = MOI.get(model, MOI.ConstraintIndex, "c1")
    c2 = MOI.get(model, MOI.ConstraintIndex, "c2")
    c3 = MOI.get(model, MOI.ConstraintIndex, "c3")
    @test MOI.get(model, MOI.ConstraintDual(), c1) ≈ 1.0 atol = 1e-6
    @test MOI.get(model, MOI.ConstraintDual(), c2) ≈ 0.0 atol = 1e-6
    @test MOI.get(model, MOI.ConstraintDual(), c3) ≈ 0.0 atol = 1e-6
    return
end

function test_QCPDual_default()
    model = Gurobi.Optimizer(GRB_ENV)
    MOI.set(model, MOI.Silent(), true)
    MOI.Utilities.loadfromstring!(
        model,
        """
variables: x, y, z
minobjective: 1.0 * x + 1.0 * y + 1.0 * z
c1: x + y == 2.0
c2: x + y + z >= 0.0
c3: 1.0 * x * x + -1.0 * y * y + -1.0 * z * z >= 0.0
x >= 0.0
y >= 0.0
z >= 0.0
""",
    )
    MOI.optimize!(model)
    @test MOI.get(model, MOI.TerminationStatus()) == MOI.OPTIMAL
    @test MOI.get(model, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT
    @test MOI.get(model, MOI.DualStatus()) == MOI.NO_SOLUTION
    c1 = MOI.get(model, MOI.ConstraintIndex, "c1")
    c2 = MOI.get(model, MOI.ConstraintIndex, "c2")
    c3 = MOI.get(model, MOI.ConstraintIndex, "c3")
    @test_throws ErrorException MOI.get(model, MOI.ConstraintDual(), c1)
    @test_throws ErrorException MOI.get(model, MOI.ConstraintDual(), c2)
    @test_throws ErrorException MOI.get(model, MOI.ConstraintDual(), c3)
    return
end

function test_Buffered_deletion_test()
    # Check if the VarHintVal of variables (could be any attribute that is not
    # cached, the most common as: the bounds, name, and the start, are all
    # cached) are right if we remove some variables and add new ones without
    # updating the model (and that they stay correct after updating the model,
    # here done by `MOI.optimize!`, as the update needs to be done inside MOI
    # code, that is aware of the lazy updates).
    model = Gurobi.Optimizer(GRB_ENV)
    MOI.set(model, MOI.Silent(), true)
    vars = MOI.add_variables(model, 3)
    vars_obj = [7.0, 42.0, -0.5]
    obj_attr = Gurobi.VariableAttribute("Obj")
    MOI.set.(model, obj_attr, vars, vars_obj)
    @test all(MOI.is_valid.(model, vars))
    @test MOI.get.(model, obj_attr, vars) == vars_obj
    MOI.delete(model, vars[[1, 3]])
    @test !MOI.is_valid(model, vars[1])
    @test !MOI.is_valid(model, vars[3])
    fourth_var = MOI.add_variable(model)
    fourth_var_obj = -77.0
    MOI.set(model, obj_attr, fourth_var, fourth_var_obj)
    # Check before updating the model and after the delete+insert.
    @test !MOI.is_valid(model, vars[1])
    @test !MOI.is_valid(model, vars[3])
    @test MOI.is_valid(model, vars[2])
    @test vars_obj[2] == MOI.get(model, obj_attr, vars[2])
    @test MOI.is_valid(model, fourth_var)
    @test fourth_var_obj == MOI.get(model, obj_attr, fourth_var)
    # Then optimize to force an update.
    MOI.optimize!(model)
    # And check again.
    @test MOI.is_valid(model, vars[2])
    @test vars_obj[2] == MOI.get(model, obj_attr, vars[2])
    @test MOI.is_valid(model, fourth_var)
    @test fourth_var_obj == MOI.get(model, obj_attr, fourth_var)
    return
end

function test_ConstraintAttribute()
    model = Gurobi.Optimizer(GRB_ENV)
    MOI.Utilities.loadfromstring!(
        model,
        """
variables: x
minobjective: x
x >= 0.0
c: 2x >= 1.0
x in Integer()
""",
    )
    c2 = MOI.get(model, MOI.ConstraintIndex, "c")
    # Linear constraints are supported - one test for each different type.
    # Integer attribute
    MOI.set(model, Gurobi.ConstraintAttribute("Lazy"), c2, 2)
    @test MOI.get(model, Gurobi.ConstraintAttribute("Lazy"), c2) == 2
    # Real attribute
    MOI.set(model, Gurobi.ConstraintAttribute("RHS"), c2, 2.0)
    @test MOI.get(model, Gurobi.ConstraintAttribute("RHS"), c2) == 2.0
    # Char attribute
    MOI.set(model, Gurobi.ConstraintAttribute("Sense"), c2, '<')
    @test MOI.get(model, Gurobi.ConstraintAttribute("Sense"), c2) == '<'
    # String attribute
    MOI.set(model, Gurobi.ConstraintAttribute("ConstrName"), c2, "c4")
    @test MOI.get(model, Gurobi.ConstraintAttribute("ConstrName"), c2) == "c4"
    # Getting/setting a non-existing attribute.
    attr = Gurobi.ConstraintAttribute("Non-existing")
    err = ErrorException("Gurobi Error 10004: Unknown attribute 'Non-existing'")
    @test_throws err MOI.set(model, attr, c2, 1)
    @test_throws err MOI.get(model, attr, c2)
    # Setting an attribute to a value of the wrong type.
    @test_throws(
        ArgumentError(
            "Attribute Lazy requires Int64 arguments. Provided argument was of type Float64.",
        ),
        MOI.set(model, Gurobi.ConstraintAttribute("Lazy"), c2, 1.5)
    )
    return
end

function test_VariableAttribute()
    model = Gurobi.Optimizer(GRB_ENV)
    MOI.Utilities.loadfromstring!(
        model,
        """
variables: x
minobjective: x
c1: x >= 0.0
c2: 2x >= 1.0
c3: x in Integer()
""",
    )
    x = MOI.get(model, MOI.VariableIndex, "x")
    # Setting attributes of each type
    # Integer attribute
    MOI.set(model, Gurobi.VariableAttribute("VarHintPri"), x, 2)
    @test MOI.get(model, Gurobi.VariableAttribute("VarHintPri"), x) == 2
    # Real Attribute
    MOI.set(model, Gurobi.VariableAttribute("LB"), x, 2.0)
    @test MOI.get(model, Gurobi.VariableAttribute("LB"), x) == 2.0
    # Char Attribute
    MOI.set(model, Gurobi.VariableAttribute("VType"), x, 'B')
    @test MOI.get(model, Gurobi.VariableAttribute("VType"), x) == 'B'
    # String Attribute
    MOI.set(model, Gurobi.VariableAttribute("VarName"), x, "my_var")
    @test MOI.get(model, Gurobi.VariableAttribute("VarName"), x) == "my_var"
    # Things that should fail follow.
    # Getting/setting a non-existing attribute.
    attr = Gurobi.VariableAttribute("Non-existing")
    err = ErrorException("Gurobi Error 10004: Unknown attribute 'Non-existing'")
    @test_throws err MOI.set(model, attr, x, 1)
    @test_throws err MOI.get(model, attr, x)
    # Setting an attribute to a value of the wrong type.
    @test_throws(
        ArgumentError(
            "Attribute BranchPriority requires Int64 arguments. Provided argument was of type Float64.",
        ),
        MOI.set(model, Gurobi.VariableAttribute("BranchPriority"), x, 1.5)
    )
    return
end

function test_ModelAttribute()
    model = Gurobi.Optimizer(GRB_ENV)
    MOI.Utilities.loadfromstring!(
        model,
        """
variables: x
minobjective: x
c1: x >= 0.0
c2: 2x >= 1.0
c3: x in Integer()
""",
    )
    # Setting attributes of each type
    # Integer attribute
    MOI.set(model, Gurobi.ModelAttribute("ModelSense"), -1)
    @test MOI.get(model, Gurobi.ModelAttribute("ModelSense")) == -1
    # Real Attribute
    MOI.set(model, Gurobi.ModelAttribute("ObjCon"), 3.0)
    @test MOI.get(model, Gurobi.ModelAttribute("ObjCon")) == 3.0
    # String Attribute
    MOI.set(model, Gurobi.ModelAttribute("ModelName"), "My model")
    @test MOI.get(model, Gurobi.ModelAttribute("ModelName")) == "My model"
    # Things that should fail follow.
    # Getting/setting a non-existing attribute.
    attr = Gurobi.ModelAttribute("Non-existing")
    err = ErrorException("Gurobi Error 10004: Unknown attribute 'Non-existing'")
    @test_throws err MOI.set(model, attr, 1)
    @test_throws err MOI.get(model, attr)
    # Setting an attribute to a value of the wrong type.
    @test_throws(
        ArgumentError(
            "Attribute NumStart requires $(Int) arguments. Provided argument was of type Float64.",
        ),
        MOI.set(model, Gurobi.ModelAttribute("NumStart"), 4.0)
    )
    return
end

function test_Attributes()
    model = Gurobi.Optimizer(GRB_ENV)
    MOI.set(model, MOI.Silent(), true)
    MOI.optimize!(model)
    @test MOI.get(model, MOI.NodeCount()) == 0
    @test MOI.get(model, MOI.BarrierIterations()) == 0
    @test MOI.get(model, MOI.SimplexIterations()) == 0
end

function test_GRBterminate()
    model = Gurobi.Optimizer(GRB_ENV)
    MOI.set(model, MOI.Silent(), true)
    x = MOI.add_variable(model)
    MOI.set(
        model,
        Gurobi.CallbackFunction(),
        (cb_data, cb_where) -> GRBterminate(model),
    )
    MOI.optimize!(model)
    @test MOI.get(model, MOI.TerminationStatus()) == MOI.INTERRUPTED
    return
end

"""
    test_InterruptException()

This test simulates an InterruptException being thrown. It is a little
complicated due to the delayed handling of GRBterminate, which _schedules_ a
request for termination, rather than terminating immediately. This means Gurobi
may continue to call the callback after the interruption.

First, we must ensure that InterruptException() is only thrown once. Double
interrupting would interrupt our handling of the first interrupt!

Second, if the model is too simplisitic, Gurobi may be able to prove optimality
after we have interrupted, but before it has decided to actually exit the solve.
"""
function test_InterruptException()
    model = Gurobi.Optimizer(GRB_ENV)
    MOI.set(model, MOI.Silent(), true)
    x = MOI.add_variable(model)
    MOI.add_constraint(model, x, MOI.Integer())
    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)
    MOI.set(model, MOI.ObjectiveFunction{MOI.VariableIndex}(), x)
    MOI.set(model, MOI.RawOptimizerAttribute("LazyConstraints"), 1)
    interrupt_thrown = false
    i = 0.0
    MOI.set(
        model,
        Gurobi.CallbackFunction(),
        (cb_data, cb_where) -> begin
            if cb_where != Gurobi.GRB_CB_MIPSOL
                return
            end
            MOI.submit(
                model,
                MOI.LazyConstraint(cb_data),
                MOI.ScalarAffineFunction{Float64}(
                    [MOI.ScalarAffineTerm(1.0, x)],
                    0.0,
                ),
                MOI.GreaterThan{Float64}(i),
            )
            i += 1
            if !interrupt_thrown
                interrupt_thrown = true
                throw(InterruptException())
            end
        end,
    )
    MOI.optimize!(model)
    @test MOI.get(model, MOI.TerminationStatus()) == MOI.INTERRUPTED
    return
end

function _build_basis_model()
    T = Float64
    model = Gurobi.Optimizer(GRB_ENV)
    MOI.set(model, MOI.Silent(), true)
    # Min -x
    # s.t. x + y <= 1
    # x, y >= 0

    x = MOI.add_variable(model)
    y = MOI.add_variable(model)

    cf = MOI.ScalarAffineFunction{T}(
        MOI.ScalarAffineTerm{T}.([one(T), one(T)], [x, y]),
        zero(T),
    )
    c = MOI.add_constraint(model, cf, MOI.LessThan(one(T)))

    vc1 = MOI.add_constraint(model, x, MOI.GreaterThan(zero(T)))
    vc2 = MOI.add_constraint(model, y, MOI.GreaterThan(zero(T)))
    objf = MOI.ScalarAffineFunction{T}(
        MOI.ScalarAffineTerm{T}.([-one(T), zero(T)], [x, y]),
        zero(T),
    )
    MOI.set(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{T}}(), objf)
    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)
    return model, x, y, c
end

function test_set_basis()
    # The following is an indirect, brittle way to test the desired behavior.
    # Gurobi appears to not allow you to query VBasis/CBasis attributes without
    # optimizing the problem
    # (https://support.gurobi.com/hc/en-us/community/posts/360075729911-Getting-VBasis-CBasis-attributes-without-optimizing).
    # The following just verifies that the problem is solved with 0 simplex iterations
    # when we seed the optimal solution and 1 simplex iteration when we seed a
    # suboptimal feasible basis. Ideally, we would remove the call to MOI.optimize!
    # and just verify that, after setting the basis, we can get the same basis
    # statuses back.
    let
        model, x, y, c = _build_basis_model()
        MOI.set(model, MOI.VariableBasisStatus(), x, MOI.BASIC)
        MOI.set(model, MOI.VariableBasisStatus(), y, MOI.NONBASIC_AT_LOWER)
        MOI.set(model, MOI.ConstraintBasisStatus(), c, MOI.NONBASIC)
        MOI.optimize!(model)
        @test MOI.get(model, MOI.SimplexIterations()) == 0
    end
    let
        model, x, y, c = _build_basis_model()
        MOI.set(model, MOI.VariableBasisStatus(), x, MOI.NONBASIC_AT_LOWER)
        MOI.set(model, MOI.VariableBasisStatus(), y, MOI.NONBASIC_AT_LOWER)
        MOI.set(model, MOI.ConstraintBasisStatus(), c, MOI.BASIC)
        MOI.optimize!(model)
        @test MOI.get(model, MOI.SimplexIterations()) == 1
    end
    return
end

function test_add_constrained_variables()
    model = Gurobi.Optimizer(GRB_ENV)
    MOI.set(model, MOI.Silent(), true)
    set = MOI.Interval{Float64}(-1.2, 3.4)
    vi, ci = MOI.add_constrained_variable(model, set)
    @test MOI.get(model, MOI.NumberOfVariables()) == 1
    @test MOI.get(model, MOI.ListOfConstraintTypesPresent()) ==
          [(MOI.VariableIndex, MOI.Interval{Float64})]
    @test MOI.get(model, MOI.ConstraintFunction(), ci) == vi
    @test MOI.get(model, MOI.ConstraintSet(), ci) == set
    return
end

function _is_binary(x; atol = 1e-6)
    return isapprox(x, 0; atol = atol) || isapprox(x, 1; atol = atol)
end

function test_multiple_solutions()
    model = Gurobi.Optimizer(GRB_ENV)
    MOI.set(model, MOI.RawOptimizerAttribute("OutputFlag"), 0)
    MOI.set(model, MOI.RawOptimizerAttribute("Cuts"), 0)
    MOI.set(model, MOI.RawOptimizerAttribute("Presolve"), 0)
    MOI.set(model, MOI.RawOptimizerAttribute("PreCrush"), 0)
    MOI.set(model, MOI.RawOptimizerAttribute("Heuristics"), 0)
    N = 30
    x = MOI.add_variables(model, N)
    MOI.add_constraints(model, x, MOI.ZeroOne())
    MOI.set.(model, MOI.VariablePrimalStart(), x, 0.0)
    Random.seed!(1)
    item_weights, item_values = rand(N), rand(N)
    MOI.add_constraint(
        model,
        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(item_weights, x), 0.0),
        MOI.LessThan(10.0),
    )
    MOI.set(
        model,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(item_values, x), 0.0),
    )
    MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE)
    MOI.optimize!(model)
    RC = MOI.get(model, MOI.ResultCount())
    @test RC > 1
    for n in [0, RC + 1]
        @test_throws(
            MOI.ResultIndexBoundsError,
            MOI.get(model, MOI.VariablePrimal(n), x),
        )
        @test_throws(
            MOI.ResultIndexBoundsError,
            MOI.get(model, MOI.ObjectiveValue(n)),
        )
        @test MOI.get(model, MOI.PrimalStatus(n)) == MOI.NO_SOLUTION
        @test MOI.get(model, MOI.DualStatus(n)) == MOI.NO_SOLUTION
    end
    for n in 1:RC
        xn = MOI.get(model, MOI.VariablePrimal(n), x)
        @test all(_is_binary, xn)
        @test isapprox(
            MOI.get(model, MOI.ObjectiveValue(n)),
            item_values' * xn,
            atol = 1e-6,
        )
        @test MOI.get(model, MOI.PrimalStatus(n)) == MOI.FEASIBLE_POINT
        @test MOI.get(model, MOI.DualStatus(n)) == MOI.NO_SOLUTION
    end
    return
end

function test_log_file()
    model = Gurobi.Optimizer()
    @test MOI.get(model, MOI.RawOptimizerAttribute("LogFile")) == ""
    MOI.set(model, MOI.RawOptimizerAttribute("LogFile"), "test.log")
    @test MOI.get(model, MOI.RawOptimizerAttribute("LogFile")) == "test.log"
    MOI.set(model, MOI.RawOptimizerAttribute("LogFile"), "")
    rm("test.log")
    return
end

end

TestMOIWrapper.runtests()


text:
# Infrastructure for value transformation

"""
    clamp01(x) -> y

Produce a value `y` that lies between 0 and 1, and equal to `x` when
`x` is already in this range. Equivalent to `clamp(x, 0, 1)` for
numeric values. For colors, this function is applied to each color
channel separately.

See also: [`clamp01!`](@ref), [`clamp01nan`](@ref).
"""
clamp01(x::Union{N0f8,N0f16}) = x
clamp01(x::Number) = clamp(x, zero(x), one(x))
clamp01(c::Colorant) = mapc(clamp01, c)

"""
    clamp01!(array::AbstractArray)

Restrict values in array to [0, 1], in-place. See also [`clamp01`](@ref).
"""
function clamp01!(img::AbstractArray)
    # slgihtly faster than map!(clamp01, img, img)
    @inbounds for i in eachindex(img)
        img[i] = clamp01(img[i])
    end
    img
end

"""
    clamp01nan(x) -> y

Similar to `clamp01`, except that any `NaN` values are changed to 0.

See also: [`clamp01nan!`](@ref), [`clamp01`](@ref).
"""
clamp01nan(x) = clamp01(x)
clamp01nan(x::AbstractFloat) = ifelse(isnan(x), zero(x), clamp01(x))
clamp01nan(c::Colorant) = mapc(clamp01nan, c)

"""
    clamp01nan!(array::AbstractArray)

Similar to `clamp01!`, except that any `NaN` values are changed to 0.

See also: [`clamp01!`](@ref), [`clamp01nan`](@ref)
"""
function clamp01nan!(img::GenericImage)
    # slgihtly faster than map!(clamp01nan, img, img)
    @inbounds for i in eachindex(img)
        img[i] = clamp01nan(img[i])
    end
    img
end

"""
    scaleminmax(min, max) -> f
    scaleminmax(T, min, max) -> f

Return a function `f` which maps values less than or equal to `min` to
0, values greater than or equal to `max` to 1, and uses a linear scale
in between. `min` and `max` should be real values.

Optionally specify the return type `T`. If `T` is a colorant (e.g.,
RGB), then scaling is applied to each color channel.

# Examples

## Example 1

```julia
julia> f = scaleminmax(-10, 10)
(::#9) (generic function with 1 method)

julia> f(10)
1.0

julia> f(-10)
0.0

julia> f(5)
0.75
```

## Example 2

```julia
julia> c = RGB(255.0,128.0,0.0)
RGB{Float64}(255.0,128.0,0.0)

julia> f = scaleminmax(RGB, 0, 255)
(::#13) (generic function with 1 method)

julia> f(c)
RGB{Float64}(1.0,0.5019607843137255,0.0)
```

See also: [`takemap`](@ref).
"""
function scaleminmax(min::T, max::T) where T
    @inline function(x)
        xp, minp, maxp = promote(x, min, max)  # improves performance to promote explicitly
        y = clamp(xp, minp, maxp)
        (y-minp)/(maxp-minp)
    end
end
function scaleminmax(::Type{Tout}, min::T, max::T) where {Tout,T}
    @inline function(x)
        xp, minp, maxp = promote(x, min, max)
        y = clamp(xp, minp, maxp)
        smmconvert(Tout, (y-minp)/(maxp-minp))
    end
end
@inline smmconvert(::Type{Tout}, x) where {Tout} = convert(Tout, x)
# since we know the result will be between 0 and 1, we can use rem to save a check
@inline smmconvert(::Type{Tout}, x) where {Tout<:Normed} = rem(x, Tout)

scaleminmax(min, max) = scaleminmax(promote(min, max)...)
scaleminmax(::Type{T}, min, max) where {T} = scaleminmax(T, promote(min, max)...)

# TODO: use triangular dispatch when we can count on Julia 0.6+
function scaleminmax(::Type{C}, min::T, max::T) where {C<:Colorant, T<:Real}
    _scaleminmax(C, eltype(C), min, max)
end
function _scaleminmax(::Type{C}, ::Type{TC}, min::T, max::T) where {C<:Colorant, TC<:Real, T<:Real}
    freal = scaleminmax(TC, min, max)
    @inline function(c)
        C(mapc(freal, c))
    end
end
function _scaleminmax(::Type{C}, ::Type{Any}, min::T, max::T) where {C<:Colorant, T<:Real}
    freal = scaleminmax(min, max)
    @inline function(c)
        C(mapc(freal, c))
    end
end

function takemap(::typeof(scaleminmax), A::AbstractArray{T}) where T<:Real
    min, max = extrema(A)
    scaleminmax(min, max)
end
function takemap(::typeof(scaleminmax), A::AbstractArray{C}) where C<:Colorant
    min, max = extrema(channelview(A))
    scaleminmax(C, min, max)
end
function takemap(::typeof(scaleminmax), ::Type{Tout}, A::AbstractArray{T}) where {Tout,T<:Real}
    min, max = extrema(A)
    scaleminmax(Tout, min, max)
end
function takemap(::typeof(scaleminmax), ::Type{Cout}, A::AbstractArray{C}) where {Cout<:Colorant,C<:Colorant}
    min, max = extrema(channelview(A))
    scaleminmax(Cout, min, max)
end

"""
    scalesigned(maxabs) -> f

Return a function `f` which scales values in the range `[-maxabs,
maxabs]` (clamping values that lie outside this range) to the range
`[-1, 1]`.

See also: [`colorsigned`](@ref).
"""
function scalesigned(maxabs::Real)
    maxabs > 0 || throw(ArgumentError("maxabs must be positive, got $maxabs"))
    x -> clamp(x, -maxabs, maxabs)/maxabs
end

"""
    scalesigned(min, center, max) -> f

Return a function `f` which scales values in the range `[min, center]`
to `[-1,0]` and `[center,max]` to `[0,1]`. Values smaller than
`min`/`max` get clamped to `min`/`max`, respectively.

See also: [`colorsigned`](@ref).
"""
function scalesigned(min::T, center::T, max::T) where T<:Real
    min <= center <= max || throw(ArgumentError("values must be ordered, got $min, $center, $max"))
    sneg, spos = 1/(center-min), 1/(max-center)
    function(x)
        Δy = clamp(x, min, max) - center
        ifelse(Δy < 0, sneg*Δy, spos*Δy)
    end
end
scalesigned(min::Real, center::Real, max::Real) = scalesigned(promote(min, center, max)...)

function takemap(::typeof(scalesigned), A::AbstractArray{T}) where T<:Real
    mn, mx = extrema(A)
    scalesigned(max(abs(mn), abs(mx)))
end

"""
    colorsigned()
    colorsigned(colorneg, colorpos) -> f
    colorsigned(colorneg, colorcenter, colorpos) -> f

Define a function that maps negative values (in the range [-1,0]) to
the linear colormap between `colorneg` and `colorcenter`, and positive
values (in the range [0,1]) to the linear colormap between
`colorcenter` and `colorpos`.

The default colors are:

- `colorcenter`: white
- `colorneg`: green1
- `colorpos`: magenta

See also: [`scalesigned`](@ref).
"""
colorsigned(neg::C, center::C, pos::C) where {C<:Colorant} = function(x)
    y = clamp(x, -one(x), one(x))
    yabs = abs(y)
    C(ifelse(y>0, weighted_color_mean(yabs, pos, center),
                  weighted_color_mean(yabs, neg, center)))
end

function colorsigned(colorneg::C, colorpos::C) where C<:Colorant
    colorsigned(colorneg, C(colorant"white"), colorpos)
end

colorsigned() = colorsigned(colorant"green1", colorant"magenta")

"""
    takemap(f, A) -> fnew
    takemap(f, T, A) -> fnew

Given a value-mapping function `f` and an array `A`, return a
"concrete" mapping function `fnew`. When applied to elements of `A`,
`fnew` should return valid values for storage or display, for example
in the range from 0 to 1 (for grayscale) or valid colorants. `fnew`
may be adapted to the actual values present in `A`, and may not
produce valid values for any inputs not in `A`.

Optionally one can specify the output type `T` that `fnew` should produce.

# Example:
```julia
julia> A = [0, 1, 1000];

julia> f = takemap(scaleminmax, A)
(::#7) (generic function with 1 method)

julia> f.(A)
3-element Array{Float64,1}:
 0.0
 0.001
 1.0
```
"""
takemap

takemap(f, A) = f
takemap(f, ::Type{T}, A) where {T} = x->T(f(x))


text:
using OpenQuantumTools, OrdinaryDiffEq, Plots, LaTeXStrings
using Optim, Printf
# schedules.jl defines the annealing schedules
include("../../modules/3_qubit_entanglement/schedules.jl")
# hamiltonians.jl contains functions to construct the Hamiltonians
include("../../modules/3_qubit_entanglement/hamiltonians.jl")
# couplings.jl contains functions to construct the system-bath coupling operators
include("../../modules/3_qubit_entanglement/couplings.jl")
# evolutions.jl contains functions to construct the `Annealing` object
include("../../modules/3_qubit_entanglement/evolutions.jl")

s_axis = range(0, 1, length=200)
plot(
    s_axis,
    [A.(s_axis) B.(s_axis)] / 2 / π,
    linewidth=2,
    label=[L"a(s)" L"b(s)"],
    ylabel=L"\mathrm{Schedule}\ (\mathrm{GHz})",
    xlabel=L"s",
    legend=(0.5, 0.9);
)

# hp value for the instance
hp = 1.2
# construct the Hamiltonian
# the second argument is the stage number
H = build_H(hp, 1)
# get the energy level of the initial state
ilvl = init_lvl(hp)
s_list = range(0, 1, length=201)
ev = []
for s in s_list
    w, = eigen_decomp(H, s, lvl=4)
    push!(ev, w)
end
gap_lvl = ilvl == 1 ? 1 : 2
sm, vm = min_gap(H, gap_lvl, gap_lvl+1)
# create a plot
p1 = plot(s_list, hcat(ev...)', label=[L"E_0" L"E_1" L"E_2" L"E_3"], linewidth=2, xlabel=L"s", ylabel=L"E_i\ (\mathrm{GHz})",
    legend=:bottomright, title=latexstring(@sprintf("\\mathrm{Stage\\ 1}; h_p=%.2f", hp)))
# add a tracking line
sm <= 0.5 && !(hp == 1.8) ? vline!(p1, [sm], label="crossing", linewidth=2) : nothing
plot!(p1, [0.05, 0.0], [ev[1][ilvl], ev[1][ilvl]], arrow=0.4, color=:red, label="initial level")

dH = build_dH(hp, 1)
coupling = ConstantCouplings(["ZII", "IZI", "IIZ"], unit=:ħ)
s_axis = range(0, 1, length=1000)

proj = project_to_lowlevel(H, s_axis, coupling, dH, lvl=5, direction=:backward)
p1 = plot(s_axis, abs.(get_dθ(proj, 1, 2)), legend=:topleft, linewidth=2, label=L"|d\theta_{12}|", ylims=[0,65], ylabel=L"|d\theta|", xlabel=L"s")
plot!(p1, s_axis, abs.(get_dθ(proj, 1, 3)), linewidth=2, label=L"|d\theta_{13}|")
plot!(p1, s_axis, abs.(get_dθ(proj, 2, 3)), linewidth=2, label=L"|d\theta_{23}|")
plot!(p1, s_axis, abs.(get_dθ(proj, 2, 4)), linewidth=2, label=L"|d\theta_{24}|")
plot!(p1, s_axis, abs.(get_dθ(proj, 3, 4)), linewidth=2, label=L"|d\theta_{34}|")
plot!(p1, s_axis, abs.(get_dθ(proj, 1, 4)), linewidth=2, label=L"|d\theta_{14}|")

title!(p1, latexstring(@sprintf "h_p = %.3f" hp))

# define Ohmic bath
bath = Ohmic(1.2732e-4, 4, 12.5)
# construct the Hamiltonian H, system-bath couplings C and
# instantaneous pulse control callback cb for each stage of
# the experiment
# solve stage 1
H, C, cb = build_adiabatic_frame_evolution(hp, 1)
C = C[1]
u1 = build_adiabatic_frame_u0(hp)
annealing = Annealing(H, u1, coupling=C, bath=bath)
sol1 = solve_ame(annealing, 2 * τ1, ω_hint=range(-15, 15, length=200), alg=Tsit5(), reltol=1e-6, callback=cb, saveat=range(0,2*τ1,length=100))

# solve stage 2
# with τ2 given by
τ2 = 1e5
u2 = sol1[end]
H, C, cb = build_adiabatic_frame_evolution(hp, 2)
C = C[1]
annealing = Annealing(H, u2, coupling=C, bath=bath)
sol2 =  solve_ame(annealing, τ2, ω_hint = range(-15, 15, length = 200), alg=Tsit5(), reltol=1e-8, abstol=1e-8, saveat=range(0,τ2,length=100))

# solve stage 3
u3 = sol2[end]
H, C, cb = build_adiabatic_frame_evolution(hp, 3)
C = C[1]
annealing = Annealing(H, u3, coupling=C, bath=bath)
sol3 = solve_ame(annealing, 2 * τ1, ω_hint=range(-15, 15, length=200), alg=Tsit5(), reltol=1e-6, callback=cb, saveat=range(0,2*τ1,length=100))

# combine 3 stages of the evolution together for plotting
t_axis = [sol1.t; 2 * τ1 .+ sol2.t[2:end]; 2 * τ1 .+ τ2 .+ sol3.t[2:end]]
# in the adiabatic frame, the diagonal elements of the density matrix are
# the instantaneous eigen-populations
s1 = [real.(diag(u)) for u in sol1.u]
s2 = [real.(diag(u)) for u in sol2.u[2:end]]
s3 = [real.(diag(u)) for u in sol3.u[2:end]]
s_axis = hcat(s1..., s2..., s3...)
plot(t_axis/1e3, s_axis', xlabel=L"t\ (\mu s)", ylabel=L"P_i(t)", label=[L"E_0" L"E_1" L"E_2" L"E_3"], linewidth=2)
title!(latexstring(@sprintf "h_p = %.3f" hp))

# pre-calculate the spectral density of the
# hybrid Ohmic bath to speed up the computation
function hybrid_ohmic(ω, W, η, fc, T)
    ω = 2*π*ω
    bath = HybridOhmic(W, η, fc, T)
    y = [spectrum(w, bath) for w in ω]
    sfun = construct_interpolations(ω, y)
    CustomBath(spectrum = sfun)
end

η = 1.2732e-4; fc = 4; T = 12.5; W = 3;
bath_ohmic = Ohmic(η, fc, T)
bath_hybrid = hybrid_ohmic(range(-15, 15, length=1001), W, 100*η, fc, T)

H1, C1, cb1 = build_adiabatic_frame_evolution(hp, 1, polaron=true)
C1p = C1[1]
C1s = C1[2]
u1 = build_adiabatic_frame_u0(hp)
# For the adiabatic PTRE, we only perform the polaron transformation
# on the probe qubit. So there are two types of interactions:
interaction_1 = Interaction(C1p, bath_hybrid)
interaction_2 = Interaction(C1s, bath_ohmic)
interaction1 = InteractionSet(interaction_1, interaction_2)
annealing = Annealing(H1, u1, interactions = interaction1)
# we ignore the Lamb shift for the simulation;
# if the user would like to include the Lamb shift in the simulation
# it is better to provide an interpolated version of the `S`
# function via the keyword argument `lambshift_S`
sol1 = solve_ame(annealing, 2 * τ1, lambshift = false, alg=Tsit5(), reltol=1e-6, callback=cb1, saveat=range(0,2*τ1,length=100))

τ2 = 1e5
u2 = sol1[end]
H2, C2, cb2 = build_adiabatic_frame_evolution(hp, 2, polaron=true)
C2p = C2[1]
C2s = C2[2]
interaction_1 = Interaction(C2p, bath_hybrid)
interaction_2 = Interaction(C2s, bath_ohmic)
interaction2 = InteractionSet(interaction_1, interaction_2)
annealing = Annealing(H2, u2, interactions = interaction2)
sol2 = solve_ame(annealing, τ2, lambshift = false, alg=Tsit5(), reltol=1e-6, callback=cb2, saveat=range(0,τ2,length=100))

u3 = sol2[end]
H3, C3, cb3 = build_adiabatic_frame_evolution(hp, 3, polaron=true)
C3p = C3[1]
C3s = C3[2]
interaction_1 = Interaction(C3p, bath_hybrid)
interaction_2 = Interaction(C3s, bath_ohmic)
interaction3 = InteractionSet(interaction_1, interaction_2)
annealing = Annealing(H3, u3, interactions = interaction3)
sol3 = solve_ame(annealing, 2 * τ1, lambshift=false, alg=Tsit5(), reltol=1e-6, callback=cb3, saveat=range(0,2*τ1,length=100))

# combine 3 stages of the evolution for plotting
t_axis = [sol1.t; 2 * τ1 .+ sol2.t[2:end]; 2 * τ1 .+ τ2 .+ sol3.t[2:end]]
# in adiabatic frame, the diagonal elements of the density matrix are
# the instantaneous eigen-populations
s1 = [real.(diag(u)) for u in sol1.u]
s2 = [real.(diag(u)) for u in sol2.u[2:end]]
s3 = [real.(diag(u)) for u in sol3.u[2:end]]
s_axis = hcat(s1..., s2..., s3...)
plot(t_axis/1e3, s_axis', xlabel=L"t\ (\mu s)", ylabel=L"P_i(t)", label=[L"E_0" L"E_1" L"E_2" L"E_3"], linewidth=2)
title!(latexstring(@sprintf "h_p = %.3f" hp))


text:
using NumCME
using Test

@testset "CME model" begin
    include("test_propensity.jl")
    include("test_autodiff.jl")
    include("test_catalyst_interface.jl")
end

@testset "State Space" begin
    include("test_statespace.jl")
end

@testset "Transient CME" begin
    include("test_fspvec.jl")
    include("test_fspmat.jl")
    include("test_solver.jl")
end

@testset "Forward Sensitivity" begin
    include("test_sensmat.jl")
    include("test_sensfsp.jl")
end


text:
using YaoBlocks, LuxurySparse, SparseArrays, LinearAlgebra
using YaoBlocks.AD
using Test, Random

@testset "zero rand!" begin
    T = ComplexF64
    D = 10
    Random.seed!(2)
    for pm in [pmrand(T, D), Diagonal(randn(T, D)), sprand(T, D, D, 0.4)]
        z = zero(pm)
        @test z == zeros(D, D)
    end
end


text:
using Dolfin
using Base.Test

import Dolfin:dot, Function
import PyCall: PyDict

mesh = UnitSquareMesh(80, 80)
V = FunctionSpace(mesh, "Lagrange", 1)

# Define boundary condition
u0 =  Function(V)
bc =  DirichletBC(V, u0, "x[0] < DOLFIN_EPS || x[0] > 1.0 - DOLFIN_EPS")

# Define variational problem
u =  TrialFunction(V)
v =  TestFunction(V)
f =  Expression("10*exp(-(pow(x[0] - 0.5, 2) + pow(x[1] - 0.5, 2)) / 0.02)",
               degree=1)
g = Expression("sin(5*x[0])", degree=1)
a = inner(grad(u), grad(v))*dx
L = f*v*dx+ g*v*ds
# Define function for the solution
u = Function(V)
# Define error tolerance

# create a PyObject for parameters
solver_parameters=PyDict()

solver_parameters["linear_solver"]="lu"
solve(a==L, u, bc, solver_parameters=solver_parameters)
# Solve equation a = L with respect to u and the given boundary
# conditions, such that the estimated error (measured in M) is less
# than tol
plot(u, interactive=true)


text:
@testset "mapreduce" begin
  function maximum_avx(x)
    s = typemin(eltype(x))
    @turbo for i in eachindex(x)
      s = max(s, x[i])
    end
    s
  end
  function minimum_avx(x)
    s = typemax(eltype(x))
    @turbo for i in eachindex(x)
      s = min(s, x[i])
    end
    s
  end
  for T ∈ (Int32, Int64, Float32, Float64)
    @show T, @__LINE__
    if T <: Integer
      R = T(1):T(100)
      x7 = rand(R, 7)
      y7 = rand(R, 7)
      x = rand(R, 127, 7, 7)
      y = rand(R, 127, 7, 7)
      z = rand(R, 3, 3, 3, 3, 3, 3)
    else
      x7 = rand(T, 7)
      y7 = rand(T, 7)
      x = rand(T, 127, 7, 7)
      y = rand(T, 127, 7, 7)
      if VERSION ≥ v"1.4"
        @test vmapreduce(hypot, +, x, y) ≈ mapreduce(hypot, +, x, y)
        @test vmapreduce(^, (a, b) -> a + b, x7, y7) ≈ mapreduce(^, +, x7, y7)
      else
        @test vmapreduce(hypot, +, x, y) ≈ sum(hypot.(x, y))
        @test vmapreduce(^, (a, b) -> a + b, x7, y7) ≈ sum(x7 .^ y7)
      end
      z = rand(T, 3, 3, 3, 3, 3, 3)
    end
    @test vreduce(+, x7) ≈ sum(x7)
    @test vreduce(+, x) ≈ sum(x)
    if T === Int32
      @test vreduce(*, x7) == (prod(x7) % Int32)
      @test vreduce(*, x) == (prod(x) % Int32)
    else
      @test vreduce(*, x7) ≈ prod(x7)
      @test vreduce(*, x) ≈ prod(x)
    end
    @test vmapreduce(abs2, max, x) ≈ mapreduce(abs2, max, x)
    @test vmapreduce(abs2, min, x) ≈ mapreduce(abs2, min, x)
    @test vmapreduce(sqrt, *, x) ≈ mapreduce(sqrt, *, x)
    @test_throws AssertionError vmapreduce(hypot, +, x7, x)
    if VERSION ≥ v"1.4"
      @test vmapreduce(a -> 2a, *, x) ≈ mapreduce(a -> 2a, *, x)
      @test vmapreduce(sin, +, x7) ≈ mapreduce(sin, +, x7)
    else
      @test vmapreduce(a -> 2a, *, x) ≈ prod(2 .* x)
      @test vmapreduce(sin, +, x7) ≈ sum(sin.(x7))
    end
    @test vmapreduce(log, +, x) ≈ sum(log, x)
    @test vmapreduce(abs2, +, x) ≈ sum(abs2, x)
    @test maximum(x) == vreduce(max, x) == maximum_avx(x)
    @test minimum(x) == vreduce(min, x) == minimum_avx(x)

    @test vreduce(max, vec(x); dims = 1) == maximum(vec(x); dims = 1)
    @test vreduce(min, vec(x); dims = 1) == minimum(vec(x); dims = 1)
    @test vreduce(+, vec(x); dims = 1) ≈ sum(vec(x); dims = 1)
    for d = 1:ndims(x)
      @test vreduce(max, x; dims = d) == maximum(x; dims = d)
      @test vreduce(min, x; dims = d) == minimum(x; dims = d)
      @test vreduce(+, x; dims = d) ≈ sum(x; dims = d)
    end
    for d = 1:ndims(z)
      vreduce(max, z, dims = d) == maximum(z, dims = d)
      vreduce(min, z, dims = d) == minimum(z, dims = d)
      vreduce(+, z, dims = d) ≈ sum(z, dims = d)
    end
  end

end


text:
# basic C stuff

puts(s::Union{String,SubString{String}}) = ccall(:puts, Cint, (Ptr{Cchar},), s)

jl_malloc(n::Integer) = ccall(:jl_malloc, Ptr{Cvoid}, (Csize_t,), n)

# check if a function or C call failed

function check(ex::Expr, lock::Bool)
    ex.head == :call ||
        error("@check: not a call: $ex")
    if ex.args[1] == :ccall
        ex.args[2] isa QuoteNode ||
            error("@check: ccallee must be a symbol")
        f = ex.args[2].value :: Symbol
    else
        f = ex.args[1] :: Symbol
    end
    prefix = "$f: "
    ex = esc(ex)
    if lock
        ex = quote
            Base.iolock_begin()
            value = $ex
            Base.iolock_end()
            value
        end
    end
    quote
        r = $ex
        iszero(r) || @async @error($prefix * string(r))
        nothing
    end
end

macro check(ex::Expr) check(ex, false) end
macro check_iolock(ex::Expr) check(ex, true) end

# some libuv wrappers

const UV_READABLE = 1
const UV_WRITABLE = 2

function uv_poll_alloc()
    # allocate memory for: uv_poll_t struct + extra for curl_socket_t
    jl_malloc(Base._sizeof_uv_poll + sizeof(curl_socket_t))
end

function uv_poll_init(p::Ptr{Cvoid}, sock::curl_socket_t)
    @check_iolock ccall(:uv_poll_init, Cint,
        (Ptr{Cvoid}, Ptr{Cvoid}, curl_socket_t), Base.eventloop(), p, sock)
end

function uv_poll_start(p::Ptr{Cvoid}, events::Integer, cb::Ptr{Cvoid})
    @check_iolock ccall(:uv_poll_start, Cint,
        (Ptr{Cvoid}, Cint, Ptr{Cvoid}), p, events, cb)
end

function uv_poll_stop(p::Ptr{Cvoid})
    @check_iolock ccall(:uv_poll_stop, Cint, (Ptr{Cvoid},), p)
end

function uv_close(p::Ptr{Cvoid}, cb::Ptr{Cvoid})
    Base.iolock_begin()
    ccall(:uv_close, Cvoid, (Ptr{Cvoid}, Ptr{Cvoid}), p, cb)
    Base.iolock_end()
end

function uv_timer_init(p::Ptr{Cvoid})
    @check_iolock ccall(:uv_timer_init, Cint,
        (Ptr{Cvoid}, Ptr{Cvoid}), Base.eventloop(), p)
end

function uv_timer_start(p::Ptr{Cvoid}, cb::Ptr{Cvoid}, t::Integer, r::Integer)
    @check_iolock ccall(:uv_timer_start, Cint,
        (Ptr{Cvoid}, Ptr{Cvoid}, UInt64, UInt64), p, cb, t, r)
end

function uv_timer_stop(p::Ptr{Cvoid})
    @check_iolock ccall(:uv_timer_stop, Cint, (Ptr{Cvoid},), p)
end

# additional libcurl methods

import LibCURL: curl_multi_socket_action

function curl_multi_socket_action(multi_handle, s, ev_bitmask)
    curl_multi_socket_action(multi_handle, s, ev_bitmask, Ref{Cint}())
end

# curl string list structure

struct curl_slist_t
    data::Ptr{Cchar}
    next::Ptr{curl_slist_t}
end


text:
using ScenTreesMakie
using Test

@testset "Alois.jl" begin
    @testset "Predefined tree - Tree 402" begin
        a = Tree(402)
        @test length(a.parent) == 15
        @test length(a.state) == length(a.probability) == length(a.parent) == 15
        @test sum(a.probability) == 8.0
        @test length(a.children) == 8
        @test length(root(a)) == 1
    end
    @testset "Initial Trees" begin
        init = Tree([1,2,2,2])
        @test length(init.parent) == 15
        @test length(init.state) == length(init.probability) == length(init.parent) == 15
        @test length(init.children) == 8
        @test length(stage(init)) == 15
        @test height(init) == 3
        @test length(leaves(init)) == 3
        @test nodes(init) == 1:15
        @test length(nodes(init)) == 15
        @test length(root(init)) == 1
    end
    @testset "A sample of a Scenario Tree 1D" begin
        x = Tree([1,3,3,3,3])
        @test length(x.parent) == 121
        @test length(x.state) == length(x.probability) == length(x.parent) == 121
        @test length(x.children) == 41
        @test length(root(x)) == 1
        @test length(leaves(x)) == 3
    end
    @testset "Sample stochastic functions" begin
        a = gaussian_path()
        b = running_maximum()
        @test length(a) == 4
        @test length(b) == 4
    end
    @testset "ScenTrees.jl - Tree Approximation 1D" begin
        paths = [gaussian_path,running_maximum]
        trees = [Tree([1,2,2,2]),Tree([1,3,3,3])]
        samplesize = 1001
        p = 2
        r = 2
        plot(trees[1])
        for path in paths
            for newtree in trees
                tree_approximation!(newtree,path,samplesize,p,r)
                @test length(newtree.parent) == length(newtree.state)
                @test length(newtree.parent) == length(newtree.probability)
                @test length(stage(newtree)) == length(newtree.parent)
                @test height(newtree) == maximum(stage(newtree))
                @test round(sum(leaves(newtree)[3]),digits=1) == 1.0   #sum of unconditional probabilities of the leaves
                @test length(root(newtree)) == 1
            end
        end
    end

    @testset "ScenTrees.jl - Lattice Approximation" begin
        tstLat = lattice_approximation([1,2,3,4],gaussian_path,20001,2,1)
        @test length(tstLat.state) == length(tstLat.probability)
        @test round.(sum.(tstLat.probability), digits = 1)  == [1.0, 1.0, 2.0, 3.0] #sum of probs at every stage
    end

end

@testset "LatticeStructure" begin
    lat = Lattice(302)
    @test length(lat.state) == 3

    lat = Lattice(303)
    @test length(lat.state) == 3

    lat = Lattice(304)
    @test length(lat.state) == 4

    name = "Lattice 1x2x2"
    state = [[0.0],[-0.856;0.933],[-1.27;1.307]]
    probability = [[1.0],[0.52;0.48],[0.75 0.25;0.26 0.74]]

    lat = Lattice(name,state,probability)
    @test lat.name == Lattice(302).name
    @test lat.state == Lattice(302).state

    @test structure(lat)[1] == [1,2,2]
    @test structure(lat)[2] == 5
    @test structure(lat)[3] == 6
    @test structure(lat)[4] == 4

end

@testset "Approx" begin
    lat = Lattice(304);

    @test length(path_ident(lat, gaussian_path, 1)[1]) == 4

    trr = Tree(404)
    @test length(path_ident(trr,gaussian_path, 1)[1]) == 4
end

@testset "Kernel" begin
    using Distributions
    gsdata = Array{Float64}(undef,1000,4)
    for i = 1:1000
           gsdata[i,:] = gaussian_path()
    end
    @test length(kernel_scenarios(gsdata,Logistic; Markovian = true)()) == 4
    @test length(kernel_scenarios(gsdata,Logistic; Markovian = false)()) == 4
end

@testset "helper" begin
    probs = [1.0, 2.0, 3.0, 4.0];
    ontoSimplex!(probs)
    @test sum(probs) == 1.0
    probs = [1.0, -2.0, 3.0, 4.0];
    ontoSimplex!(probs)
    @test sum(probs) == 1.0

    trr = Tree(404);
    ScenTreesMakie.summary(trr)

end

@testset "TreeStructure" begin
    @test Tree(0).parent == zeros(Int64 , 0)
    @test Tree(302).parent ==[0,1,1,2,2,3,3]
    @test Tree(303).parent == [0,1,2,2,2,2]
    @test Tree(304).parent == [0,1,2,0,4,5,0,7,8,0,10,11]
    @test Tree(305).parent == [0,1,2,2,2,2]
    @test Tree(306).parent == [0,1,1,2,2,3,3]
    @test Tree(307).parent == [0,1,1,1,1,2,3,4,5]
    @test Tree(401).parent == [0,1,2,2,3,3,4,4]
    @test Tree(402).parent == [0,1,1,2,2,3,3,4,4,5,5,6,6,7,7]
    @test Tree(404).parent == [0,1,1,2,2,3,3,4,4,5,5,6,6,7,7]
    @test Tree(405).parent == [0,1,2,3,3,2,6,6,6,1,10,10,11,12,12,12]

    @test children(2) == [[1]]

    trr = Tree(404)
    @test stage(trr, 2) == [1]

    @test length(leaves(trr,2)[1]) == 4

end

@testset "WassersteinSkript" begin

p1= [.3, 0.7]
s1= [-.7, 1.1]
p2= [0.2, 0.3, 0.2, 0.1, 0.2]
s2= [-.3, .9, 1.7, 3.2, 1.9]
d= distFunction(s1, s2)

@show A= Wasserstein(p1, p2, distFunction(s1, s2), 2.)

@test A.distance ≈ 0.9679876032264049
@test A.π ≈ [0.2 0.1 0.0 0.0 0.0; 0.0 0.2 0.2 0.1 0.2]

println()
@show A= Sinkhorn(p1, p2, distFunction(s1, s2), 2., 4.3)
@test A.distance ≈ 0.967988 atol= 1e-4
@test A.π ≈ [0.2 0.1 0.0 0.0 0.0; 0.0 0.2 0.2 0.1 0.2] atol = 1e-6

println()
@show A= Sinkhorn(p1, p2, distFunction(s1, s2), 2., 1.)
@test A.distance ≈ 0.9923876770989377 atol = .01
@test A.π ≈ [0.2 0.1 0.0 0.0 0.0; 0.0 0.2 0.2 0.1 0.2] atol = 0.02

end

@testset "nested distance" begin
    trr1 = Tree([1,2,2,2]); trr2 = Tree([1,2,2,2]);
    trr1.state = [0,0,0,0,0,0,0,1,3,20,40,2,4,21,41];
    trr2.state = [0,0,0,0,0,0,0,1,21,3,41,2,20,4,40];

    trr1.probability = [1.0,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5];
    trr2.probability = trr1.probability;

    @test nestedWasserstein(trr1,trr2,1.0) == 8.75

end


text:
function volume_sim(config)
    (;n_shapes,n_sim,n_dims) = config
    (;n_dims,n_part,n_cells, n_start) = config
    # partition boundaries
    bounds = fill((0,1), n_dims)
    indices = set_indices(n_dims, n_part, n_cells)
    p_bounds = range(0, 1, length = n_part + 1)

    odd_shape = OddShape(indices, p_bounds, n_cells)

    sample(bounds) = map(b -> rand(Uniform(b...)), bounds)

    init_parms = map(_ -> sample(bounds), 1:n_start)

    options = Options(;
        radius = .10,
        bounds,
        n_iters = 20_000,
        init_parms
    )

    results = find_partitions(
        model,
        p_fun,
        options,
        odd_shape
    )

    df = DataFrame(results)

    parm_names = Symbol.("p",1:n_dims)
    transform!(
        df,
        :parms => identity => parm_names
    )

    transform!(df, :pattern => denserank => :pattern_id)

    target_pattern = 1
    df_p = filter(x-> x.pattern == target_pattern, df)
    points = Array(df_p[:,parm_names])
    v1 = estimate_volume(
        model,
        p_fun,
        points,
        target_pattern,
        bounds,
        odd_shape;
        n_sim,
    )

    target_pattern = 2
    df_p = filter(x-> x.pattern == target_pattern, df)
    points = Array(df_p[:,parm_names])
    v2 = estimate_volume(
        model,
        p_fun,
        points,
        target_pattern,
        bounds,
        odd_shape;
        n_sim
    )
    return v1 / v2
end

max_stars_repo_name: itsdfish/ParameterSpacePartitions

max_forks_repo_name: itsdfish/ParameterSpacePartitions


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

mutable struct IoK8sApiCoreV1WindowsSecurityContextOptions <: SwaggerModel
    gmsaCredentialSpec::Any # spec type: Union{ Nothing, String } # spec name: gmsaCredentialSpec
    gmsaCredentialSpecName::Any # spec type: Union{ Nothing, String } # spec name: gmsaCredentialSpecName
    runAsUserName::Any # spec type: Union{ Nothing, String } # spec name: runAsUserName

    function IoK8sApiCoreV1WindowsSecurityContextOptions(;gmsaCredentialSpec=nothing, gmsaCredentialSpecName=nothing, runAsUserName=nothing)
        o = new()
        validate_property(IoK8sApiCoreV1WindowsSecurityContextOptions, Symbol("gmsaCredentialSpec"), gmsaCredentialSpec)
        setfield!(o, Symbol("gmsaCredentialSpec"), gmsaCredentialSpec)
        validate_property(IoK8sApiCoreV1WindowsSecurityContextOptions, Symbol("gmsaCredentialSpecName"), gmsaCredentialSpecName)
        setfield!(o, Symbol("gmsaCredentialSpecName"), gmsaCredentialSpecName)
        validate_property(IoK8sApiCoreV1WindowsSecurityContextOptions, Symbol("runAsUserName"), runAsUserName)
        setfield!(o, Symbol("runAsUserName"), runAsUserName)
        o
    end
end # type IoK8sApiCoreV1WindowsSecurityContextOptions

const _property_map_IoK8sApiCoreV1WindowsSecurityContextOptions = Dict{Symbol,Symbol}(Symbol("gmsaCredentialSpec")=>Symbol("gmsaCredentialSpec"), Symbol("gmsaCredentialSpecName")=>Symbol("gmsaCredentialSpecName"), Symbol("runAsUserName")=>Symbol("runAsUserName"))
const _property_types_IoK8sApiCoreV1WindowsSecurityContextOptions = Dict{Symbol,String}(Symbol("gmsaCredentialSpec")=>"String", Symbol("gmsaCredentialSpecName")=>"String", Symbol("runAsUserName")=>"String")
Base.propertynames(::Type{ IoK8sApiCoreV1WindowsSecurityContextOptions }) = collect(keys(_property_map_IoK8sApiCoreV1WindowsSecurityContextOptions))
Swagger.property_type(::Type{ IoK8sApiCoreV1WindowsSecurityContextOptions }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiCoreV1WindowsSecurityContextOptions[name]))}
Swagger.field_name(::Type{ IoK8sApiCoreV1WindowsSecurityContextOptions }, property_name::Symbol) =  _property_map_IoK8sApiCoreV1WindowsSecurityContextOptions[property_name]

function check_required(o::IoK8sApiCoreV1WindowsSecurityContextOptions)
    true
end

function validate_property(::Type{ IoK8sApiCoreV1WindowsSecurityContextOptions }, name::Symbol, val)
end


text:
type KalmanFiltered{T}
    filtered::Array{T}
    predicted::Array{T}
    error_cov::Array{T}
    pred_error_cov::Array{T}
    model::StateSpaceModel
    y::Array{T}
    u::Array{T}
    loglik::T
end

function Base.show{T}(io::IO, filt::KalmanFiltered{T})
    n = size(filt.y, 1)
    dx, dy = filt.model.nx, filt.model.ny
    println("KalmanFiltered{$T}")
    println("$n observations, $dx-D process x $dy-D observations")
    println("Negative log-likelihood: $(filt.loglik)")
end

function kalman_filter{T}(y::Array{T}, model::StateSpaceModel{T}; u::Array{T}=zeros(size(y,1), model.nu))

    @assert size(u,1) == size(y,1)
    @assert size(y,2) == model.ny
    @assert size(u,2) == model.nu

    function kalman_recursions(y_i::Vector{T}, u_i::Vector{T},
                                  C_i::Matrix{T}, D_i::Matrix{T}, W_i::Matrix{T},
                                  x_pred_i::Vector{T}, P_pred_i::Matrix{T})
        if !any(isnan(y_i))
            innov =  y_i - C_i * x_pred_i - D_i * u_i
            S = C_i * P_pred_i * C_i' + W_i  # Innovation covariance
            K = P_pred_i * C_i' / S # Kalman gain
            x_filt_i = x_pred_i + K * innov
            P_filt_i = (I - K * C_i) * P_pred_i
            dll = (dot(innov,S\innov) + logdet(S))/2
        else
            x_filt_i = x_pred_i
            P_filt_i = P_pred_i
            dll = 0
        end
        return x_filt_i, P_filt_i, dll
    end #kalman_recursions

    y = y'
    u = u'
    n = size(y, 2)
    x_pred = zeros(model.nx, n)
    x_filt = zeros(x_pred)
    P_pred = zeros(model.nx, model.nx, n)
    P_filt = zeros(P_pred)
    log_likelihood = n*model.ny*log(2pi)/2

    # first iteration
    A_1 = model.A(1)
    x_pred[:, 1] = model.x1
    P_pred[:, :, 1] = model.P1
    x_filt[:, 1], P_filt[:,:,1], dll = kalman_recursions(y[:, 1], u[:, 1],
                                            model.C(1), model.D(1), model.W(1),
                                            model.x1, model.P1)
    log_likelihood += dll

    for i=2:n
        A_i1 = model.A(i-1)
        x_pred[:, i] =  A_i1 * x_filt[:, i-1] + model.B(i-1) * u[:, i-1]
        P_pred[:, :, i] = A_i1 * P_filt[:, :, i-1] * A_i1' + model.V(i-1)
        x_filt[:, i], P_filt[:,:,i], dll = kalman_recursions(y[:, i], u[:, i],
                                                model.C(i), model.D(i), model.W(i),
                                                x_pred[:,i], P_pred[:,:,i])
        log_likelihood += dll
    end

    return KalmanFiltered(x_filt', x_pred', P_filt, P_pred, model, y', u', log_likelihood)
end


text:
using DSP
"""
    Filter(filters)

Apply filter to the signal

Inputs
------
* `filters` array of DSP filter objects.

Example
-------
```julia
using DSP
responsetype = Bandpass(500, 4000; fs=48000)
designmethod = Butterworth(4)
zpg = digitalfilter(responsetype, designmethod)
f_left = DSP.Filters.DF2TFilter(zpg)
f_right = DSP.Filters.DF2TFilter(zpg)

bandpass = AuditoryStimuli.Filter([f_left, f_right])
filtered_sound = modify(bandpass, original_sound)
```
"""
mutable struct Filter
    filter::Array{DSP.DF2TFilter, 1}  # Filter object
    enable::Bool
    num_filters::Int
end
function Filter(filters::Any)
    AuditoryStimuli.Filter(filters, true, length(filters))
end
function modify(sink::Filter, buf)
    if sink.enable
        for idx in 1:sink.num_filters
            buf.data[:, idx] = DSP.filt(sink.filter[idx], buf.data[:, idx])
        end
    end
    return buf
end


text:
include("SimulationsUtils.jl")
include("MethodUtils.jl")
include("WassersteinVariance.jl")
using Serialization

function benchmark(nsim, numS, ne=500, noisetype="normal", numvar=18, seed=123)
    Random.seed!(seed)

    ######################## SIMULATION ##############################
    # nsim = 100
    # These will store a four tuple for each simulation: nfp, nfn, length(S), length(Ŝ)
    ols_errs = []
    lingam_errs = []
    gies_errs = []
    icp_errs = []
    wvm_errs = []

    simtype = "numS-$(numS)/ne-$(ne)/$(noisetype)"
    for i in 1:nsim
        ############### DATA SIM #################
        # collect simulation data
        X_full, S, nₑs = deserialize("../data/Simulations/$(simtype)/sim-$(i).jls") #"Simulations/$(noisetype)/numS-$(numS)/sim-$(i)$(addedinfo).jls"
        n, p = size(X_full,1), size(X_full,2)
        X, y = X_full[:, 1:end-1], X_full[:, end]
        E = length(nₑs)

        ################# OLS #################
        mm = lm(hcat(X, ones(size(X,1))),y)
        Ŝ_ols = sig_pvals_ols(mm, 0.1/(p-1))
        push!(ols_errs, error_count(S, Ŝ_ols))

        ################# LinGam #################
        Ŝ_lingam = lingam(X_full)
        push!(lingam_errs, error_count(S, Ŝ_lingam))

        ################# GIES #################
        Ŝ_gies = gies(X_full, E, nₑs[1])
        push!(gies_errs, error_count(S, Ŝ_gies))

        ################# ICP #################
        α = 0.1
        icp_res = icp(X, y, nₑs, α, true)
        Ŝ_icp = findall_nonzero(isnothing(icp_res) ? zeros(p-1) : icp_res)
        push!(icp_errs, error_count(S, Ŝ_icp))

        ################# lasso #################
        cf = select_lasso_by_num_var(X, y, numvar)
        select = 1:p-1
        select = filter(x -> cf[x] != 0, select)

        ################# WVM #################
        wv = DirectWassersteinVariance(X[:,select], y, nₑs, nₑs / n)
        res = WV_screening(wv, α=α, solver=LBFGS(m=100), B=50)
        Ŝ_wvm = select[findall_nonzero(res[2])]
        push!(wvm_errs, error_count(S, Ŝ_wvm))

        println("###### SIMULATION $i IS DONE ! ######")
    end
    run(`mkdir -p ../data/results/$(simtype)`)
    serialize("../data/results/$(simtype)/ols.jls", ols_errs)
    serialize("../data/results/$(simtype)/lingam.jls", lingam_errs)
    serialize("../data/results/$(simtype)/gies.jls", gies_errs)
    serialize("../data/results/$(simtype)/icp.jls", icp_errs)
    serialize("../data/results/$(simtype)/wvm.jls", wvm_errs)
end
benchmark(100, 6, 500, "normal")
benchmark(100, 6, 500, "mixed")


text:
function deleteat(M::Array, i1::Integer, i2::Integer)
    return [M[1:i1 - 1,1:i2 - 1] M[1:i1 - 1,i2 + 1:end];M[i1 + 1:end,1:i2 - 1] M[i1 + 1:end,i2 + 1:end]]
end

deleteat(M::Array,i::Integer) = deleteat(M, i, i)

function orthogonalrows(axis::AbstractVector)
    if norm(axis) > 0
        axis = normalize(axis)
    end
    A = svd(skew(axis)).Vt
    inds = SA[1; 2; 3]
    V1 = A[1,inds]'
    V2 = A[2,inds]'
    V3 = axis' # instead of A[3,:] for correct sign: abs(axis) = abs(A[3,:])

    return V1, V2, V3
end

function orthogonalcols(axis::AbstractVector)
    V1, V2, V3 = orthogonalrows(axis)
    return V1', V2', V3'
end

@inline function skewplusdiag(v::AbstractVector{T},w::T) where T
    SA[
         w    -v[3]  v[2]
         v[3]  w    -v[1]
        -v[2]  v[1]  w
    ]
end

function getfieldnumber(obj)
    i = 1
    while true
        !isdefined(obj, i) ? break : (i+=1)
    end
    return i-1
end

@inline offsetrange(offset, length) = (offset-1)*length+1:offset*length
@inline offsetrange(offset, length, totallength, inneroffset) = (offset-1)*totallength+(inneroffset-1)*length+1:(offset-1)*totallength+inneroffset*length


text:
using OptimPack
using Base.Test

include("rosenbrock.jl")