text:
import NLPModelsIpopt: ipopt

"""
    `ipopt(stp::NLPStopping; subsolver_verbose::Int = 0, kwargs...)`

Stopping-version of the `ipopt` function from NLPModelsIpopt.jl.
This function calls `fill_in!` (doesn't update hessian) and `stop!` after `ipopt` call,
if the problem is a success (`:first_order` or `:acceptable`) and `fill_in_on_success` is true or if it failed and `fill_in_on_failure` is true.

`subsolver_verbose` corresponds to `print_level` argument in `ipopt`.
Other keyword arguments are passed to the `ipopt` call.
Selection of possible [options](https://coin-or.github.io/Ipopt/OPTIONS.html#OPTIONS_REF).
"""
function NLPModelsIpopt.ipopt(
  stp::NLPStopping;
  subsolver_verbose::Int = 0,
  fill_in_on_success = true,
  fill_in_on_failure = true,
  kwargs...,
)
  stp.meta.start_time = time()
  #xk = solveIpopt(stop.pb, stop.current_state.x)
  nlp = stp.pb
  stats = ipopt(
    nlp;
    print_level = subsolver_verbose,
    tol = stp.meta.rtol,
    x0 = stp.current_state.x,
    max_iter = stp.meta.max_iter,
    max_cpu_time = stp.meta.max_time,
    dual_inf_tol = stp.meta.atol,
    constr_viol_tol = stp.meta.atol,
    compl_inf_tol = stp.meta.atol,
    kwargs...,
  )

  #Update the meta boolean with the output message
  #=
  if stats.status == :first_order stp.meta.suboptimal      = true end
  if stats.status == :acceptable  stp.meta.suboptimal      = true end
  if stats.status == :infeasible  stp.meta.infeasible      = true end
  if stats.status == :small_step  stp.meta.stalled         = true end
  if stats.status == :max_eval    stp.meta.max_eval        = true end
  if stats.status == :max_iter    stp.meta.iteration_limit = true end
  if stats.status == :max_time    stp.meta.tired           = true end
  if stats.status ∈ [:neg_pred,
                     :not_desc]   stp.meta.fail_sub_pb     = true end
  if stats.status == :unbounded   stp.meta.unbounded       = true end
  if stats.status == :user        stp.meta.stopbyuser      = true end
  if stats.status ∈ [:stalled,
                     :small_residual,
                     :small_step]   stp.meta.stalled       = true end
  #if stats.status == :exception   stp.meta.exception      = true end #available ≥ 0.2.6
  =#
  stp = stats_status_to_meta!(stp, stats)

  if status(stp) == :Unknown
    @warn "Error in StoppingInterface statuses: return status is $(stats.status)"
    @show stats.solver_specific
  end

  stp.meta.nb_of_stop = stats.iter
  Stopping._update_time!(stp.current_state, time()) # stats.elapsed_time

  success = stats.status ∈ [:first_order, :acceptable]
  if (success && fill_in_on_success) || (!success && fill_in_on_failure)
    x = stats.solution
    # Not mandatory, but in case some entries of the State are used to stop
    fill_in!(stp, x, Hx = stp.current_state.Hx)

    stop!(stp)
  end

  return stp
end


text:
"""
Subspace of particles of states
"""
mutable struct MCVISubspace{S,A}
    particles::Vector{S}
    imm_rewards::Vector{Reward}
    next_state::Dict{A,MCVISubspace{S,A}}
end
"""
MCVISubspace constructor
"""
MCVISubspace(p::Vector{S}, r, ns::Dict{A,MCVISubspace{S,A}}) where {S,A} = MCVISubspace(p, r, ns)

"""
Number of particles in the subspace
"""
Base.length(ss::MCVISubspace) = length(ss.particles)

"""
    particle(ss::MCVISubspace, i::Int64)

Returns i-th particle for the given subspace
"""
particle(ss::MCVISubspace, i::Int64) = ss.particles[i]

"""
    weights(ss::MCVISubspace{S,A}, obs::O), pomdp::POMDPs.POMDP

Returns an observation's weight for each particle in the subspace
"""
function weights(ss::MCVISubspace{S,A}, obs::O, pomdp::POMDPs.POMDP) where {S,A,O}
    wts = zeros(length(ss))
    for i in 1:length(wts)
        wts[i] = obs_weight(pomdp, particle(ss, i), obs)
    end
    return wts
end

"""
    create_next(ss::MCVISubspace{S,A}, act::A, pomdp::POMDPs.POMDP{S,A,O})

Creates the next subspace for an action
"""
function create_next(ss::MCVISubspace{S,A}, act::A, pomdp::POMDPs.POMDP, rng::AbstractRNG) where {S,A}
    next_particles = Vector{S}(undef, length(ss))
    imm_rs = zeros(length(ss))

    for (i, s) in enumerate(ss.particles)
        (next_particles[i], imm_rs[i]) = gen(DDNOut(:sp,:r), pomdp, s, act, rng)
    end
    return MCVISubspace(next_particles, imm_rs, Dict{A, MCVISubspace{S,A}}())
end

"""
    next(ss::MCVISubspace{S,A}, act::A, pomdp::POMDPs.POMDP{S,A,O})

Returns the next subspace according to the action. If it does not exist, creates one.
"""
function next(ss::MCVISubspace{S,A}, act::A, pomdp::POMDPs.POMDP, rng::AbstractRNG) where {S,A}
    local next_ss::MCVISubspace{S,A}
    if haskey(ss.next_state, act)
        next_ss = ss.next_state[act]
    else
        next_ss = create_next(ss, act, pomdp, rng)
        ss.next_state[act] = next_ss
    end
    return next_ss
end

"""
    reward(ss::MCVISubspace{S,A}, weights::Vector{Float64})

Returns the weighted immediate reward from a subspace
"""
function reward(ss::MCVISubspace, weights::Vector{Float64})
    if isempty(ss.imm_rewards)
        return 0
    end
    return dot(weights, ss.imm_rewards)
end

ext: jl


text:
using SimpleTraits

import Base: display

type RotoloDisplay <: Display ; end
pushdisplay(RotoloDisplay())

@traitdef ShouldRedirect{X}
@traitimpl ShouldRedirect{X} <- isredirected(X)

@traitfn display{X; ShouldRedirect{X}}(::RotoloDisplay, x::X)  = showmsg(x)
@traitfn display{X; !ShouldRedirect{X}}(::RotoloDisplay, x::X) = display(x)

macro redirect(args...)
  currentSession == nothing &&
    error("A session needs to be created before calling @redirect")

  for a in args
    t = try
          eval(current_module(), a)
        catch e
          error("can't evaluate $a, error $e")
        end
    isa(t, Type) || error("$a does not evaluate to a type")

    push!(currentSession.redirected_types, t)
    loadmsg(t)

    # @eval display(::RotoloDisplay, x::($t)) = showmsg(x)
  end
  nothing
end


text:
##########################
# FunctionCell Interface #
##########################
"""
	loss([x, f_targets=identity, f_output=identity, losstype=LossFunctions.L2DistLoss(), avgmode=LossFunctions.AvgMode.Mean()])

Computes the loss of `x`. Here `x` is assumend to contain information on both the output of a classifier or regressor as well as
the target information e.g. labels or regression values.

# Arguments
  * `x` is the data. Can be a labeled `DataCell`, `Tuple` of the form `(<output data>, <target data>)` or an `Array` of outputs; in this latter case
`f_targets` needs to be specified
  * `f_targets` is a function that is applied to the labels of `x`, before applying the loss function
  * `f_output` is a function that is applied to the data of `x`, before applying the loss function
  * `losstype`::LeanBase.Loss is the loss object describing the type of loss applied
  * `avgmode`::LossFunctions.AverageMode is the loss aggregation object which describes how individual observation losses are combined

It is assumed that the inputs of `f_targets` and `f_output` are known, as well as their outputs. Generally, only one of the functions is
actually needed, transforming either targets or output to a common format.

Check the `LossFunctions.jl` and `MLLabelUtils.jl` documentation for further details.

# Examples
```
# Train a small classifier and check the mean squared error
julia> using j4pr; tr, ts = splitobs(shuffleobs(DataGenerator.iris())); LC = lindisc(); LCT = tr |> LC;

julia> LC
Linear discriminant classifier: r1=0.0, r2=0.0, no I/O size information, untrained

julia> LCT
Linear discriminant classifier: r1=0.0, r2=0.0, 4 -> 3, trained

julia> ff1=x->MLLabelUtils.convertlabel(MLLabelUtils.LabelEnc.OneOfK,x,LCT.y["labels"]) # ff1 processes existing labels
(::#43) (generic function with 1 method)

julia> ff2=x->MLLabelUtils.convertlabel(MLLabelUtils.LabelEnc.OneOfK, targets(indmax,x)) #ff2 processes the outputs
(::#47) (generic function with 1 method)

julia> ff1(-ts)
3×45 Array{Int64,2}:
 1  1  0  1  1  0  0  0  1  0  1  0  0  1  0  1  0  0  1  0  1  1  0  1  1  0  1  0  0  0  0  0  0  0  0  1  0  0  0  1  1  0  1  1  1
 0  0  1  0  0  1  1  1  0  1  0  0  0  0  0  0  1  0  0  0  0  0  0  0  0  0  0  0  0  1  0  0  0  0  1  0  0  1  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  1  1  0  1  0  0  1  0  1  0  0  1  0  0  1  0  1  1  0  1  1  1  1  0  0  1  0  1  0  0  1  0  0  0

julia> ff2(+ts|>LCT)
3×45 Array{Int64,2}:
 1  1  0  1  1  0  0  0  1  0  1  0  0  1  0  1  0  0  1  0  1  1  0  1  1  0  1  0  0  0  0  0  0  0  0  1  0  0  0  1  1  0  1  1  1
 0  0  1  0  0  1  1  1  0  1  0  0  0  0  0  0  1  0  0  0  0  0  0  0  0  0  0  1  0  1  0  0  0  0  1  0  0  1  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  1  1  0  1  0  0  1  0  1  0  0  1  0  0  1  0  0  1  0  1  1  1  1  0  0  1  0  1  0  0  1  0  0  0

# Calculating the error:
julia> j4pr.loss(ts|> LCT, ff1, ff2)
0.014814814814814815
```
"""
# Create fixed function cell (no `x`)
loss(f_targets::F=identity, f_output::G=identity,
	losstype::L=LossFunctions.L2DistLoss(), avgmode::A=LossFunctions.AvgMode.Mean()) where {F<:Function, G<:Function, L<:LearnBase.Loss, A<:LossFunctions.AverageMode} =
begin

	FunctionCell(loss, (f_targets, f_output, losstype, avgmode), "Loss: type=$(losstype), avgmode=$(avgmode)")
end

# Work with labeled data cells
loss(x::T where T<:CellDataL, f_targets::F=identity, f_output::G=identity,
     losstype::L=LossFunctions.L2DistLoss(), avgmode::A=LossFunctions.AvgMode.Mean()) where {F<:Function, G<:Function, L<:LearnBase.Loss, A<:LossFunctions.AverageMode} =
begin
	return loss(f_targets(gety!(x)), f_output((getx!(x))), losstype, avgmode)
	#       	y         	     ŷ
end

# Work with Tuples
loss(x::T where T<:Tuple{U,V} where {U,V}, f_targets::F=identity, f_output::G=identity,
     losstype::L=LossFunctions.L2DistLoss(), avgmode::A=LossFunctions.AvgMode.Mean()) where {F<:Function, G<:Function, L<:LearnBase.Loss, A<:LossFunctions.AverageMode} =
begin
	return loss(f_targets(x[2]), f_output(x[1]), losstype, avgmode)
	#                y               ŷ
end

# Work with Array input (output); In this case labels need to be specified through f_targets (e.g. f_targets=()->Labels::Array{...})
loss(x::T where T<:AbstractArray, f_targets::F, f_output::G=identity,
     losstype::L=LossFunctions.L2DistLoss(), avgmode::A=LossFunctions.AvgMode.Mean()) where {F<:Function, G<:Function, L<:LearnBase.Loss, A<:LossFunctions.AverageMode} =
begin
	return loss(f_targets(), f_output(x), losstype, avgmode)
	#                y               ŷ
end

#Define low level loss function
loss(y::T where T<:AbstractArray, ŷ::S where S<:AbstractArray, losstype::U where U<:LearnBase.Loss, avgmode::V where V<:LossFunctions.AverageMode ) = LossFunctions.value(losstype, y, ŷ, avgmode)


text:
# Tests based on the Vector API tutorial found at http://gdal.org/ogr_apitut.html

# Reading From OGR

# TODO: provide a convenient way to construct nOpenFlags
# nOpenFlags: a combination of GDAL_OF_ flags that may be combined through logical or operator
# Driver kind: GDAL_OF_RASTER for raster drivers, GDAL_OF_VECTOR for vector drivers. If none of the value is specified, both kinds are implied.
# Access mode: GDAL_OF_READONLY (exclusive)or GDAL_OF_UPDATE.
# Shared mode: GDAL_OF_SHARED. If set, it allows the sharing of GDALDataset handles for a dataset with other callers that have set GDAL_OF_SHARED. In particular, GDALOpenEx() will first consult its list of currently open and shared GDALDataset's, and if the GetDescription() name for one exactly matches the pszFilename passed to GDALOpenEx() it will be referenced and returned, if GDALOpenEx() is called from the same thread.
# Verbose error: GDAL_OF_VERBOSE_ERROR. If set, a failed attempt to open the file will lead to an error message to be reported.

# file from GDAL SVN
dataset = GDAL.openex("data/point.geojson", GDAL.GDAL_OF_VECTOR, C_NULL, C_NULL, C_NULL)

@test GDAL.datasetgetlayercount(dataset) == 1
layer = GDAL.datasetgetlayer(dataset, 0)
# GDAL 2.2 GeoJSON driver change
layername = GDAL.GDALVERSION >= v"2.2.0" ? "point" : "OGRGeoJSON"
layerbyname = GDAL.datasetgetlayerbyname(dataset, layername)
@test layer == layerbyname
@test GDAL.getname(layer) == layername

GDAL.resetreading(layer)

featuredefn = GDAL.getlayerdefn(layer)
@test GDAL.getfieldcount(featuredefn) == 2

fielddefn = GDAL.getfielddefn(featuredefn, 0)
@test GDAL.gettype(fielddefn) == GDAL.OFTReal
fielddefn = GDAL.getfielddefn(featuredefn, 1)
@test GDAL.gettype(fielddefn) == GDAL.OFTString

feature = GDAL.getnextfeature(layer) # first feature
@test GDAL.getfieldasdouble(feature, 0) == 2.0
@test GDAL.getfieldasstring(feature, 1) == "point-a"
feature = GDAL.getnextfeature(layer) # second feature
@test GDAL.getfieldasdouble(feature, 0) == 3.0
@test GDAL.getfieldasstring(feature, 1) == "point-b"

geometry = GDAL.getgeometryref(feature)
@test GDAL.getgeometryname(geometry) == "POINT"
@test GDAL.getgeometrytype(geometry) == GDAL.wkbPoint
@test GDAL.getgeomfieldcount(featuredefn) == 1
@test GDAL.getx(geometry, 0) == 100.2785
@test GDAL.gety(geometry, 0) == 0.0893

# export to WKT
wkt_ptr = Ref(Cstring(C_NULL))
@test GDAL.exporttowkt(geometry, wkt_ptr) == GDAL.OGRERR_NONE
@test unsafe_string(wkt_ptr[]) == "POINT (100.2785 0.0893)"
GDAL.C.OGRFree(pointer(wkt_ptr[]))

GDAL.destroy(feature)
GDAL.close(dataset)

# Writing to OGR
pointshapefile = joinpath("tmp", "point_out")
driver = GDAL.getdriverbyname("ESRI Shapefile")
dataset = GDAL.create(driver, "$pointshapefile.shp", 0, 0, 0, GDAL.GDT_Unknown, C_NULL)
nosrs = convert(Ptr{GDAL.OGRSpatialReferenceH}, C_NULL)
layer = GDAL.datasetcreatelayer(dataset, "point_out", nosrs, GDAL.wkbPoint, C_NULL)
fielddefn = GDAL.fld_create("Name", GDAL.OFTString)
GDAL.setwidth(fielddefn, 32)
@test GDAL.createfield(layer, fielddefn, GDAL.TRUE) == GDAL.OGRERR_NONE
GDAL.destroy(fielddefn)

# repeat to add multiple features
featuredefn = GDAL.getlayerdefn(layer)
feature = GDAL.f_create(featuredefn)
GDAL.setfieldstring(feature, GDAL.getfieldindex(feature, "Name"), "myname")
point = GDAL.creategeometry(GDAL.wkbPoint)
GDAL.setpoint_2d(point, 0, 100.123, 0.123)
@test GDAL.setgeometry(feature, point) == GDAL.OGRERR_NONE
GDAL.destroygeometry(point)

# test getfilelist with unsafe_loadstringlist
shp_exts = Set([".shp", ".shx", ".dbf"])
fileset = map(x -> pointshapefile * x, shp_exts)
@test Set(GDAL.getfilelist(dataset)) == fileset

GDAL.close(dataset)

map(rm, fileset)


text:
using VLKeggSDK
using ProgressMeter
using DataFrames

# get matching table -
metabolite_matching_table = build_metabolite_keggid_matching_table() |> check

# compute the master compound table -
# master_compound_array = Array{KEGGCompound,1}()
# value_itr = values(metabolite_matching_table)
# test_codes = value_itr |> collect
# @showprogress for kegg_compound_code in test_codes[1:10]

#     # mod_code -
#     mod_compund_code = "cpd:$(kegg_compound_code)"

#     # grab from cloud -
#     kegg_compound_obj = get_compound_records(mod_compund_code) |> check
#     if (isnothing(kegg_compound_obj) == false)
#         push!(master_compound_array, kegg_compound_obj)
#     end
# end


text:
using SQLite_Octo

adapter = SQL.connect(db="test.db")

SQL.execute(adapter, """
DROP TABLE IF EXISTS users; """)

SQL.execute(adapter, """
CREATE TABLE users (
    id INT NOT NULL,
    name VARCHAR(20),
    PRIMARY KEY (id)
); """)

SQL.execute(adapter, "SELECT name FROM sqlite_master WHERE type='table';")

SQL.execute(adapter, "INSERT INTO users (id, name) VALUES (1, 'John');")

SQL.execute(adapter, "SELECT * FROM users;")


text:
using Printf

export AbstractOptiOutput, OptiOutput

###########################################################################
# Status
###########################################################################
const STATUSES = Dict(
        :exception      => "unhandled exception",
        :first_order    => "first-order stationary",
        :acceptable     => "solved to within acceptable tolerances",
        :infeasible     => "problem may be infeasible",
        :max_eval       => "maximum number of function evaluations",
        :max_iter       => "maximum iteration",
        :max_time       => "maximum elapsed time",
        :neg_pred       => "negative predicted reduction",
        :not_desc       => "not a descent direction",
        :small_residual => "small residual",
        :small_step     => "step too small",
        :stalled        => "stalled",
        :unbounded      => "objective function may be unbounded from below",
        :unknown        => "unknown",
        :ieee_nan_inf   => "NaN or Inf occurred",
        :user           => "user-requested stop",
       )

###########################################################################
# Output
###########################################################################

abstract type AbstractOptiOutput end

"""
    OptiOutput(status; ...)
An OptiOutput is a struct for storing output information of solvers.
It contains the following fields:
- `status::Symbol`: Indicates the output of the solver. Use `show_statuses()` for the full list;
- `x::Vector`: The final approximation returned by the solver (default: `nothing`);
- `y::Vector`: The Lagrange multiplers wrt to the constraints (default: `nothing`);
- `objective::Real`: The objective value at `x` (default: `nothing`);
- `optimality::Real`: The dual feasibility norm at `x` (and `y`) (default: `nothing`);
- `cviolation::Real`: The primal feasibility norm at `x` (default: `nothing`);
- `iterations::Int`: The number of iterations used by the solver (default: `nothing`);
- `time::Real`: The time elapsed while running the solver, in seconds [s] (default: `nothing`);
- `solver_name::String`: The solver denomination (default: `nothing`);
- `solver::Dict{Symbol,Any}`: A solver specific dictionary.

The `status` is mandatory on construction. All other variables can be input as keyword arguments.

Notice that `OptiOutput` does not compute anything, it simply stores.
"""

mutable struct OptiOutput <: AbstractOptiOutput
    status::Symbol             # solver final status
    x::Maybe{Vector}           # primal solution
    y::Maybe{Vector}           # dual solution
    objective::Maybe{Real}     # value of objective function
    optimality::Maybe{Real}    # value of dual residual at (x,y) : optimality
    cviolation::Maybe{Real}    # value of primal residual at x : constraint violation
    iterations::Maybe{Int}     # iterations
    time::Maybe{Real}          # elapsed time
    solver_name::Maybe{String} # solver name
    solver::Dict{Symbol,Any}   # solver info
end

function OptiOutput( status::Symbol;
                     x::Maybe{Vector}=nothing,
                     y::Maybe{Vector}=nothing,
                     objective::Maybe{Real}=nothing,
                     optimality::Maybe{Real}=nothing,
                     cviolation::Maybe{Real}=nothing,
                     iterations::Maybe{Int}=nothing,
                     time::Maybe{Real}=nothing,
                     solver_name::Maybe{String}=nothing,
                     solver::Dict{Symbol,T}=Dict{Symbol,Any}() ) where {T}
    if !(status in keys(STATUSES))
        @error "Invalid status $status. Use one of the following: " join(keys(STATUSES), ", ")
        throw(KeyError(status))
    end
    return OptiOutput( status, x, y, objective, optimality, cviolation, iterations, time, solver_name, solver )
end

###########################################################################
# Printing
###########################################################################
import Base.show, Base.print, Base.println

function show(io::IO, out::AbstractOptiOutput)
    show(io, "Output: $(getStatus(out))")
end

function disp_vector(io::IO, x::Vector)
    if length(x) == 0
        print(io, "∅")
    elseif length(x) <= 5
        Base.show_delim_array(io, x, "[", " ", "]", false)
    else
        Base.show_delim_array(io, x[1:4], "[", " ", "", false)
        print(io, " ⋯ $(x[end])]")
    end
end

function print(io::IO, out::OptiOutput; showvec::Function=disp_vector)
    println(io, "Output")
    println(io, "  status: " * getStatus(out))
    if out.objective === nothing
        println(io, "  objective value: n/a ")
    else
        println(io, "  objective value: ", out.objective)
    end
    if out.cviolation === nothing
        println(io, "  constr. violation: n/a ")
    else
        println(io, "  constr. violation: ", out.cviolation)
    end
    if out.optimality === nothing
        println(io, "  optimality: n/a ")
    else
        println(io, "  optimality: ", out.optimality)
    end
    if out.x === nothing
        println(io, "  primal solution: n/a ")
    else
        print(io, "  primal solution: "); showvec(io, out.x); println(io, "")
    end
    if out.y === nothing
        println(io, "  dual solution: n/a ")
    else
        print(io, "  dual solution: "); showvec(io, out.y); println(io, "")
    end
    if out.iterations === nothing
        println(io, "  iterations: n/a ")
    else
        println(io, "  iterations: ", out.iterations)
    end
    if out.time === nothing
        println(io, "  elapsed time: n/a ")
    else
        println(io, "  elapsed time: ", out.time, " s")
    end
    if out.solver_name === nothing
        println(io, "  solver name: n/a ")
    else
        println(io, "  solver name: ", out.solver_name)
    end
    if length( out.solver ) > 0
        println(io, "  solver info:")
        for (k,v) in out.solver
            @printf(io, "    %s: ", k)
            if v isa Vector
                showvec(io, v)
            else
                show(io, v)
            end
            println(io, "")
        end
    end
end

print(out::OptiOutput; showvec::Function=disp_vector) = print(Base.stdout, out, showvec=showvec)

println(out::OptiOutput; showvec::Function=disp_vector) = print(Base.stdout, out, showvec=showvec)

println(io::IO, out::OptiOutput; showvec::Function=disp_vector) = print(io, out, showvec=showvec)

function getStatus(out::AbstractOptiOutput)
    return STATUSES[out.status]
end


text: module PredictMDTestModuleB end


text:
# [[file:../fresa.org::modulestart][modulestart]]
# All code copyright © Eric S Fraga.
# Date of last change in version variable below.
module Fresa
version = "[2021-07-03 13:58]"
using Dates
using Distributed
using Printf
function __init__()
    if myid() == 1
        println("# -*- mode: org; eval: (org-content 3); -*-")
        println(": Fresa PPA last change $version")
    end
end
# modulestart ends here

# [[file:../fresa.org::pointtype][pointtype]]
"""

Point (`x`) in the search space along with objective function values
(`z[]`) and feasbility indication (`g`).  The type of `x` is problem
specific.  `z[]` and `g` hold `Float64` values.  `g` should be of
length 1.

"""
struct Point
    x :: Any                    # decision point
    z :: Vector                 # objective function values
    g :: Float64                # constraint violation
    ancestor                    # the parent of this point
end
# pointtype ends here

# [[file:../fresa.org::showpoint][showpoint]]
import Base
Base.show(io::IO, p::Fresa.Point) = print(io, "f(", p.x, ")=", p.z, " g=", p.g)
# and also an array of points
function Base.show(io::IO, p::Array{Point,1})
    np = length(p)
    if np > 0
        nz = length(p[1].z)
        println(io, "|-")
        for i=1:nz
            print(io,"| z$(i) ")
        end
        println(io, "| g | x |")
        println(io,"|-")
        for i=1:length(p)
            for j=1:nz
                print(io,"| ", p[i].z[j], " ")
            end
            print(io, "| ", p[i].g, " ")
            print(io, "| ", p[i].x, " |\n")
        end
        println(io,"|-")
    else
        print(io,"empty")
    end
end
# showpoint ends here

# [[file:../fresa.org::pointsize][pointsize]]
import Base.size
Base.size(p :: Point) = ()
# pointsize ends here

# [[file:../fresa.org::ancestortype][ancestortype]]
struct Ancestor
    point :: Point        # the actual ancestor point
    fitness :: Float64    # the fitness of the ancestor
    generation :: Int32   # the generation when this point was created
end
# ancestortype ends here

# [[file:../fresa.org::*Ancestor <<ancestor>>][Ancestor <<ancestor>>:2]]
ancestor(p :: Point) = p.ancestor :: Union{Ancestor,Nothing}
# Ancestor <<ancestor>>:2 ends here

# [[file:../fresa.org::*Domain <<domain>>][Domain <<domain>>:1]]
struct Domain
    lower                       # function which returns lower bound on search variable(s)
    upper                       # function which returns upper bound on search variable(s)
end
# Domain <<domain>>:1 ends here

# [[file:../fresa.org::createpoint][createpoint]]
function createpoint(x,f,parameters = nothing,ancestor = nothing)
    z = 0
    g = 0
    if typeof(parameters) != Nothing
        (z, g) = f(x, parameters)
    else
        (z, g) = f(x)
    end
    if typeof(g) == Int
        g = Float64(g)
    end
    p = Nothing
    if rank(z) == 1
        p = Point(x, z, g, ancestor)
    elseif rank(z) == 0
        p = Point(x, [z], g, ancestor)
    else
        error("Fresa can only handle scalar and vector criteria, not $(typeof(z)).")
    end
    return p
end
# createpoint ends here

# [[file:../fresa.org::fitness][fitness]]
function fitness(pop, fitnesstype, steepness, generation, ngen)
    l = length(pop)
    indexfeasible = (1:l)[map(p->p.g,pop) .<= 0]
    indexinfeasible = (1:l)[map(p->p.g,pop) .> 0]
    @debug "Feasible/infeasible breakdown" indexfeasible indexinfeasible maxlog=3
    fit = zeros(l)
    factor = 1              # for placement in fitness interval (0,1)
    if length(indexfeasible) > 0
        feasible = view(pop,indexfeasible)
        # use objective function value(s) for ranking
        feasiblefit = vectorfitness(map(p->p.z,feasible), fitnesstype, steepness, generation, ngen)
        if length(indexinfeasible) > 0
            feasiblefit = feasiblefit./2 .+ 0.5 # upper half of fitness interval
            factor = 2                        # have both feasible & infeasible
        end
        fit[indexfeasible] = (feasiblefit.+factor.-1)./factor
    end
    if length(indexinfeasible) > 0
        # squeeze infeasible fitness values into (0,0.5) or (0,1) depending
        # on factor, i.e. whether there are any feasible solutions as well or not
        infeasible = view(pop,indexinfeasible)
        # use constraint violation for ranking as objective function values
        # may not make any sense given that points are infeasible
        fit[indexinfeasible] = vectorfitness(map(p->p.g, infeasible),
                                             fitnesstype,
                                             steepness,
                                             generation,
                                             ngen
                                             ) / factor;
    end
    fit
end
# fitness ends here

# [[file:../fresa.org::vectorfitness][vectorfitness]]
"""
For single objective problems, the fitness is simply the normalised
objective function value.

For multi-objective cases, there are three alternative measures of
fitness ranking possible.  The first is based on the Hadamard product
of the rank of each member of population accoring to each
criterion.  The second is based on a weighted Borda ranking based on
each criterion ranking.  Finally, a measure based on dominance,
similar to that used by the popular NSGA-II genetic algorithm, is
available.

"""
function vectorfitness(v, fitnesstype, steepness, generation, ngen)
    # determine number of objectives (or pseudo-objectives) to consider in
    # ranking
    l = length(v)
    if l == 1
        # no point in doing much as there is only one solution
        fit = [0.5]
    else
        m = length(v[1])
        # println("VF: v=$v")
        # println("  : of size $(size(v))")
        if m == 1                   # single objective
            fitness = [v[i][1] for i=1:l]
        else                  # multi-objective
            rank = ones(m,l); #rank of each solution for each objective function
            if fitnesstype == :hadamard
                for i=1:m
                    rank[i,sortperm([v[j][i] for j=1:l])] = 1:l
                end
                # hadamard product of ranks
                fitness = map(x->prod(x), rank[:,i] for i=1:l)
            elseif fitnesstype == :borda
                for i=1:m
                    rank[i,sortperm([v[j][i] for j=1:l])] = 1:l
                end
                # borda sum of ranks
                fitness = map(x->sum(x), rank[:,i] for i=1:l)
            elseif fitnesstype == :nondominated
                # similar to that used by NSGA-II (Deb 2000)
                fitness = zeros(l)
                maxl = assigndominancefitness!(fitness,v,1)
                # println("Resulting fitness: $fitness")
            else
                throw(ArgumentError("Type of fitness evaluation must be either :borda, :nondominated, or :hadamard, not $(repr(fitnesstype))."))
            end
        end
        # normalise (1=best, 0=worst) while avoiding
        # extreme 0,1 values using the hyperbolic tangent
        fit = adjustfitness(fitness, steepness, generation, ngen)
        # println(":  scaled fitness: $fit")
        @debug "Fitness calculations" v[1][1] v[2][1] v[l][1] fitness[1] fitness[2] fitness[l] fit[1] fit[2] fit[l] maxlog=3
    end
    fit
end
# vectorfitness ends here

# [[file:../fresa.org::adjustfitness][adjustfitness]]
function adjustfitness(fitness, steepness, generation, ngen)
    if (maximum(fitness)-minimum(fitness)) > eps()
        s = steepness
        if isa(steepness, Tuple)
            a = (2*steepness[1]-2*steepness[2])/3
            b = - (3*steepness[1] - 3*steepness[2])/ngen^2
            d = steepness[1]
            s = a*generation^3 + b*generation^2 + c*generation + d
            @debug "Steepness " s "at generation" g
        end
        fit = 0.5*(tanh.(4*s*(maximum(fitness) .- fitness)
                         / (maximum(fitness)-minimum(fitness))
                         .- 2*s) .+ 1)
    else
        # only one solution (or all solutions the same) in population
        fit = 0.5*ones(length(fitness))
    end
    fit
end
# adjustfitness ends here

# [[file:../fresa.org::assigndominancefitness][assigndominancefitness]]
function assigndominancefitness!(f,v,l)
    # assign value l to all members of v which dominate rest and then
    # recurse on those which are dominated
    (p, d) = paretoindices(v)
    # println("Assigning fitness $l to $p")
    f[p] .= l
    if !isempty(d)
        assigndominancefitness!(view(f,d),v[d],l+1)
    else
        l
    end
end
# assigndominancefitness ends here

# [[file:../fresa.org::neighbourarray][neighbourarray]]
function neighbour(x :: Array{Float64,1},
                   a :: Array{Float64,1},
                   b :: Array{Float64,1},
                   f :: Float64
                   ) :: Array{Float64,1}
    xnew = x .+ (1.0 .- f) .* 2(rand(length(x)).-0.5) .* (b.-a)
    xnew[xnew.<a] = a[xnew.<a];
    xnew[xnew.>b] = b[xnew.>b];
    return xnew
end
# neighbourarray ends here

# [[file:../fresa.org::neighbourfloat][neighbourfloat]]
function neighbour(x :: Float64,
                   a :: Float64,
                   b :: Float64,
                   f :: Float64
                   ) :: Float64
    # allow movements both up and down
    # in the domain for this variable
    newx = x + (b-a)*(2*rand()-1)/2.0 * (1-f)
    if newx < a
        newx = a
    elseif newx > b
        newx = b
    end
    newx
end
# neighbourfloat ends here

# [[file:../fresa.org::dominates][dominates]]
function dominates(a, b)
    all(a .<= b) && any(a .< b)
end
≻(a,b) = dominates(a,b)
# dominates ends here

# [[file:../fresa.org::*find Pareto set][find Pareto set:1]]
function paretoindices(z)
    n = length(z)
    dominance = [reduce(&, [!(z[i] ≻ z[j]) for i ∈ 1:n]) for j ∈ 1:n]
    paretoindices = filter(j -> dominance[j], 1:n)
    dominatedindices = filter(j -> !dominance[j], 1:n)
    (paretoindices, dominatedindices)
end
# find Pareto set:1 ends here

# [[file:../fresa.org::pareto][pareto]]
# indices of non-dominated and dominated points from the population of
# Point objects
function pareto(pop :: Vector{Point})
    l = length(pop)
    indexfeasible = (1:l)[map(p->p.g,pop) .<= 0]
    indexinfeasible = (1:l)[map(p->p.g,pop) .> 0]
    if length(indexfeasible) > 0
        subset = view(pop,indexfeasible)
        indices = indexfeasible
    else
        println(": Fresa.pareto warning: no feasible solutions.  Pareto set meaningless?")
        subset = pop
        indices = 1:l
    end
    z = map(p->p.z, subset)
    # use function below to return indices of non-dominated and
    # dominated from objective function values alone in the subset of
    # feasible solutions
    (p, d) = paretoindices(z)
    (indices[p], indices[d])
end
# pareto ends here

# [[file:../fresa.org::printhistorytrace][printhistorytrace]]
function printHistoryTrace(p :: Point)
    a = p.ancestor
    while typeof(a) != Nothing
        println("| $(a.generation) | $(a.fitness) |")
        a = a.point.ancestor
    end
end
# printhistorytrace ends here

# [[file:../fresa.org::prune][prune]]
function prune(pop :: AbstractArray, tolerance)
    npruned = 0
    z = map(p->p.z, pop)
    @show z[1]
    println("typeof(z)=$(typeof(z))")
    l = length(z)
    println("typeof(z[1])=$(typeof(z[1]))")
    n = length(z[1])
    @show n
    zmin = zeros(n)
    zmax = zeros(n)
    try
        for i=1:n
            row = [z[j][i] for j=1:l]
            zmin[i] = minimum(row)
            zmax[i] = maximum(row)
            if zmax[i] - zmin[i] < 100*eps()
                zmax[i] = zmin[i]+100*eps()
            end
        end
        pruned = [pop[1]]
        for i=2:l
            similar = false
            for j=1:length(pruned)
                if all(abs.(z[i]-pruned[j].z) .< tolerance*(zmax-zmin))
                    similar = true;
                    break;
                end
            end
            if !similar
                push!(pruned,pop[i])
            else
                npruned += 1
            end
        end
        (pruned, npruned)
    catch e
        println("prune method error: $e")
        if isa(e, MethodError)
            # probably (possibly) due to objective function type not
            # being a number.  In this case, we try again but looking
            # at the decision variable values instead.
            x = map(p->p.x, pop)
            # println("typeof(z)=$(typeof(z))")
            l = length(x)
            # start building up the population that remains after
            # pruning.  The first entry will always be there as any
            # similar solutions will not be included by the search
            # that follows.
            pruned = [pop[1]]
            try
                for i=2:l
                    similar = false
                    # now check this solution against all those already in
                    # the list we are collating
                    for j=1:length(pruned)
                        if all(Float64.(abs.(x[i]-pruned[j].x)) .< tolerance)
                            similar = true;
                            break;
                        end
                    end
                    if !similar
                        push!(pruned,pop[i])
                    else
                        npruned += 1
                    end
                end
                (pruned, npruned)
            catch e
                println("prune method second error: $e")
                if isa(e, MethodError)
                    # this is now probably/possibly due to not being
                    # to find the difference between two decision
                    # points.  In that case, return the whole
                    # original population
                    (pop, 0)
                else
                    @error "Unexpected error in prune method for points" e
                end             # if method error
            end                 # try pruning on points
        else
            @error "Unexpected error in prune method for objective function values" e
        end                     # if is a method error
    end                         # try pruning on objective function values
end                             # function
# prune ends here

# [[file:../fresa.org::randompopulation][randompopulation]]
function randompopulation(n,f,parameters,a,b)
    p = Point[]                 # population object
    for j in 1:n
        push!(p, createpoint(randompoint(a,b), f, parameters))
    end
    p
end
# randompopulation ends here

# [[file:../fresa.org::randompoint][randompoint]]
function randompoint(a,b)
    x = a + rand(length(a)).*b
end
# randompoint ends here

# [[file:../fresa.org::select][select]]
function select(f)
    l = length(f)
    ind1 = rand(1:l)
    if ind1 == 0
        ind1 = 1
    end
    ind2 = rand(1:l)
    # println("Comparing $ind1 to $ind2")
    if f[ind1] > f[ind2]
        return ind1
    else
        return ind2
    end
end
# select ends here

# [[file:../fresa.org::solve][solve]]
"""

Solve an optimisation problem, defined as the minimization of the
values returned by the objective function, `f`.  `f` returns not only
the objective function values, an array of `Float64` values, but also
a measure of feasibility (≤0) or infeasibility (>0).  The problem is
solved using the Fresa algorithm.  `p0` is the initial population
which has to have at least one member, a `Point`, and `a` and `b` are
*bounds* on the search space.

The return values for the solution of a single criterion problem are
the best point and the full population at the end of the search.

For a multi-objective problem, the returned values are the set of
indices for the points within the full population (the second returned
value) approximating the *Pareto* front.

The population will consist of an array of `Fresa.Point` objects, each
of which will have the point in the search space, the objective
function value and the feasibility measure.

"""
function solve(f, p0, domain;        # required arguments
               parameters = nothing, # allow parameters for objective function
               archiveelite = false, # save thinned out elite members
               elite = true,    # elitism by default
               fitnesstype = :hadamard, # how to rank solutions in multi-objective case
               multithreading = false, # use multiple threads for objective function evaluation
               ngen = 100,           # number of generations
               npop = 10,            # population size: fixed (single value) or dynamic (tuple)
               nrmax = 5,            # number of runners maximum
               ns = 100,             # number of stable solutions for stopping
               output = 1,           # how often to output information
               plotvectors = false,  # generate output file for search plot
               populationoutput = false, # output population every generation?
               steepness = 1.0,      # show steep is the adjustment shape for fitness
               tolerance = 0.001,    # tolerance for similarity detection
               usemultiproc = false) # parallel processing by Fresa itself?
    output > 0 && println("** solve $f $(orgtimestamp(now()))")
    tstart = time()
    nf = 1                   # number of function evaluations
    npruned = 0              # number solutions pruned from population
    nz = length(p0[1].z)     # number of criteria
    pop = copy(p0);          # create/initialise the population object
    if archiveelite
        archive = Point[]
    end
    if output > 0
        println("#+name: $(f)settings")
        println("| variable | value |")
        println("|-")
        println("| ngen | $ngen |")
        println("| npop | $npop |")
        println("| nrmax | $nrmax |")
        println("| ns | $ns |")
        println("| elite | $elite |")
        println("| archive | $archiveelite |")
        println("| fitness | $fitnesstype |")
        println("| steepness | $steepness |")
        println("|-")
        # output != 0 && println(": solving with ngen=$ngen npop=$npop nrmax=$nrmax ns=$ns")
        # output != 0 && println(": elite=$elite archive elite=$archiveelite fitness type=$fitnesstype")
    end
    if plotvectors
        plotvectorio = open("fresa-vectors-$(orgtimestamp(now())).data", create=true, write=true)
        output > 0 && println(": output of vectors for subsequent plotting")
    end
    # if npop was given as a tuple, we are to have a dynamic
    # population size.  This only makes sense for multi-objective
    # optimization problems so a warning will be given otherwise.
    npopmin = npop
    npopmax = npop
    if isa(npop, Tuple)
        if nz > 1
            npopmin = npop[1]
            npopmax = npop[2]
            if npopmin > npopmax
                error("Dynamic population sizing requires min <= max; you specified $npop")
            end
            npop = npopmin      # start with minimum possible
        else
            println("*Warning*: you have specified a tuple for population size: $npop")
            println("This only makes sense for multi-objective optimization problems.")
            println("npop will be set to $(npop[1]).")
            npop = npop[1]      # be optimistic and use minimum given
        end
    end
    # we use multithreading if asked for *and* if we have more than
    # one thread available
    multithreading = multithreading && Threads.nthreads() > 1
    # we use parallel computing if we have more than one processor
    parallel = usemultiproc && nprocs() > 1
    # parallel = false
    if output > 0
        println(": function evaluations performed ",
                parallel
                ? "in parallel with $(nprocs()) processors."
                : (multithreading
                   ? "in parallel with $(Threads.nthreads()) threads."
                   : "sequentially."))
        println("*** initial population")
        println("#+name: $(f)initial")
        println(pop)
    end
    if output > 0
        println("*** evolution")
        println("#+name: $(f)evolution")
        println("#+plot: ind:1 deps:(6) with:\"points pt 7\" set:\"logscale x\"")
        @printf("| %9s | %9s | %9s | %9s | %9s |", "gen", "npop",
                (elite && nz > 1) ? "pareto" : "nf", "pruned", "t (s)")
        for i in 1:nz
            @printf(" z%-8d |", i)
        end
        @printf(" %9s |", "g")
        @printf("\n|-\n")
    end
    # now evolve the population for a predetermined number of generations
    for gen in 1:ngen
        # evaluate fitness which is adjusted depending on value of
        # steepness, a value that may depend on the generation
        fit = fitness(pop, fitnesstype, steepness, gen, ngen)
        if gen == 1
            @debug "Initial fitness" f=fit
        end
        # sort
        index = sortperm(fit)
        if populationoutput
            println("\nGeneration $gen full population is:")
            println(pop)
            println("Fitness vector: $fit")
        end
        # and remember best which really only makes sense in single
        # criterion problems but is best in multi-objective case in
        # the ranking measure used by Fresa
        best = pop[index[end]]
        # if elitism is used
        if elite
            if nz > 1
                # elite set is whole pareto set unless it is too
                # big. Recall that the pareto function returns the set
                # of indices into the population
                wholepareto = pareto(pop)[1]
                # if using dynamic population sizing, adjust the population
                npop = 2 * length(wholepareto)
                if npop < npopmin
                    npop = npopmin
                end
                if npop > npopmax
                    npop = npopmax
                end
                # now check that the pareto is not too big.  if it is, thin it out
                if length(wholepareto) > ceil(npop/2)
                    newpop, removed = thinout(pop, fit, wholepareto, ceil(Int,npop/2))
                    if archiveelite
                        archive = prune(append!(archive, removed), tolerance)[1]
                        archive = archive[pareto(archive)[1]]
                    end
                else
                    newpop = pop[wholepareto]
                end
            else
                # elite set is single element only
                newpop = [best]
            end
            # if plotting vectors for the search, include elitism
            if plotvectors
                for p in newpop
                    write(plotvectorio, "$(gen-1) $(p.x)\n$gen $(p.x)\n\n")
                end
            end
        else
            newpop = Point[]
        end
        if output >= 0
            print(stderr, ": $gen np=$(length(newpop))/$npop",
                  archiveelite ? " na=$(length(archive))" : "",
              " with most fit z=$(best.z)           \r")
            # if output has been requested, check to see if output is
            # required now and then also check to see if the frequency
            # needs to be reduced.
            if output > 0
                if gen%output == 0
                    @printf("| %9d | %9d | %9d | %9d | %9.2f |", gen, length(fit),
                            (elite && nz > 1) ? length(newpop) : nf, npruned, time()-tstart)
                    for i = 1:length(best.z)
                        print(" $(best.z[i]) |")
                    end
                    print(" $(best.g) |")
                    println()
                end
                if 10^(floor(log10(gen))) > output
                    output = 10^(Int(floor(log10(gen))))
                end
            end
        end
        # if we are using any form of multiprocessing, either threads
        # or multiple cores, create an array to store all new points
        # which we evaluate later in parallel.  Ideally, also keep
        # track of the points from which new points are derived to
        # provide the backward link through the evolution but this is
        # currently disabled as the creation of the Ancestor object
        # requires more information than I am currently storing away.
        if multithreading || parallel
            x = Any[] # typeof(newpop[1].x)[]
            # points = Point[]
        end
        # now loop through population, applying selection and then
        # generating neighbours
        l = length(pop)
        for i in 1:min(l,npop)
            s = select(fit)
            # println(": selection $i is $s")
            # println(": size of pop is $(size(pop))")
            selected = pop[s]
            if !elite
                # if no elitism, we ensure selected members remain in population
                push!(newpop, selected)
                if plotvectors
                    write(plotvectorio, "$(gen-1) $(selected.x)\n$gen $(selected.x)\n\n")
                end
            end
            # number of runners to generate, function of fitness
            nr = ceil(fit[s]*nrmax*rand())
            if nr < 1
                nr = 1
            end
            # println(": generating $nr runners")
            for r in 1:nr
                # create a neighbour, also function of fitness,
                # passing through the lower and upper bounds
                # appropriate for the particular solution point
                newx = neighbour(pop[s].x, domain.lower(pop[s].x), domain.upper(pop[s].x), fit[s])
                # for parallel evaluation, we store the neighbours and
                # evaluate them later; otherwise, we evaluate
                # immediately and save the resulting point
                if multithreading || parallel
                    push!(x, newx)
                    # push!(points, pop[s])
                else
                    push!(newpop, createpoint(newx, f, parameters, Ancestor(pop[s],fit[s],gen)))
                    if plotvectors
                        write(plotvectorio, "$(gen-1) $(pop[s].x)\n$gen $newx\n\n")
                    end
                    nf += 1
                end
            end
            # remove selected member from the original population so
            # it is not selected again
            splice!(fit, s)
            splice!(pop, s)
        end
        # if we are making use of parallel computing, we evaluate all
        # points generated in previous loop.  Parallel processing is
        # done either via multithreading or with multiple
        # processors.  The former is easier as it's based on shared
        # memory.
        if multithreading       # using threads and shared memory
            results = Array{Point}(undef,length(x))
            Threads.@threads for i ∈ 1:length(x)
                results[i] = createpoint(x[i],f,parameters)
            end
            append!(newpop, results)
        elseif parallel        # using multiple processors with remote calls
            # will be used to collect results from worker processors
            results = Array{Future,1}(undef, nprocs())
            i = 0;
            while i < length(x)
                # issue remote evaluation call
                for j=1:nprocs()
                    if i+j <= length(x)
                        # TODO: the information about the ancestor is
                        # not available; this needs to be stored above
                        results[j] = @spawn createpoint(x[i+j],f,parameters)
                        nf += 1
                    end
                end
                # now wait for results
                for j=1:nprocs()
                    if i+j <= length(x)
                        push!(newpop, fetch(results[j]))
                    end
                end
                i += nprocs()
            end
        end
        # and finally, if we have elitism, remove any duplicate points
        # in the new population and make it the current population for
        # the next generation; otherwise, simply copy over
        if elite
            (pop, nn) = prune(newpop, tolerance)
            npruned += nn
        else
            pop = newpop
        end
    end
    output > 0 && println("*** Fresa run finished\n: nf=$nf npruned=$npruned", archiveelite ? " archived=$(length(archive))" : "")
    if plotvectors
        close(plotvectorio)
    end
    if nz == 1
        fit = fitness(pop, fitnesstype, steepness, ngen, ngen)
        index = sortperm(fit)
        best = pop[index[end]]
        return best, pop
    else
        return pareto(archiveelite ? append!(pop,archive) : pop)[1], pop
    end
end
# solve ends here

# [[file:../fresa.org::thinout][thinout]]
function thinout(pop, fit, pareto, n::Int)
    indices = sortperm(fit[pareto])
    return pop[pareto[indices[end-n+1:end]]], pop[pareto[indices[1:end-n]]]
end
# thinout ends here

# [[file:../fresa.org::orgtimestamp][orgtimestamp]]
function orgtimestamp(dt::DateTime)
    return @sprintf("[%d-%02d-%02d %02d:%02d]",
                    Dates.year(dt),
                    Dates.month(dt),
                    Dates.day(dt),
                    Dates.hour(dt),
                    Dates.minute(dt))
end
# orgtimestamp ends here

# [[file:../fresa.org::rank][rank]]
rank(x :: Any) = length(size(x))
# rank ends here

# [[file:../fresa.org::moduleend][moduleend]]
end
# moduleend ends here


text:
using Test, CompilerBridge

@testset "CompilerBridge.jl" begin
    include("test_bridges.jl")
end

max_stars_repo_name: aviatesk/CompilerBridge

max_issues_repo_name: aviatesk/CompilerBridge

max_forks_repo_name: aviatesk/CompilerBridge


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   Tests of formatters.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

@testset "Formatters" begin
    ft = (v,i,j)->v === missing ? "\$\emptyset\$" : v

    matrix = deepcopy(data)
    matrix[3,3] = missing

    expected = """
\\begin{table}
  \\begin{tabular}{rrrr}
    \\hline\\hline
    \\textbf{Col. 1} & \\textbf{Col. 2} & \\textbf{Col. 3} & \\textbf{Col. 4} \\\\\\hline
    1 & false & 1.0 & 1 \\\\
    2 & true & 2.0 & 2 \\\\
    3 & false & \$\\emptyset\$ & 3 \\\\
    4 & true & 4.0 & 4 \\\\
    5 & false & 5.0 & 5 \\\\
    6 & true & 6.0 & 6 \\\\\\hline\\hline
  \\end{tabular}
\\end{table}
"""

    result = pretty_table(String, matrix, backend = :latex, formatters = ft)
    @test result == expected
end

@testset "Pre-defined formatters" begin
    matrix = hcat(1:1:11,
                  1:1.0:11,
                  [10^i for i = -5:1.:5],
                  [ i == 5 ? nothing : "Teste" for i = 1:11])

    expected = """
\\begin{table}
  \\begin{tabular}{rrrr}
    \\hline\\hline
    \\textbf{Col. 1} & \\textbf{Col. 2} & \\textbf{Col. 3} & \\textbf{Col. 4} \\\\\\hline
    1 & 1 & \$1 \\cdot 10^{-5}\$ & Teste \\\\
    2 & 2 & 0.0001 & Teste \\\\
    3 & 3 & 0.001 & Teste \\\\
    4 & 4 & 0.01 & Teste \\\\
    5 & 5 & 0.1 & nothing \\\\
    6 & 6 & 1 & Teste \\\\
    7 & 7 & \$1 \\cdot 10^{1}\$ & Teste \\\\
    8 & 8 & \$1 \\cdot 10^{2}\$ & Teste \\\\
    9 & 9 & \$1 \\cdot 10^{3}\$ & Teste \\\\
    \$1 \\cdot 10^{1}\$ & \$1 \\cdot 10^{1}\$ & \$1 \\cdot 10^{4}\$ & Teste \\\\
    \$1 \\cdot 10^{1}\$ & \$1 \\cdot 10^{1}\$ & \$1 \\cdot 10^{5}\$ & Teste \\\\\\hline\\hline
  \\end{tabular}
\\end{table}
"""

    result = pretty_table(String, matrix, backend = :latex,
                          formatters = ft_latex_sn(1))
    @test result == expected

    expected = """
\\begin{table}
  \\begin{tabular}{rrrr}
    \\hline\\hline
    \\textbf{Col. 1} & \\textbf{Col. 2} & \\textbf{Col. 3} & \\textbf{Col. 4} \\\\\\hline
    1 & 1 & \$1 \\cdot 10^{-5}\$ & Teste \\\\
    2 & 2 & 0.0001 & Teste \\\\
    3 & 3 & 0.001 & Teste \\\\
    4 & 4 & 0.01 & Teste \\\\
    5 & 5 & 0.1 & nothing \\\\
    6 & 6 & 1 & Teste \\\\
    7 & 7 & 10 & Teste \\\\
    8 & 8 & \$1 \\cdot 10^{2}\$ & Teste \\\\
    9 & 9 & \$1 \\cdot 10^{3}\$ & Teste \\\\
    10 & 10 & \$1 \\cdot 10^{4}\$ & Teste \\\\
    11 & 11 & \$1 \\cdot 10^{5}\$ & Teste \\\\\\hline\\hline
  \\end{tabular}
\\end{table}
"""

    result = pretty_table(String, matrix, backend = :latex,
                          formatters = ft_latex_sn(2))
    @test result == expected

    expected = """
\\begin{table}
  \\begin{tabular}{rrrr}
    \\hline\\hline
    \\textbf{Col. 1} & \\textbf{Col. 2} & \\textbf{Col. 3} & \\textbf{Col. 4} \\\\\\hline
    1 & 1.0 & \$1 \\cdot 10^{-5}\$ & Teste \\\\
    2 & 2.0 & 0.0001 & Teste \\\\
    3 & 3.0 & 0.001 & Teste \\\\
    4 & 4.0 & 0.01 & Teste \\\\
    5 & 5.0 & 0.1 & nothing \\\\
    6 & 6.0 & 1 & Teste \\\\
    7 & 7.0 & \$1 \\cdot 10^{1}\$ & Teste \\\\
    8 & 8.0 & \$1 \\cdot 10^{2}\$ & Teste \\\\
    9 & 9.0 & \$1 \\cdot 10^{3}\$ & Teste \\\\
    \$1 \\cdot 10^{1}\$ & 10.0 & \$1 \\cdot 10^{4}\$ & Teste \\\\
    \$1 \\cdot 10^{1}\$ & 11.0 & \$1 \\cdot 10^{5}\$ & Teste \\\\\\hline\\hline
  \\end{tabular}
\\end{table}
"""

    result = pretty_table(String, matrix, backend = :latex,
                          formatters = ft_latex_sn(1, [1, 3]))
    @test result == expected

    expected = """
\\begin{table}
  \\begin{tabular}{rrrr}
    \\hline\\hline
    \\textbf{Col. 1} & \\textbf{Col. 2} & \\textbf{Col. 3} & \\textbf{Col. 4} \\\\\\hline
    1 & 1.0 & \$1 \\cdot 10^{-5}\$ & Teste \\\\
    2 & 2.0 & 0.0001 & Teste \\\\
    3 & 3.0 & 0.001 & Teste \\\\
    4 & 4.0 & 0.01 & Teste \\\\
    5 & 5.0 & 0.1 & nothing \\\\
    6 & 6.0 & 1 & Teste \\\\
    7 & 7.0 & 10 & Teste \\\\
    8 & 8.0 & 100 & Teste \\\\
    9 & 9.0 & \$1 \\cdot 10^{3}\$ & Teste \\\\
    \$1 \\cdot 10^{1}\$ & 10.0 & \$1 \\cdot 10^{4}\$ & Teste \\\\
    \$1 \\cdot 10^{1}\$ & 11.0 & \$1 \\cdot 10^{5}\$ & Teste \\\\\\hline\\hline
  \\end{tabular}
\\end{table}
"""

    result = pretty_table(String, matrix, backend = :latex,
                          formatters = ft_latex_sn([1, 3], [1, 3]))
    @test result == expected
end


text:
using ProtoGrad
using ProtoGrad: Model, Linear, Conv, ReLU, relu, softmax, Dropout, Compose, SupervisedObjective, mse
using LinearAlgebra
using Serialization
using Statistics
using Test

@testset "Basic operations" begin
    @testset "$(T)" for T in [Float16, Float32, Float64]
        @testset "$(name)" for (name, m) in [
            ("Linear", Linear(T, 3=>2)),
            ("Compose", Compose(Linear(T, 4=>3), ReLU(), Linear(T, 3=>2), x -> relu.(x)))
        ]
            vec_m = vec(m)

            @test eltype(vec_m) == T

            @test length(ProtoGrad.overlap(m, m)) == length(ProtoGrad.allparams(m))

            res = m + m
            @test typeof(res) == typeof(m)
            vec_res = vec(res)
            @test eltype(vec_res) == T
            @test vec_res == vec_m + vec_m

            res .= m .+ m
            @test typeof(res) == typeof(m)
            vec_res = vec(res)
            @test eltype(vec_res) == T
            @test vec_res == vec_m + vec_m

            res = m - m
            @test typeof(res) == typeof(m)
            vec_res = vec(res)
            @test eltype(vec_res) == T
            @test vec_res == vec_m - vec_m

            res .= m .- m
            @test typeof(res) == typeof(m)
            vec_res = vec(res)
            @test eltype(vec_res) == T
            @test vec_res == vec_m - vec_m

            res = 2 * m
            @test typeof(res) == typeof(m)
            vec_res = vec(res)
            @test eltype(vec_res) == T
            @test vec_res == 2 * vec_m

            res .= 2 .* m
            @test typeof(res) == typeof(m)
            vec_res = vec(res)
            @test eltype(vec_res) == T
            @test vec_res == 2 * vec_m

            res = m / 3
            @test typeof(res) == typeof(m)
            vec_res = vec(res)
            @test eltype(vec_res) == T
            @test vec_res == vec_m / 3

            res .= m ./ 3
            @test typeof(res) == typeof(m)
            vec_res = vec(res)
            @test eltype(vec_res) == T
            @test vec_res == vec_m / 3

            mp1 = m .+ 1
            dot_m_mp1 = dot(m, mp1)
            @test typeof(dot_m_mp1) == T
            @test dot_m_mp1 ≈ dot(vec_m, vec(mp1))

            mcopy = copy(m)
            @test typeof(mcopy) == typeof(m)
            vec_mcopy = vec(mcopy)
            @test eltype(vec_mcopy) == T
            @test vec_mcopy == vec_m
            @test length(ProtoGrad.overlap(m, mcopy)) == 0

            msimilar = similar(m)
            @test typeof(msimilar) == typeof(m)
            vec_msimilar = vec(msimilar)
            @test eltype(vec_msimilar) == T

            mzero = zero(m)
            @test typeof(mzero) == typeof(m)
            vec_mzero = vec(mzero)
            @test eltype(vec_mzero) == T
            @test all(vec_mzero .== T(0))

            msum = sum((k * m for k in 1:3))
            @test typeof(msum) == typeof(m)
            vec_msum = vec(msum)
            @test eltype(vec_msum) == T
            @test vec_msum == 6 * vec_m

            mmean = mean((k * m for k in [1, 2, 4]))
            @test typeof(mmean) == typeof(m)
            vec_mmean = vec(mmean)
            @test eltype(vec_mmean) == T
            @test vec_mmean ≈ (vec_m * 7) / 3

            coll_m = collect(m)
            @test_skip eltype(coll_m) == T
            @test coll_m == vec_m

            length_m = length(m)
            @test length_m == length(vec_m)

            v = m[end-1]
            @test vec_m[end-1] == v
            m[end-1] = T(42)
            @test vec(m)[end-1] == T(42)
        end
    end
end

@testset "Gradient" begin
    @testset "Linear ($(T))" for T in [Float32, Float64]
        input_size, output_size = 1000, 1
        batch_size = 50

        W = randn(T, output_size, input_size)
        b = randn(T, output_size)
        m = Linear(W, b)

        x = randn(T, input_size, batch_size)
        y = W * x .+ b + randn(T, output_size, batch_size)

        data_iter = Iterators.repeated((x, y))
        f = SupervisedObjective(mse, data_iter)

        grad, out = ProtoGrad.gradient(f, m)

        @test typeof(out) <: Number  # TODO improve this assertion
        @test typeof(grad) == typeof(m)

        @test out ≈ (norm(m.W * x .+ m.b - y) ^ 2) / batch_size
        @test grad.W ≈ (2 / batch_size) * ((m.W * x .+ m.b - y) * x')
        @test grad.b ≈ (2 / batch_size) * (m.W * x .+ m.b - y) * ones(batch_size)

        m .= m .- 3.0 .* grad
    end

    @testset "Compose ($(T))" for T in [Float32, Float64]
        input_size, hidden_size, output_size = 1000, 50, 1
        batch_size = 64

        W1 = randn(T, hidden_size, input_size)
        b1 = randn(T, hidden_size)
        m1 = Linear(W1, b1)

        W2 = randn(T, output_size, hidden_size)
        b2 = randn(T, output_size)
        m2 = Linear(W2, b2)

        m = Compose(m1, ReLU(), m2, x -> relu.(x))

        W_true = randn(T, output_size, input_size)
        b_true = randn(T, output_size)

        x = randn(T, input_size, batch_size)
        y = W_true * x .+ b_true + randn(T, output_size, batch_size)

        data_iter = Iterators.repeated((x, y))
        f = SupervisedObjective(mse, data_iter)

        grad, out = ProtoGrad.gradient(f, m)

        @test typeof(out) <: Number  # TODO improve this assertion
        @test typeof(grad) == typeof(m)

        m .= m .- 3.0 .* grad
    end

    @testset "Custom ($(T))" for T in [Float32, Float64]
        struct CustomModel <: Model
            m1
            act1
            m2
            act2
        end

        (m::CustomModel)(x) = (m.act2 ∘ m.m2 ∘ m.act1 ∘ m.m1)(x)

        input_size, hidden_size, output_size = 1000, 50, 1
        batch_size = 64

        W1 = randn(T, hidden_size, input_size)
        b1 = randn(T, hidden_size)
        m1 = Linear(W1, b1)

        W2 = randn(T, output_size, hidden_size)
        b2 = randn(T, output_size)
        m2 = Linear(W2, b2)

        m = CustomModel(m1, ReLU(), m2, x -> relu.(x))

        W_true = randn(T, output_size, input_size)
        b_true = randn(T, output_size)

        x = randn(T, input_size, batch_size)
        y = W_true * x .+ b_true + randn(T, output_size, batch_size)

        data_iter = Iterators.repeated((x, y))
        f = SupervisedObjective(mse, data_iter)

        grad, out = ProtoGrad.gradient(f, m)

        @test typeof(out) <: Number  # TODO improve this assertion
        @test typeof(grad) == typeof(m)

        m .= m .- 3.0 .* grad
    end

end

@testset "Dropout" begin
    for T in [Float32, Float64]
        input_size, hidden_size, output_size = 1000, 50, 1
        batch_size = 64

        m = ProtoGrad.Dropout(0.1)

        x = randn(T, input_size, batch_size)
        y = randn(T, output_size, batch_size)

        data_iter = Iterators.repeated((x, y))
        f = model -> ProtoGrad.mse(model(x), y)

        grad, out = ProtoGrad.within_training_mode() do
            ProtoGrad.gradient(f, m)
        end

        @test typeof(out) <: Number  # TODO improve this assertion
        @test typeof(grad) == typeof(m)

        m .= m .- 3.0 .* grad
    end
end

@testset "Serde" begin
    for (m, x) in [
        (Linear(Float32, 20=>10), randn(Float32, 20, 4)),
        (Linear(Float64, 20=>10), randn(Float64, 20, 4)),
        (Conv(Float32, 3=>5, (3, 3)), randn(Float32, 10, 10, 3, 4)),
        (Conv(Float64, 3=>5, (3, 3)), randn(Float64, 10, 10, 3, 4)),
        (
            Compose(
                Linear(20=>10),
                ReLU(),
                Linear(10=>5),
                x -> relu.(x),
                softmax
            ),
            randn(Float32, 20, 4)
        )
    ]
        out = m(x)

        (path, io) = mktemp()
        Serialization.serialize(io, m)
        close(io)

        io = open(path, "r")
        m_copy = Serialization.deserialize(io)
        close(io)

        # NOTE this is not guaranteed to work
        # e.g. when local definitions (such as functions) are serialized
        # @test typeof(m_copy) == typeof(m)

        out_copy = m_copy(x)

        @test out_copy == out
    end
end


text:
# Main module
module SimBioReactor
    # Function to call the GUI
    export SimBioReactorGUI

    # Include the main file .fl
    include("mainGUI.jl")
end


text:
function get!push!(d; kwargs...)
    for (k, val) in kwargs
        push!(get!(d, k, []), val)
    end
    d
end


text:
"""
Hartree term: for a decaying potential V the energy would be

1/2 ∫ρ(x)ρ(y)V(x-y) dxdy

with the integral on x in the unit cell and of y in the whole space.
For the Coulomb potential with periodic boundary conditions, this is rather

1/2 ∫ρ(x)ρ(y) G(x-y) dx dy

where G is the Green's function of the periodic Laplacian with zero
mean (-Δ G = sum_{R} 4π δ_R, integral of G zero on a unit cell).
"""
struct Hartree
    scaling_factor::Real  # to scale by an arbitrary factor (useful for exploration)
end
Hartree(; scaling_factor=1) = Hartree(scaling_factor)
(hartree::Hartree)(basis) = TermHartree(basis, hartree.scaling_factor)

struct TermHartree <: Term
    basis::PlaneWaveBasis
    scaling_factor::Real  # scaling factor, absorbed into poisson_green_coeffs
    # Fourier coefficients of the Green's function of the periodic Poisson equation
    poisson_green_coeffs
end
function TermHartree(basis::PlaneWaveBasis{T}, scaling_factor) where T
    # Solving the Poisson equation ΔV = -4π ρ in Fourier space
    # is multiplying elementwise by 4π / |G|^2.
    poisson_green_coeffs = 4T(π) ./ [sum(abs2, G)
                                     for G in G_vectors_cart(basis)]

    # Zero the DC component (i.e. assume a compensating charge background)
    poisson_green_coeffs[1] = 0
    TermHartree(basis, scaling_factor, scaling_factor .* poisson_green_coeffs)
end

@timing "ene_ops: hartree" function ene_ops(term::TermHartree, ψ, occ; ρ, kwargs...)
    basis = term.basis
    T = eltype(basis)
    pot_fourier = term.poisson_green_coeffs .* ρ.fourier
    potential = real(G_to_r(basis, pot_fourier))  # TODO optimize this
    E = real(dot(pot_fourier, ρ.fourier) / 2)

    ops = [RealSpaceMultiplication(basis, kpoint, potential) for kpoint in basis.kpoints]
    (E=E, ops=ops)
end

function compute_kernel(term::TermHartree; kwargs...)
    @assert term.basis.model.spin_polarization in (:none, :spinless, :collinear)
    vc_G = term.poisson_green_coeffs
    # Note that `real` here: if omitted, will result in high-frequency noise of even FFT grids
    K = real(G_to_r_matrix(term.basis) * Diagonal(vec(vc_G)) * r_to_G_matrix(term.basis))

    # Hartree kernel is independent of spin, so to apply it to (ρtot, ρspin)^T
    # and obtain the same contribution to Vα and Vβ the operator has the block structure
    #     ( K 0 )
    #     ( K 0 )
    n_spin = term.basis.model.n_spin_components
    n_spin == 1 ? K : [K 0I; K 0I]
end

function apply_kernel(term::TermHartree, dρ::RealFourierArray, dρspin=nothing; kwargs...)
    @assert term.basis.model.spin_polarization in (:none, :spinless, :collinear)
    kernel = from_fourier(dρ.basis, term.poisson_green_coeffs .* dρ.fourier)
    n_spin = term.basis.model.n_spin_components
    n_spin == 1 ? (kernel, ) : (kernel, kernel)  # Hartree term does not care about spin
end


text:
## Functions
# Pendulum equations of motion
function pendulum!(dx, x, p, t; M=0)
	θ, ω = x
	g, m, L, c, I = p

	dx[1] = ω
	dx[2] = (M + m*g*sin(θ)*L/2 - c*ω) / I
	return nothing
end
function pendulum!(D, vars::ComponentArray, params, t; M=0)
	@unpack θ, ω = vars
	@unpack g, m, L, c, I = params

	D.θ = ω
	D.ω = (M + m*g*sin(θ)*L/2 - c*ω) / I
	return nothing
end

pendulum_with_control(u) = (dx,x,p,t) -> pendulum!(dx,x,p,t; M=u(x,p.control,t))


text:
# Forward methods from generic domains to concrete interpreter #

satisfy(domain::GenericDomain, state::State, term::Term) =
    satisfy(ConcreteInterpreter(), domain, state, term)
satisfy(domain::GenericDomain, state::State, terms::AbstractVector{<:Term}) =
    satisfy(ConcreteInterpreter(), domain, state, terms)

satisfiers(domain::GenericDomain, state::State, term::Term) =
    satisfiers(ConcreteInterpreter(), domain, state, term)
satisfiers(domain::GenericDomain, state::State, terms::AbstractVector{<:Term}) =
    satisfiers(ConcreteInterpreter(), domain, state, terms)

evaluate(domain::GenericDomain, state::State, term::Term) =
    evaluate(ConcreteInterpreter(), domain, state, term)

initstate(domain::GenericDomain, problem::Problem) =
    initstate(ConcreteInterpreter(), domain, problem)
initstate(domain::GenericDomain, objtypes::AbstractDict) =
    initstate(ConcreteInterpreter(), domain, objtypes)
initstate(domain::GenericDomain, objtypes::AbstractDict, fluents) =
    initstate(ConcreteInterpreter(), domain, objtypes, fluents)

goalstate(domain::GenericDomain, problem::Problem) =
    goalstate(ConcreteInterpreter(), domain, problem)
goalstate(domain::GenericDomain, objtypes::AbstractDict, terms) =
    goalstate(ConcreteInterpreter(), domain, objtypes, terms)

transition(domain::GenericDomain, state::State, action::Term; options...) =
    transition(ConcreteInterpreter(), domain, state, action; options...)
transition!(domain::GenericDomain, state::State, action::Term; options...) =
    transition!(ConcreteInterpreter(), domain, state, action; options...)

available(domain::GenericDomain, state::State, action::Action, args) =
    available(ConcreteInterpreter(), domain, state, action, args)
available(domain::GenericDomain, state::State) =
    available(ConcreteInterpreter(), domain, state)

execute(domain::GenericDomain, state::State, action::Action, args; options...) =
    execute(ConcreteInterpreter(), domain, state, action, args; options...)
execute!(domain::GenericDomain, state::State, action::Action, args; options...) =
    execute!(ConcreteInterpreter(), domain, state, action, args; options...)

relevant(domain::GenericDomain, state::State, action::Action, args) =
    relevant(ConcreteInterpreter(), domain, state, action, args)
relevant(domain::GenericDomain, state::State) =
    relevant(ConcreteInterpreter(), domain, state)

regress(domain::GenericDomain, state::State, action::Action, args; options...) =
    regress(ConcreteInterpreter(), domain, state, action, args; options...)
regress!(domain::GenericDomain, state::State, action::Action, args; options...) =
    regress!(ConcreteInterpreter(), domain, state, action, args; options...)

update!(domain::GenericDomain, state::State, diff::Diff) =
    update!(ConcreteInterpreter(), domain, state, diff)
update(domain::GenericDomain, state::State, diff::Diff) =
    update(ConcreteInterpreter(), domain, state, diff)


text:
"""
Collection of utilities for working with Requests data
"""
module Requests

import Genie, Genie.Router, Genie.Input
import HTTP, Reexport

export jsonpayload, rawpayload, filespayload, postpayload, getpayload
export request, getrequest, matchedroute, matchedchannel, wsclient
export infilespayload, filename, payload, peer, currenturl

"""
    jsonpayload()

Processes an `application/json` `POST` request.
If it fails to successfully parse the `JSON` data it returns `nothing`. The original payload can still be accessed invoking `rawpayload()`
"""
function jsonpayload()
  haskey(Router.params(), Genie.PARAMS_JSON_PAYLOAD) ? Router.params(Genie.PARAMS_JSON_PAYLOAD) : nothing
end

"""
    jsonpayload(v)

Processes an `application/json` `POST` request attempting to return value corresponding to key v.
"""
function jsonpayload(v)
  jsonpayload()[v]
end

"""
    rawpayload() :: String

Returns the raw `POST` payload as a `String`.
"""
function rawpayload() :: String
  haskey(Router.params(), Genie.PARAMS_RAW_PAYLOAD) ? Router.params(Genie.PARAMS_RAW_PAYLOAD) : ""
end

"""
    filespayload() :: Dict{String,HttpFile}

Collection of form uploaded files.
"""
function filespayload() :: Dict{String,Input.HttpFile}
  haskey(Router.params(), Genie.PARAMS_FILES) ? Router.params(Genie.PARAMS_FILES) : Dict{String,Input.HttpFile}()
end

"""
    filespayload(filename::Union{String,Symbol}) :: HttpFile

Returns the `HttpFile` uploaded through the `key` input name.
"""
function filespayload(key::Union{String,Symbol}) :: Input.HttpFile
  Router.params(Genie.PARAMS_FILES)[string(key)]
end

"""
    infilespayload(key::Union{String,Symbol}) :: Bool

Checks if the collection of uploaded files contains a file stored under the `key` name.
"""
function infilespayload(key::Union{String,Symbol}) :: Bool
  haskey(filespayload(), string(key))
end

"""
    Base.write(filename::String = file.name; file::Input.HttpFile) :: IntBase.write(file::HttpFile, filename::String = file.name)

Writes uploaded `HttpFile` `file` to local storage under the `name` filename.
Returns number of bytes written.
"""
function Base.write(file::Input.HttpFile; filename::String = file.name) :: Int
  write(filename, IOBuffer(file.data))
end

"""
    read(file::HttpFile)

Returns the content of `file` as string.
"""
function Base.read(file::Input.HttpFile, ::Type{String}) :: String
  file.data |> String
end

"""
    filename(file::HttpFile) :: String

Original filename of the uploaded `HttpFile` `file`.
"""
function filename(file::Input.HttpFile) :: String
  file.name
end

"""
    postpayload() :: Dict{Symbol,Any}

A dict representing the POST variables payload of the request (corresponding to a `form-data` request)
"""
function postpayload() :: Dict{Symbol,Any}
  haskey(Router.params(), Genie.PARAMS_POST_KEY) ? Router.params(Genie.PARAMS_POST_KEY) : Dict{Symbol,Any}()
end

"""
    postpayload(key::Symbol) :: Any

Returns the value of the POST variables `key`.
"""
function postpayload(key::Symbol)
  postpayload()[key]
end

"""
    postpayload(key::Symbol, default::Any)

Returns the value of the POST variables `key` or the `default` value if `key` is not defined.
"""
function postpayload(key::Symbol, default::Any)
  haskey(postpayload(), key) ? postpayload(key) : default
end

"""
    getpayload() :: Dict{Symbol,Any}

A dict representing the GET/query variables payload of the request (the part correspoding to `?foo=bar&baz=moo`)
"""
function getpayload() :: Dict{Symbol,Any}
  haskey(Router.params(), Genie.PARAMS_GET_KEY) ? Router.params(Genie.PARAMS_GET_KEY) : Dict{Symbol,Any}()
end

"""
    getpayload(key::Symbol) :: Any

The value of the GET/query variable `key`, as in `?key=value`
"""
function getpayload(key::Symbol) :: Any
  getpayload()[key]
end

"""
    getpayload(key::Symbol, default::Any) :: Any

The value of the GET/query variable `key`, as in `?key=value`. If `key` is not defined, `default` is returned.
"""
function getpayload(key::Symbol, default::Any) :: Any
  haskey(getpayload(), key) ? getpayload(key) : default
end

"""
    request() :: HTTP.Request

Returns the raw HTTP.Request object associated with the request. If no request is available (not within a
request/response cycle) returns `nothing`.
"""
function request() :: Union{HTTP.Request,Nothing}
  haskey(Router.params(), Genie.PARAMS_REQUEST_KEY) ? Router.params(Genie.PARAMS_REQUEST_KEY) : nothing
end

const getrequest = request

"""
    payload() :: Any

Utility function for accessing the `params` collection, which holds the request variables.
"""
function payload() :: Dict
  Router.params()
end

"""
    payload(key::Symbol) :: Any

Utility function for accessing the `key` value within the `params` collection of request variables.
"""
function payload(key::Symbol) :: Any
  Router.params()[key]
end

"""
    payload(key::Symbol, default_value::T) :: Any

Utility function for accessing the `key` value within the `params` collection of request variables.
If `key` is not defined, `default_value` is returned.
"""
function payload(key::Symbol, default_value::T)::T where {T}
  haskey(Router.params(), key) ? Router.params()[key] : default_value
end

"""
    matchedroute() :: Route

Returns the `Route` object which was matched for the current request or `noting` if no route is available.
"""
function matchedroute() :: Union{Genie.Router.Route,Nothing}
  haskey(Router.params(), Genie.PARAMS_ROUTE_KEY) ? Router.params(Genie.PARAMS_ROUTE_KEY) : nothing
end

"""
    matchedchannel() :: Channel

Returns the `Channel` object which was matched for the current request or `nothing` if no channel is available.
"""
function matchedchannel() :: Union{Genie.Router.Channel,Nothing}
  haskey(Router.params(), Genie.PARAMS_CHANNELS_KEY) ? Router.params(Genie.PARAMS_CHANNELS_KEY) : nothing
end

"""
    wsclient() :: HTTP.WebSockets.WebSocket

The web sockets client for the current request or nothing if not available.
"""
function wsclient() :: Union{HTTP.WebSockets.WebSocket,Nothing}
  haskey(Router.params(), Genie.PARAMS_WS_CLIENT) ? Router.params(Genie.PARAMS_WS_CLIENT) : nothing
end

"""
    wtclient() :: HTTP.WebSockets.WebSocket

The web sockets client for the current request.
"""
function wtclient() :: UInt
  Router.params(:wtclient) |> hash
end

function getheaders(req::HTTP.Request) :: Dict{String,String}
  Dict{String,String}(req.headers)
end
function getheaders() :: Dict{String,String}
  getheaders(getrequest())
end

"""
    peer()

Returns information about the requesting client's IP address as a NamedTuple{(:ip,), Tuple{String}}
If the client IP address can not be retrieved, the `ip` field will return an empty string `""`.
"""
function peer() :: NamedTuple{(:ip,:port), Tuple{String,String}}
  unset_peer = (ip = "", port = "")

  if haskey(task_local_storage(), :peer)
    try
      (ip = string(task_local_storage(:peer)[1]), port = string(task_local_storage(:peer)[2]))
    catch ex
      @error ex
      unset_peer
    end
  else
    unset_peer
  end
end

"""
    isajax(req::HTTP.Request = getrequest()) :: Bool

Attempts to determine if a request is Ajax by sniffing the headers.
"""
function isajax(req::HTTP.Request = getrequest()) :: Bool
  for (k,v) in getheaders(req)
    k = replace(k, r"_|-"=>"") |> lowercase
    occursin("requestedwith", k) && occursin("xmlhttp", lowercase(v)) && return true
  end

  return false
end

"""
    currenturl() :: String

Returns the URL of the current page/request, starting from the path and including the query string.

### Example

```julia
currenturl()
"/products?promotions=yes"
```
"""
function currenturl(req::HTTP.Request = getrequest()) :: String
  req.target
end

end


text:
using GLVisualize, GLAbstraction

if !isdefined(:runtests)
    window = glscreen()
    timesignal = loop(linspace(0f0,1f0,360))
end

function juliadata(max_iterations, imgx, imgy)
    scalex, scaley = 4.0/imgx, 4.0/imgy

    # initialize our Float32 heightfield
    heightfield = zeros(Float32, imgx, imgy)

    # do julia set magic!
    for x=1:imgx, y=1:imgy
        cy = y * scaley - 2.0
        cx = x * scalex - 2.0
        z = Complex(cx, cy)
        c = Complex(-0.4, 0.6)
        i = 0
        for t in 0:max_iterations
            if norm(z) > 2.0
                break
            end
            z = z * z + c
            i = t
        end
        heightfield[x,y] = -(i/512f0)
    end
    heightfield
end

rotation_angle = const_lift(*, timesignal, 2f0*pi)
rotation = map(rotationmatrix_z, rotation_angle)

heightfield = juliadata(256, 700, 700)

# visualize the heightfield as a surface
vis = visualize(
    heightfield, :surface,
    model=rotation
)

# display it on the window
_view(vis, window)

if !isdefined(:runtests)
    renderloop(window) # render!
end


text:
module Readers

using Avro.DataFile
using Avro.DataFile.Codecs
using Avro.Generic
using Avro.Io
using Avro.Schemas

export wrap

struct Reader
    input_decoder::BinaryDecoder
    schema::Schemas.Schema
    codec::Codec
    sync_marker::Vector{UInt8}
end

function wrap(input::IO)
    input_decoder = BinaryDecoder(input)

    header = _read_header(input_decoder)
    magic = getindex(header, 1)
    meta = getindex(header, 2)
    sync_marker = getindex(header, 3).bytes

    # The magic header is invalid
    if magic.bytes != OBJECT_CONTAINER_FILE_MAGIC
        throw("Not an Avro data file.")
    end

    # Parse the schema from the metadata
    schema = Schemas.parse(meta[META_SCHEMA_KEY])
    codec = Codecs.create(get(meta, META_CODEC_KEY, "null"))

    Reader(input_decoder, schema, codec, sync_marker)
end

function Base.close(file_reader::Reader)
    close(file_reader.input_decoder.stream)
end

function start(file_reader::Reader)
    _read_block_header(file_reader)
end

# The invariant maintained by the iterator is that block_count = 0 only at the end of the
# stream. If block_count does drop to 0, we verify the 16 sync bytes and then read block
# headers until one with size > 0 is found.
function next(file_reader::Reader, state)
    buffer_decoder, block_count = state
    @assert block_count != 0
    item = read(buffer_decoder, file_reader.schema)
    block_count -= 1
    if block_count == 0
        sync = read(file_reader.input_decoder, SYNC_SCHEMA)
        @assert file_reader.sync_marker == sync.bytes
        buffer_decoder, block_count = _read_block_header(file_reader)
    end
    item, (buffer_decoder, block_count)
end

function done(state)
    decoder, block_count = state
    block_count == 0
end

function Base.iterate(file_reader::Reader, state)
    done(state) ? nothing : next(file_reader, state)
end
Base.iterate(file_reader::Reader) = iterate(file_reader, start(file_reader))

Base.IteratorSize(reader::Reader) = Base.SizeUnknown()

"""
Read the header directly the input decoder.
"""
function _read_header(input_decoder::BinaryDecoder)
    read(input_decoder, METADATA_SCHEMA)
end

# Always returns block_count > 0, except in the case of eof, when we return (nothing, 0)
# This allows us to keep the invariant that block_count is always > 0 unless at end of stream.
function _read_block_header(file_reader::Reader)
    block_count = 0
    block_data = nothing
    while block_count == 0
        if (eof(file_reader.input_decoder.stream))
            return (nothing, 0)
        end
        input_decoder = file_reader.input_decoder
        block_count = decode_long(input_decoder)
        num_bytes = decode_long(input_decoder)
        block_data = decode_fixed(input_decoder, num_bytes)
        block_data = Codecs.decompress(file_reader.codec, block_data)
        if block_count == 0
            sync = read(input_decoder, SYNC_SCHEMA)
            @assert file_reader.sync_marker == sync.bytes
        end
    end

    buffer_decoder = BinaryDecoder(IOBuffer(block_data))
    buffer_decoder, block_count
end

end


text:
function update_particles_noFEMcl!(mesh,u,v,pList,wList,param,data,cl)
    pUarr,pVarr = compute_particle_velocity_noFEMcl(mesh,u,v,pList,wList,param,data,cl)
    for ti=1:data.n_particles
        pList[ti].pos.x += pUarr[ti]*param.Δt
        pList[ti].pos.y += pVarr[ti]*param.Δt
    end

    nothing
end

# does not use cell lists
function compute_particle_velocity_noFEMcl(mesh,u,v,pList,wList,param,data,cl)
    StokesParticles.compute_gravity_force!(data.gfX,data.gfY,data.n_particles,param.G)
    compute_seepage_force_nocl!(data.sfX,data.sfY,mesh,u,v,pList,data)
    StokesParticles.compute_cohesion_force_cl!(data.cfX,data.cfY,cl,pList,param.sC,param.ϵ)
    StokesParticles.compute_adhesion_force!(data.afX,data.afY,pList,wList,data.pointOnWall,param.sA,param.ϵ)

    pUarr = data.gfX + data.cfX + data.afX + data.sfX
    pVarr = data.gfY + data.cfY + data.afY + data.sfY

    return pUarr,pVarr
end

function update_particles!(mesh,u,v,pList,wList,param,data,cl,femap)
    pUarr,pVarr = compute_particle_velocity(mesh,u,v,pList,wList,param,data,cl,femap)
    for ti=1:data.n_particles
        pList[ti].pos.x += pUarr[ti]*param.Δt
        pList[ti].pos.y += pVarr[ti]*param.Δt
    end

    nothing
end

# does not use cell lists
function compute_particle_velocity(mesh,u,v,pList,wList,param,data,cl,femap)
    StokesParticles.compute_gravity_force!(data.gfX,data.gfY,data.n_particles,param.G)
    compute_seepage_force!(data.sfX,data.sfY,mesh,u,v,pList,data,cl,femap)
    StokesParticles.compute_cohesion_force_cl!(data.cfX,data.cfY,cl,pList,param.sC,param.ϵ)
    StokesParticles.compute_adhesion_force!(data.afX,data.afY,pList,wList,data.pointOnWall,param.sA,param.ϵ)

    pUarr = data.gfX + data.cfX + data.afX + data.sfX
    pVarr = data.gfY + data.cfY + data.afY + data.sfY

    return pUarr,pVarr
end

max_stars_repo_name: pseastham/eFEMpart


text:
#=

  Firing squad example (from Judea Pearl)

  court_order: a court order to shoot a prisoner is given
  captain_signals: captain signal to a and b to shoot the prisoner
  a_shoots: person a shoots at the prisoner
  b_shoots: person b shoots at the prisoner
  death: the prisoner is dead (due to the shooting of a and/or b)

  Cf the ProbLog model firing_squad.pl

  See ~/cplint/firing_squad.pl
      ~/blog/firing_squad.blog
      ~/psi/firing_squad.psi
      ~/webppl/firing_squad.wppl

=#

import Turing, StatsPlots, DataFrames
include("jl_utils.jl")

@model firing_squad(court_order=missing,
                    death=missing,
                    captain_signals=missing,
                    a_shoots=missing,
                    b_shoots=missing) = begin

    court_order ~ Bernoulli(0.8)
    captain_signals ~ court_order ? Bernoulli(0.9) : Bernoulli(0.1)
    a_shoots ~ captain_signals ? Bernoulli(0.95) : Bernoulli(0.1)
    b_shoots ~ captain_signals ? Bernoulli(0.95) : Bernoulli(0.1)

    if a_shoots || b_shoots
        death ~ Bernoulli(0.9)
    else
        death ~ Bernoulli(0) # the prisoner don't die if not neither a or b shoots
    end

    # true ~ Dirac(court_order == false)
    # true ~ Dirac(captain_signals == false)
    # true ~ Dirac(a_shoots == true)
    # true ~ Dirac(b_shoots == true)
    true ~ Dirac(death == true)

end

court_order=missing
captain_signals=missing
a_shoots=missing # 'missing' is perhaps a bit confusing here :-)
b_shoots=missing
death=missing # missing
println("""
court_order:$court_order
captain_signals:$captain_signals
a_shoots:$a_shoots
b_shoots:$b_shoots
death:$death
""")
model = firing_squad(court_order,death,captain_signals,a_shoots,b_shoots)

num_chains = 4
# chains = sample(model, Prior(), MCMCThreads(), 1000, num_chains)

chains = sample(model, MH(), MCMCThreads(), 10_000, num_chains) # 1.47s
# chains = sample(model, IS(), MCMCThreads(), 10_000, num_chains) # 1.48s # don't handle Turing.@addlogprob! correct!

# chains = sample(model, PG(20), MCMCThreads(), 10_000, num_chains) # Too slow!
# chains = sample(model, SMC(), MCMCThreads(), 10_000, num_chains) # 4.1s, don't handle observes!

# chains = sample(model, HMC(0.05,10), MCMCThreads(), 10_000, num_chains) # Nope!
# chains = sample(model, NUTS(), MCMCThreads(), 10_000, num_chains) # Nope!

display(chains)

max_stars_repo_name: Wikunia/hakank

max_issues_repo_name: Wikunia/hakank

max_forks_repo_name: Wikunia/hakank


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

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

                     @variable(m, -16.957145096415093 <= q <= 17.94995800016443)

                     add_NL_constraint(m, :(swish(0.510608367415756 + 0.13576703913191768*swish(-0.023366318069798275 + -0.6972645804796413*swish(-0.6865389686816434 + 0.7515314314281274*$(x[1]) + -0.6334263786214325*$(x[2]) + -0.023277189078302563*$(x[3])) + 0.623708383159082*swish(0.5286761878087685 + -0.9557352407199864*$(x[1]) + -0.7485274707208927*$(x[2]) + -0.13145677322460392*$(x[3])) + -0.5526763214522168*swish(0.7131777426657733 + -0.00917591460146161*$(x[1]) + -0.5192676904312425*$(x[2]) + 0.597536428647615*$(x[3]))) + 0.3585964327428499*swish(-0.17574028037161238 + -0.933468976339102*swish(-0.6865389686816434 + 0.7515314314281274*$(x[1]) + -0.6334263786214325*$(x[2]) + -0.023277189078302563*$(x[3])) + -0.824606977579692*swish(0.5286761878087685 + -0.9557352407199864*$(x[1]) + -0.7485274707208927*$(x[2]) + -0.13145677322460392*$(x[3])) + -0.3761458682864913*swish(0.7131777426657733 + -0.00917591460146161*$(x[1]) + -0.5192676904312425*$(x[2]) + 0.597536428647615*$(x[3]))) + 0.2511872166134532*swish(0.9643857894710246 + 0.9348568947078797*swish(-0.6865389686816434 + 0.7515314314281274*$(x[1]) + -0.6334263786214325*$(x[2]) + -0.023277189078302563*$(x[3])) + -0.4027229019076284*swish(0.5286761878087685 + -0.9557352407199864*$(x[1]) + -0.7485274707208927*$(x[2]) + -0.13145677322460392*$(x[3])) + -0.6455757662896615*swish(0.7131777426657733 + -0.00917591460146161*$(x[1]) + -0.5192676904312425*$(x[2]) + 0.597536428647615*$(x[3])))) + swish(-0.7493006135688143 + 0.2483293387413994*swish(-0.023366318069798275 + -0.6972645804796413*swish(-0.6865389686816434 + 0.7515314314281274*$(x[1]) + -0.6334263786214325*$(x[2]) + -0.023277189078302563*$(x[3])) + 0.623708383159082*swish(0.5286761878087685 + -0.9557352407199864*$(x[1]) + -0.7485274707208927*$(x[2]) + -0.13145677322460392*$(x[3])) + -0.5526763214522168*swish(0.7131777426657733 + -0.00917591460146161*$(x[1]) + -0.5192676904312425*$(x[2]) + 0.597536428647615*$(x[3]))) + -0.7456041150082044*swish(-0.17574028037161238 + -0.933468976339102*swish(-0.6865389686816434 + 0.7515314314281274*$(x[1]) + -0.6334263786214325*$(x[2]) + -0.023277189078302563*$(x[3])) + -0.824606977579692*swish(0.5286761878087685 + -0.9557352407199864*$(x[1]) + -0.7485274707208927*$(x[2]) + -0.13145677322460392*$(x[3])) + -0.3761458682864913*swish(0.7131777426657733 + -0.00917591460146161*$(x[1]) + -0.5192676904312425*$(x[2]) + 0.597536428647615*$(x[3]))) + 0.49415154189847854*swish(0.9643857894710246 + 0.9348568947078797*swish(-0.6865389686816434 + 0.7515314314281274*$(x[1]) + -0.6334263786214325*$(x[2]) + -0.023277189078302563*$(x[3])) + -0.4027229019076284*swish(0.5286761878087685 + -0.9557352407199864*$(x[1]) + -0.7485274707208927*$(x[2]) + -0.13145677322460392*$(x[3])) + -0.6455757662896615*swish(0.7131777426657733 + -0.00917591460146161*$(x[1]) + -0.5192676904312425*$(x[2]) + 0.597536428647615*$(x[3])))) + swish(0.3105901915386422 + 0.5384893864274902*swish(-0.023366318069798275 + -0.6972645804796413*swish(-0.6865389686816434 + 0.7515314314281274*$(x[1]) + -0.6334263786214325*$(x[2]) + -0.023277189078302563*$(x[3])) + 0.623708383159082*swish(0.5286761878087685 + -0.9557352407199864*$(x[1]) + -0.7485274707208927*$(x[2]) + -0.13145677322460392*$(x[3])) + -0.5526763214522168*swish(0.7131777426657733 + -0.00917591460146161*$(x[1]) + -0.5192676904312425*$(x[2]) + 0.597536428647615*$(x[3]))) + 0.5405413918049242*swish(-0.17574028037161238 + -0.933468976339102*swish(-0.6865389686816434 + 0.7515314314281274*$(x[1]) + -0.6334263786214325*$(x[2]) + -0.023277189078302563*$(x[3])) + -0.824606977579692*swish(0.5286761878087685 + -0.9557352407199864*$(x[1]) + -0.7485274707208927*$(x[2]) + -0.13145677322460392*$(x[3])) + -0.3761458682864913*swish(0.7131777426657733 + -0.00917591460146161*$(x[1]) + -0.5192676904312425*$(x[2]) + 0.597536428647615*$(x[3]))) + -0.18067211353899282*swish(0.9643857894710246 + 0.9348568947078797*swish(-0.6865389686816434 + 0.7515314314281274*$(x[1]) + -0.6334263786214325*$(x[2]) + -0.023277189078302563*$(x[3])) + -0.4027229019076284*swish(0.5286761878087685 + -0.9557352407199864*$(x[1]) + -0.7485274707208927*$(x[2]) + -0.13145677322460392*$(x[3])) + -0.6455757662896615*swish(0.7131777426657733 + -0.00917591460146161*$(x[1]) + -0.5192676904312425*$(x[2]) + 0.597536428647615*$(x[3])))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
using Gadfly, DataFrames

xs = 0:0.1:20

df = DataFrame(
    x=xs,
    y=cos(xs),
    ymin=cos(xs) .- 0.5,
    ymax=cos(xs) .+ 0.5,
)

plot(df, x=:x, y=:y, ymin=:ymin, ymax=:ymax, Geom.line, Geom.ribbon)


text:
# ranges.jl

(::Colon)(start::QDate, stop::QDate) = StepRange(start, Day(1), stop)

# Given a start and end date, how many steps/periods are in between
Dates.guess(a::QDate, b::QDate, c) = Int64(div(value(b - a), days(c)))

function Dates.len(a::QDate, b::QDate, c)
    lo, hi, st = min(a, b), max(a, b), abs(c)
    i = Dates.guess(a, b, c) - 1
    try
        while lo + st * i <= hi
            i += 1
        end
    catch ex
        if !isa(ex, ArgumentError) && !isa(ex, BoundsError)
            rethrow()
        # else
        #     #pass
        end
    end
    return i - 1
end


text:
# Some key events are "one-shot" either triggered on the falling edge
# (i.e. key-down) or rising edge (i.e. key-up)
#
# Some key events are repeating (i.e. as long as a key is down the event occurs)
#
mutable struct KeyState
    thrusting::Bool
    firing::Bool
    turning::Int8
end

function firing(state::KeyState)
    s = state.firing
    state.firing = false
    s
end

function thrusting(state::KeyState)
    state.thrusting
end

function turning(state::KeyState)
    state.turning
end

function set_state!(state::KeyState, event::Events.KeyboardEvent)
    if (event.keycode == Events.KEY_A && event.state == Events.KEY_PRESSED)   # thrust
        state.thrusting = true
        return
    end

    if (event.keycode == Events.KEY_A && event.state == Events.KEY_RELEASED)   # thrust
        state.thrusting = false
        return
    end

    if (event.repeat != Events.KEY_REPEATING && event.keycode == Events.KEY_S && event.state == Events.KEY_PRESSED)   # fire
        state.firing = true
        return
    end

    # turn left
    if (event.keycode == Events.KEY_LEFT && event.state == Events.KEY_PRESSED)   # thrust
        state.turning = -1
        return
    end

    if (event.keycode == Events.KEY_LEFT && event.state == Events.KEY_RELEASED)   # thrust
        state.turning = 0
        return
    end

    # turn right
    if (event.keycode == Events.KEY_RIGHT && event.state == Events.KEY_PRESSED)   # thrust
        state.turning = 1
        return
    end

    if (event.keycode == Events.KEY_RIGHT && event.state == Events.KEY_RELEASED)   # thrust
        state.turning = 0
        return
    end

end

max_stars_repo_name: wdevore/Ranger-Julia-SDL

max_issues_repo_name: wdevore/Ranger-Julia-SDL

max_forks_repo_name: wdevore/Ranger-Julia-SDL


text:
using Test, YaoBlocks, YaoArrayRegister

@test shift(1) isa ShiftGate{Float64}

@testset "test constructor" for T in [Float16, Float32, Float64]
    @test ShiftGate(T(0.1)) isa ShiftGate{T}
    @test shift(T(0.1)) isa ShiftGate{T}
    @test adjoint(shift(0.1)) == shift(-0.1)
end

@testset "test matrix" begin
    g = ShiftGate{Float64}(pi)
    @test mat(g) ≈ exp(im * pi / 2) * [exp(-im * pi / 2) 0; 0 exp(im * pi / 2)]
end

@testset "test apply" begin
    g = ShiftGate{Float64}(pi)
    reg = rand_state(1)
    @test mat(g) * state(reg) ≈ state(apply!(reg, g))
end

@testset "test compare" begin
    @test (ShiftGate(2.0) == ShiftGate(2.0)) == true
    @test (ShiftGate(2.0) == ShiftGate(1.0)) == false
end

@testset "test copy" begin
    g = ShiftGate(0.1)
    cg = copy(g)
    cg.theta = 0.2
    @test g.theta == 0.1
end

@testset "test properties" begin
    g = ShiftGate{Float64}(0.1)
    @test nqubits(g) == 1
    @test isreflexive(g) == false
    @test isunitary(g) == true
    @test ishermitian(g) == false

    g = ShiftGate{ComplexF64}(0.1 + 0im)
    @test @test_nowarn isunitary(g) == true

    g = ShiftGate{ComplexF64}(0.1 + 1im)
    @test @test_logs (
        :warn,
        "θ in ShiftGate is not real, got θ=0.1 + 1.0im, fallback to matrix-based method",
    ) isunitary(g) == false
end

@testset "test parameters" begin
    @test nparameters(shift(0.1)) == 1
    @test parameters(shift(0.1)) == [0.1]
    @test parameters(dispatch!(shift(0.1), 0.2)) == [0.2]
    @test parameters(dispatch(shift(0.1), 2)) == [2]
end

max_stars_repo_name: QuantumBFS/YaoStdCircuit

max_issues_repo_name: QuantumBFS/YaoStdCircuit


text:
@testset "Composition" begin
    @testset "total degree" begin
        @polyvar a b c x y z u v
        e = [u + 1, v - 2]
        f = [a * b - 2, a * c - 1]
        g = [x + y, y + 3, x + 2]
        res = solve(e ∘ f ∘ g)
        @test nnonsingular(res) == 2
        @test nnonsingular(solve(e ∘ f ∘ g, system_scaling = nothing)) == 2
        @test nnonsingular(solve(e ∘ f ∘ g, system_scaling = :equations_and_variables)) == 2

        res = solve(e ∘ f ∘ g, system = SPSystem)
        @test nnonsingular(res) == 2
    end

    @testset "polyhedral" begin
        @polyvar x y u v
        result = solve(
            [2y + 3 * y^2 - x * y^3, x + 4 * x^2 - 2 * x^3 * y] ∘ [u, v];
            seed = 32241, start_system = :polyhedral,
        )
        @test nsolutions(result) == 6
    end

    @testset "parameter homotopy" begin
        @polyvar p q a b c x y z u v

        f = [a * b - 2, a * c - 1]
        g = [x + y, y + 3, x + 2]
        res = solve(f ∘ g; system = FPSystem, threading = false)
        @test nnonsingular(res) == 2
        # parameters at the end
        f2 = [a * b - q, a * c - p]
        g = [x + y, y + 3, x + 2]
        r = solve(
            f2 ∘ g,
            solutions(res);
            parameters = [p, q], p₁ = [1, 2], p₀ = [2, 3], threading = false,
        )
        @test nnonsingular(r) == 2

        # parameters at the beginning
        f = [a * b - 2, a * c - 1]
        g2 = [x + y, y + u, x + v]
        r = solve(
            f ∘ g2,
            solutions(res);
            parameters = [u, v], p₁ = [3, 2], p₀ = [-2, 3], threading = false,
        )
        @test nnonsingular(r) == 2

        # parameter in the middle
        e = [u + 1, v - 2]
        res2 = solve(e ∘ f ∘ g; system = SPSystem)
        f2 = [a * b - q, a * c - p]
        r = solve(
            e ∘ f2 ∘ g,
            solutions(res2);
            parameters = [p, q], p₁ = [1, 2], p₀ = [2, 3], threading = false,
        )
        @test nnonsingular(r) == 2
    end
end


text:
# using SimpleGraphs: vlist, elist, out_neighbors

# Type piracy, this should really go in SimpleGraphs
# OTOH there's not anything else this could really mean
# function SimpleGraphs.components(g::SimpleDigraph)
#     SimpleGraphs.components(convert(SimpleGraph, g))
# end

# function _before(g::SimpleDigraph, v)
#     parents = in_neighbors(g, v)
#     for i in parents
#         append!(parents, _before(g, i))
#     end
#     return parents
# end

# before(g::SimpleDigraph, v; inclusive = true) = inclusive ? push!(_before(g, v), v) : _before(g, v)
# function before(g::SimpleDigraph, vs...; inclusive = true)
#     parents = inclusive ? collect(vs) : Symbol[]
#     for v in vs
#         append!(parents, before(g, v, inclusive = inclusive))
#     end
#     if !inclusive
#         setdiff!(parents, vs)
#     end
#     return unique!(parents)
# end

# notbefore(g::SimpleDigraph, vs...; inclusive = false) = setdiff(vlist(g), before(g, vs...; inclusive = !inclusive))

# function _after(g::SimpleDigraph, v)
#     children = out_neighbors(g, v)
#     for i in children
#         append!(children, _after(g, i))
#     end
#     return children
# end

# after(g::SimpleDigraph, v; inclusive = true) = inclusive ? push!(_after(g, v), v) : _after(g, v)

# function after(g::SimpleDigraph, vs...; inclusive = true)
#     children = inclusive ? collect(vs) : Symbol[]
#     for v in vs
#         append!(children, after(g, v, inclusive = inclusive))
#     end
#     if !inclusive
#         setdiff!(children, vs)
#     end
#     return unique!(children)
# end

# notafter(g::SimpleDigraph, vs...; inclusive = false) = setdiff(vlist(g), after(g, vs...; inclusive = !inclusive))

# function assemblefrom(m::DAGModel, params, args)
#     theModule = getmodule(m)
#     m_init = DAGModel(theModule, args, NamedTuple(), NamedTuple(), nothing)
#     m = foldl(params; init=m_init) do m0,v
#         merge(m0, DAGModel(theModule, findStatement(m, v)))
#     end
#     return m
# end