text:
using .DifferentialEquations
const DE = DifferentialEquations

using HybridSystems: states

# extend the solve API for initial-value problems

# =====================================
# Continuous system
# =====================================

function _solve_ensemble(ivp::InitialValueProblem, args...;
                         trajectories_alg=DE.Tsit5(),
                         ensemble_alg=DE.EnsembleThreads(),
                         inplace=true,
                         initial_states=nothing,
                         kwargs...)

    # get problem's vector field
    if inplace
        field = inplace_field!(ivp)
    else
        field = outofplace_field(ivp)
    end

    # get time span
    dt = _get_tspan(args...; kwargs...)
    tspan = (inf(dt), sup(dt))

    if isnothing(initial_states)
        # sample initial states
        X0 = initial_state(ivp)
        trajectories = get(kwargs, :trajectories, 100)
        initial_states = _sample_initial(X0, trajectories; kwargs...)
        # number of trajectories may increase if vertices got included
        trajectories = length(initial_states)
    else
        trajectories = length(initial_states)
    end

    # formulate ensemble ODE problem
    ensemble_prob = ODEProblem(field, first(initial_states), tspan)
    _prob_func(prob, i, repeat) = remake(prob, u0 = initial_states[i])

    # choose tolerances
    reltol = get(kwargs, :reltol, 1e-3)
    abstol = get(kwargs, :abstol, 1e-6)

    callback = get(kwargs, :callback, nothing)
    dtmax = get(kwargs, :dtmax, Inf)

    ensemble_prob = EnsembleProblem(ensemble_prob, prob_func = _prob_func)

    if isnothing(callback)
        result = DE.solve(ensemble_prob, trajectories_alg, ensemble_alg;
                          trajectories=trajectories, reltol=reltol,
                          abstol=abstol, dtmax=dtmax)
    else
        result = DE.solve(ensemble_prob, trajectories_alg, ensemble_alg;
                          trajectories=trajectories, reltol=reltol,
                          abstol=abstol, dtmax=dtmax, callback=callback)
    end
    return result
end

function _sample_initial(X0, trajectories; kwargs...)
    sampler = get(kwargs, :sampler, LazySets._default_sampler(X0))
    rng = get(kwargs, :rng, LazySets.GLOBAL_RNG)
    seed = get(kwargs, :seed, nothing)
    include_vertices = get(kwargs, :include_vertices, false)
    return sample(X0, trajectories; sampler=sampler, rng=rng,
                  seed=seed, include_vertices=include_vertices)
end

# =====================================
# Hybrid system
# =====================================

function _solve_ensemble(ivp::InitialValueProblem{<:AbstractHybridSystem},
                         args...; kwargs...)
    H = system(ivp)

    kwargs_sim = copy(kwargs)
    if :T in keys(kwargs_sim)
        # remove `T` because it is not allowed to define several time options
        pop!(kwargs_sim, :T)
    end

    time_span = _get_tspan(args...; kwargs...)
    t0_g = tstart(time_span)
    T = tend(time_span)

    termination_action = (integrator) -> terminate!(integrator)
    use_discrete_callback = get(kwargs, :use_discrete_callback, false)
    dtmax = get(kwargs, :dtmax, Inf)
    if use_discrete_callback && isinf(dtmax)
        # use a maximum time step when using discrete callbacks
        dtmax = DTMAX_SIM_DEFAULT
        kwargs_sim[:dtmax] = dtmax
    end

    jump_probability = get(kwargs, :jump_probability, 0.9)
    max_jumps = get(kwargs, :max_jumps, 100)

    # hybrid loop
    all_trajectories = []
    for initial_state in _sample_initial(ivp; kwargs...)
        loc, x0 = initial_state
        t0 = t0_g
        jump_count = 0
        trajectory_chain = []
        while true
            S_loc = mode(H, loc)
            ivp_loc = IVP(S_loc, x0)
            I⁻ = stateset(S_loc)

            # obtain trajectory in current mode
            if use_discrete_callback
                # use a discrete callback function based on membership
                condition = (u, t, integrator) -> u ∉ I⁻
                callback = DiscreteCallback(condition, termination_action)
            else
                # use a continuous callback function based on membership
                condition = (u, t, integrator) -> u ∉ I⁻ ? -1 : 1
                callback = ContinuousCallback(condition, termination_action)
            end
            trajectory = _solve_ensemble(ivp_loc, args...; initial_states=[x0],
                                         callback=callback, tspan=t0 .. T,
                                         kwargs_sim...)[1]

            # remove states from trajectory if they violate invariant
            # observation: this only concerns the last two steps
            if use_discrete_callback
                start_idx = max(1, length(trajectory.u) - 1)
                for idx in start_idx:length(trajectory.u)
                    if trajectory.u[idx] ∉ I⁻
                        indices_to_remove = idx:length(trajectory.u)
                        deleteat!(trajectory.u, indices_to_remove)
                        deleteat!(trajectory.t, indices_to_remove)
                        deleteat!(trajectory.k, indices_to_remove)
                        break
                    end
                end
            end

            # decide whether another jump should be attempted
            if jump_count > max_jumps
                # exceeded the allowed number of jumps
                try_jump = false
            elseif trajectory.t[end] >= T
                # time horizon reached → potentially do not take another jump
                try_jump = rand() <= jump_probability
            else
                try_jump = true
            end

            if try_jump
                # collect potential successors
                successors = []
                for t in out_transitions(H, loc)
                    G = guard(H, t)
                    asgn = resetmap(H, t)
                    loc′ = target(H, t)
                    I⁺ = stateset(H, loc′)

                    # find simulation steps that are enabled
                    # collect pairs (i, t, y) of indices i when transition t is
                    # enabled and leads to state y
                    for (i, xi) in enumerate(trajectory.u)
                        # check guard intersection
                        xi ∉ G && continue

                        # apply assignment
                        yi = MathematicalSystems.apply(asgn, xi)

                        # check target invariant
                        yi ∉ I⁺ && continue

                        # can take the transition
                        # store index, transition, and state
                        storage = (i, t, yi)
                        push!(successors, storage)
                    end
                end

                if isempty(successors)
                    jump_successful = false
                else
                    # choose a random successor
                    idx, t, y = rand(successors)
                    indices_to_remove = idx+1:length(trajectory.u)
                    deleteat!(trajectory.u, indices_to_remove)
                    deleteat!(trajectory.t, indices_to_remove)
                    deleteat!(trajectory.k, indices_to_remove)
                    loc′ = target(H, t)
                    push!(trajectory.u, y)
                    push!(trajectory.t, trajectory.t[end])
                    push!(trajectory.k, trajectory.k[end])
                    jump_successful = true
                end
            else
                jump_successful = false
            end

            # append new trajectory suffix to old trajectory prefix
            push!(trajectory_chain, trajectory)

            # test for termination
            if !jump_successful
                push!(all_trajectories, _merge_trajectory_chain(trajectory_chain))
                break
            end

            # continue previous simulation
            jump_count += 1
            loc = loc′
            x0 = trajectory.u[end]
            t0 = trajectory.t[end]
        end
    end
    return all_trajectories
end

# sample initial states of hybrid system given a list of pairs (loc, X0)
function _sample_initial(ivp::IVP{<:AbstractHybridSystem,
                                  <:Vector{<:Tuple{Integer, AdmissibleSet}}};
                         kwargs...)
    H = system(ivp)
    X0 = initial_state(ivp)

    # sample initial states from all possible initial regions
    trajectories = get(kwargs, :trajectories, 100)
    all_samples = []
    for (loc, X0_loc) in X0
        inv = stateset(mode(H, loc))
        X0_loc = X0_loc ∩ inv
        if isempty(X0_loc)
            continue
        end
        for x0 in _sample_initial(X0_loc, trajectories; kwargs...)
            push!(all_samples, (loc, x0))
        end
    end

    # filter for correct number of trajectories
    include_vertices = get(kwargs, :include_vertices, false)
    if !include_vertices
        # with k initial regions, we collected up to `k * trajectories` many
        # samples (some regions may be empty), so we reduce to a random
        # collection of `trajectories` samples
        rng = get(kwargs, :rng, LazySets.GLOBAL_RNG)
        all_samples = LazySets.Random.shuffle!(rng, all_samples)[1:trajectories]
    end

    return all_samples
end

# sample initial states of hybrid system given a set X0
function _sample_initial(ivp::IVP{<:AbstractHybridSystem, <:AdmissibleSet};
                         kwargs...)
    H = system(ivp)
    X0 = initial_state(ivp)
    X0_distributed = [(loc, X0) for loc in states(H)]

    # filter invalid combinations (empty intersection with invariant)
    i = 1
    while i <= length(X0_distributed)
        loc, X0 = X0_distributed[i]
        inv = stateset(H, loc)
        if isdisjoint(inv, X0)
            deleteat!(X0_distributed, i)
        else
            i += 1
        end
    end

    return _sample_initial(IVP(H, X0_distributed))
end

# merge trajectory pieces into an ODESolution object
function _merge_trajectory_chain(trajectory_chain)
    prefix = trajectory_chain[1]
    @inbounds for i in 2:length(trajectory_chain)
        infix = trajectory_chain[i]
        append!(prefix.u, infix.u)
        append!(prefix.t, infix.t)
        append!(prefix.k, infix.k)
    end
    return prefix
end


text:
#=

  Monks and doors problem in ConstraintSolver.jl

  From http://user.it.uu.se/~rolandb/LP/gammal/960615_facit.ps
  """
  There is a room with four doors and eight monks. One or more of
  the doors may be exit. Each monk is either telling a lie or the truth.

  The monks make the following statements:
  Monk 1: Door A is the exit.
  Monk 2: At least one of the doors B and C is the exit.
  Monk 3: Monk 1 and Monk 2 are telling the truth.
  Monk 4: Doors A and B are both exits.
  Monk 5: Doors A and B are both exits.
  Monk 6: Either Monk 4 or Monk 5 is telling the truth.
  Monk 7: If Monk 3 is telling the truth, so is Monk 6.
  Monk 8: If Monk 7 and Monk 8 are telling the truth, so is Monk 1.

  Which door is an exit no matter who is a liar and who is telling the
  truth.
  """

  Answer: Door A is an exit.
          And monks 1, 7, and 8 are telling the truth.

  Basic Boolean algebra, needed since ConstratinSolver.jl don't (yet)
  support more advance constraints in a reification:
     A /\ B  :  A + B == 2
     A \/ B  :  A + B >= 1
     A xor B :  A + B == 1
     A => B  :  A <= B
     A <=> B :  A == B

  Model created by Hakan Kjellerstrand, hakank@gmail.com
  See also my Julia page: http://www.hakank.org/julia/

=#

using ConstraintSolver, JuMP
using Cbc, GLPK, Ipopt
const CS = ConstraintSolver
include("constraints_utils.jl")

function monks_and_doors(print_solutions=true,all_solutions=true,timeout=6)

    cbc_optimizer = optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 0)
    glpk_optimizer = optimizer_with_attributes(GLPK.Optimizer)
    ipopt_optimizer = optimizer_with_attributes(Ipopt.Optimizer)

    model = Model(optimizer_with_attributes(CS.Optimizer,   "all_solutions"=> all_solutions,
                                                            # "all_optimal_solutions"=>all_solutions,
                                                            "logging"=>[],

                                                            "traverse_strategy"=>:BFS,
                                                            # "traverse_strategy"=>:DFS,
                                                            # "traverse_strategy"=>:DBFS,

                                                            # "branch_split"=>:Smallest,
                                                            # "branch_split"=>:Biggest,
                                                            "branch_split"=>:InHalf,

                                                            # https://wikunia.github.io/ConstraintSolver.jl/stable/options/#branch_strategy-(:Auto)
                                                            "branch_strategy" => :IMPS, # default
                                                            # "branch_strategy" => :ABS, # Activity Based Search
                                                            # "activity.decay" => 0.999, # default 0.999
                                                            # "activity.max_probes" => 10, # default, 10
                                                            # "activity.max_confidence_deviation" => 20, # default 20

                                                            # "simplify"=>false,
                                                            # "simplify"=>true, # default

                                                            "time_limit"=>timeout,

                                                            # "backtrack" => false, # default true
                                                            # "backtrack_sorting" => false, # default true

                                                            # "lp_optimizer" => cbc_optimizer,
                                                            # "lp_optimizer" => glpk_optimizer,
                                                            # "lp_optimizer" => ipopt_optimizer,
                                        ))

    num_doors = 4
    num_monks = 8
    @variable(model, doors[1:num_doors], Bin)
    da,db,dc,dd = doors
    door_names = ["A","B","C","D"]

    @variable(model, monks[1:num_monks], Bin)
    m1,m2,m3,m4,m5,m6,m7,m8 = monks

    # Monk 1: Door A is the exit.
    # M1 #= A (Picat constraint)
    @constraint(model, m1 == da)

    #  Monk 2: At least one of the doors B and C is the exit.
    # M2 #= 1 #<=> (B #\/ C)
    # @constraint(model, m2 := { db + dc >= 1})
    @constraint(model, m2 := { db == 1 || dc == 1})

    #  Monk 3: Monk 1 and Monk 2 are telling the truth.
    # M3 #= 1 #<=> (M1 #/\ M2)
    # @constraint(model, m3 := { m1 + m2 == 2})
    # @constraint(model, m3 := { m1 == 1 && m2 == 1})
    @constraint(model, m3 := { m1 == 1 && m2 == 1})

    #  Monk 4: Doors A and B are both exits.
    # M4 #= 1 #<=> (A #/\ B)
    # @constraint(model, m4 := { da + db == 2})
    @constraint(model, m4 := { da == 1 && db == 1})

    #  Monk 5: Doors A and C are both exits.
    # M5 #= 1 #<=> (A #/\ C)
    # @constraint(model, m5 := { da + dc == 2})
    @constraint(model, m5 := { da == 1 && dc == 1})

    #  Monk 6: Either Monk 4 or Monk 5 is telling the truth.
    # M6 #= 1 #<=> (M4 #\/ M5)
    # @constraint(model, m6 := { m4 + m5 == 1})
    @constraint(model, m6 := { m4 == 1|| m5 == 1})

    #  Monk 7: If Monk 3 is telling the truth, so is Monk 6.
    # M7 #= 1 #<=> (M3 #=> M6)
    @constraint(model, m7 := { m3 <= m6})
    # @constraint(model, m7 := { m3 => m6})  # This don't work!
    # @constraint(model, m7 := { m3 => {m6 == 1}}) # This don't work either.

    #  Monk 8: If Monk 7 and Monk 8 are telling the truth, so is Monk 1.
    # M8 #= 1 #<=> ((M7 #/\ M8) #=> (M1))
    b1 = @variable(model, binary=true)
    # @constraint(model, b1 := {m7 + m8 == 2})
    @constraint(model, b1 := {m7 == 1 && m8 == 1})
    @constraint(model, m8 := {b1 <= m1})

    # Exactly one door is an exit.
    # (A + B + C + D) #= 1
    @constraint(model, da + db + dc + dd == 1)

    # Solve the problem
    optimize!(model)

    status = JuMP.termination_status(model)
    # println("status:$status")
    num_sols = 0
    if status == MOI.OPTIMAL
        num_sols = MOI.get(model, MOI.ResultCount())
        println("num_sols:$num_sols\n")
        if print_solutions
            for sol in 1:num_sols
                #println("solution #$sol")
                doors_val = convert.(Integer,JuMP.value.(doors; result=sol))
                monks_val = convert.(Integer,JuMP.value.(monks; result=sol))
                println("doors:$doors_val  monks:$monks_val")

                println("the exit door          : ", [door_names[i] for i in 1:num_doors if doors_val[i] == 1 ][1] )
                println("monks telling the truth: ", [m for m in 1:num_monks if monks_val[m] == 1 ] )

            end
        end
    else
        println("status:$status")
    end

    return status, num_sols
end

@time monks_and_doors()

ext: jl

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule Qt5Location_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("Qt5Location")
JLLWrappers.@generate_main_file("Qt5Location", UUID("be351a2a-7d1f-55cd-a30a-1b8d4bc38a20"))
end  # module Qt5Location_jll


text:
### A Pluto.jl notebook ###
# v0.17.7

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local iv = try Base.loaded_modules[Base.PkgId(Base.UUID("6e696c72-6542-2067-7265-42206c756150"), "AbstractPlutoDingetjes")].Bonds.initial_value catch; b -> missing; end
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : iv(el)
        el
    end
end

# ╔═╡ f63c04d4-eefe-11eb-3fda-1729ac6de2cb
begin
	using Pkg
	Pkg.activate(joinpath(Pkg.devdir(), "MLCourse"))
    using PlutoUI
    PlutoUI.TableOfContents()
end

# ╔═╡ f63c04e8-eefe-11eb-39f6-83b31ebe73e7
using Plots, OpenML, DataFrames, CSV

# ╔═╡ f63c0540-eefe-11eb-01d6-454a05f9182a
using StatsPlots

# ╔═╡ f63c0592-eefe-11eb-0d69-91d4d396c540
# using MLJ makes the functions machine, fit!, fitted_params, predict, etc. available
# MLJLinearModel contains LinearRegressor, LogisticClassifier etc.
using MLJ, MLJLinearModels

# ╔═╡ f63c05e4-eefe-11eb-0848-9b4eb61d920c
using Random

# ╔═╡ 36a44a52-8a29-4f92-8136-2c66a81a1f5f
using Distributions

# ╔═╡ 0cfbd543-8b9b-406e-b3b4-c6cafbbec212
begin
    using MLCourse
    MLCourse.list_notebooks(@__FILE__)
end

# ╔═╡ f63c04dc-eefe-11eb-1e24-1d02a686920a
md"In this notebook we have a first look at some data sets and the MLJ machine learning package. You may want to look at its [documentation](https://alan-turing-institute.github.io/MLJ.jl/dev) or [cheatsheet](https://alan-turing-institute.github.io/MLJ.jl/dev/mlj_cheatsheet/).

# Data Sets for Supervised Learning

## MNIST Handwritten Digit Recognition"

# ╔═╡ f63c04e8-eefe-11eb-1a14-8305504a6f1c
mnist = OpenML.load(554) |> DataFrame;

# ╔═╡ f63c04f2-eefe-11eb-2562-e530426e4300
md"With the semicolon we omit the output of a cell.
"

# ╔═╡ f63c04f2-eefe-11eb-3ceb-ff1e36a2a302
mnist.class # this is a vector of 70'000 integers between 0 and 9.

# ╔═╡ f63c04fc-eefe-11eb-35b6-5345dda134e7
size(mnist)

# ╔═╡ f63c04fc-eefe-11eb-043c-7fec2fc41913
md"The input of this data set consists of 70'000 grayscale images of 28x28 pixels.

We can plot different input images with the following code.
"

# ╔═╡ f63c0506-eefe-11eb-1857-e3eaf731c152
plot(Gray.(reshape(Array(mnist[19, 1:end-1]) ./ 255, 28, 28)'))

# ╔═╡ c82fb649-e60d-43fc-b158-444b27b6e9cf
md"Explanation: `data[19, 1:end-1]` gets all pixels of the 19th image.
Next we transform the dataframe row to an `Array`. In OpenML, the inputs are
 8-bit integer; we transform to values between 0 and 1 by dividing by 255.
Then we reshape the vector of 784 entries to a matrix of size 28x28 with `reshape`.
The symbol `'` transposes the matrix (alternatively we could use `PermutedDimsArray`; look it up in the Live docs, if you are interested). `Gray.(...)` applies the `Gray` function elementwise to all the pixels, to reinterpret these numbers as grayscale values."

# ╔═╡ f63c0506-eefe-11eb-399f-455ee06548da
mnist.class[19] # the label corresponding to image 19 is a 6 :)

# ╔═╡ f63c050e-eefe-11eb-1709-0ddce2aee0ed
Markdown.parse("## Spam Detection

We load the spam data from a csv file on the harddisk to a DataFrame.
Our goal is to classify emails as spam or \"ham\".
The email texts are preprocessed. Html tags, numbers and punctuation is removed,
all characters are transformed to lowercase and words are reduced to their
[stem](https://en.wikipedia.org/wiki/Word_stem).  If you are interested in the
details of this preprocessing, you find more information about text
preprocessing
[here](https://juliahub.com/docs/TextAnalysis/5Mwet/0.7.3/documents/) and the
preprocessing code in the scripts folder of the MLCourse package
($(normpath(joinpath(dirname(pathof(MLCourse)), "..", "scripts", "preprocess_spam.jl")))).
")

# ╔═╡ f63c051a-eefe-11eb-0db2-c519292f88a2
spam = CSV.read(joinpath(@__DIR__, "..", "data", "spam.csv"), DataFrame)

# ╔═╡ f63c051a-eefe-11eb-16eb-a572fe41aa43
spam[spam.label .== "ham", :] # select only ham

# ╔═╡ f63c0524-eefe-11eb-3abd-63d677b12db9
spam[spam.label .== "spam", :] # select only spam

# ╔═╡ f63c0524-eefe-11eb-3732-258d6989e217
md"## Wind Speed Prediction

We load the weather data from a csv file on the harddisk to a DataFrame.
Our goal is to predict the wind speed in Lucerne from 5 hours old measurements.
"

# ╔═╡ f63c052e-eefe-11eb-3a14-e5e8f3d578a8
weather = CSV.read(joinpath(@__DIR__, "..", "data", "weather2015-2018.csv"),
                   DataFrame)[1:10000,:] # to speed up plotting, we use only the first 10'000 rows of the data set. Change `10000` to `end`, if you want to see all data.

# ╔═╡ f63c052e-eefe-11eb-0884-7bd433ce0c5e
y = weather.LUZ_wind_peak[6:end] # from hour 6 until the end we take all wind values

# ╔═╡ f63c0538-eefe-11eb-2a4b-d19efe6b6689
X = weather.LUZ_pressure[1:end-5] # pressure from the first hour until the fifth last hour

# ╔═╡ f63c0538-eefe-11eb-2eca-2d8bf12fae95
histogram2d(X, y, markersize = 3, xlabel = "LUZ_pressure [hPa]",
	        legend = false, ylabel = "LUZ_wind_peak [km/h]", bins = (250, 200))

# ╔═╡ f63c0538-eefe-11eb-2808-0f32a2fa84cf
md"The [StatsPlots package](https://github.com/JuliaPlots/StatsPlots.jl) has some
nice functions to explore data visually. We will use the `@df` macro, which allows
to pass columns as symbols (note that the symbols are not separated by a comma).
Have a look at the [readme of the StatsPlots package](https://github.com/JuliaPlots/StatsPlots.jl) for more examples."

# ╔═╡ f63c054c-eefe-11eb-13db-171084b417a9
@df weather corrplot([:BAS_pressure :LUG_pressure :LUZ_pressure :LUZ_wind_peak],
                     grid = false, fillcolor = cgrad(), size = (700, 600))

# ╔═╡ f63c0560-eefe-11eb-0821-73afecf00357
md"## Other Data Sets

There are many data sets collected by the [OpenML community](https://www.openml.org).
We use them as follows."

# ╔═╡ f63c0560-eefe-11eb-0251-2d0443448e66
OpenML.list_datasets(output_format = DataFrame)

# ╔═╡ f63c056a-eefe-11eb-199e-2b05b1385be4
OpenML.describe_dataset(61)

# ╔═╡ f63c056a-eefe-11eb-218c-afc51fdf1704
iris = OpenML.load(61) |> DataFrame

# ╔═╡ f63c0574-eefe-11eb-2d7c-bd74b7d83f23
md"# Function Family

In the following example we use the function family
    ``\hat y = f_\theta(x) = \theta_0 + \theta_1 x`` from the slides.
You can control the parameters ``\theta_0`` and ``\theta_1`` to find the
function that visually matches the data well. The loss is the mean squared error. This corresponds to the average squared length of the red vertical lines. On the right we track the loss of the currently defined function.
"

# ╔═╡ f63c057e-eefe-11eb-1672-19db9efd0d32
md"θ₀ = $(@bind θ₀ Slider(-1.8:.1:1.8, default = 0, show_value = true))"

# ╔═╡ f63c057e-eefe-11eb-0dd2-3333b4e2d220
md"θ₁ = $(@bind θ₁ Slider(-3:.1:3, default = 0, show_value = true))"

# ╔═╡ f63c0588-eefe-11eb-0687-a7c44f9177a4
let x = [0.62, 0.49, 0.28, 0.54, 0.0, 0.84, 0.29, 0.12, 0.49, 0.37, 0.36, 0.8, 0.63, 0.18, 0.2, 0.43, 0.68, 0.23, 0.72],
    y = [0.0, 0.35, -0.36, 0.31, -0.94, 1.0, -0.54, -1.35, -0.12, -0.45, 0.28, 0.48, 0.57, -0.48, -1.05, 0.32, 0.35, -0.39, 1.07]
    scatter(x, y,
            yrange = (-2, 2), xrange = (-.1, 1.1),
            xlabel = "x", ylabel = "ŷ", label = "data")
    p1 = plot!(x -> θ₀ + θ₁ * x, label = "function", title = "data & function")
	plot_residuals!(x, y, x -> θ₀ + θ₁ * x)
    p2 = scatter([θ₀], [θ₁], label = "current loss", color = :red)
    p2 = contour!(-1.8:.1:1.8, -3:.1:3, (θ₀, θ₁) -> mean(abs2, θ₀ .+ θ₁ * x .- y),
                 xlabel = "θ₀", ylabel = "θ₁", title = "loss", legend = :bottomright,
                 contour_labels = true, levels = collect(2 .^ -3:.5:12),
                 color = cgrad([:green, :red], [0, 10]), cbar = false,
                 )
    plot(p1, p2,
         layout = (1, 2), size = (700, 400))
end

# ╔═╡ f63c0588-eefe-11eb-1b0b-db185fcbdc4e
md"# Linear Regression

## Blackboard Example
In the following we run in code the same example we had on the blackboard."

# ╔═╡ f63c05a6-eefe-11eb-0ef4-df7e458e9ec6
md"For now we use the same data set as in the blackboard example."

# ╔═╡ f63c05a6-eefe-11eb-2927-b1e9dbbe032d
training_data = DataFrame(x = [0., 2., 2.], y = [-1., 4., 3.])

# ╔═╡ f63c05b2-eefe-11eb-03cd-ef71f9055eb3
begin
    scatter(training_data.x, training_data.y,
            xrange = (-1, 3), label = "training data", legend = :topleft)
    plot!(x -> 2x - 1, label = "mean data generator")
end

# ╔═╡ f63c05b2-eefe-11eb-21dc-1f79ae779316
md"Now we define a supervised learning machine. The function family, loss
   function and optimiser are implicitly defined by choosing `LinearRegressor()`.
   All machines for supervised learning have as arguments a method, the input
   and the output."

# ╔═╡ f63c05ba-eefe-11eb-18b5-7522b326ab65
mach = machine(LinearRegressor(),         # model
               select(training_data, :x), # input
               training_data.y);          # output

# ╔═╡ f63c05ba-eefe-11eb-1b7a-21cbd9bbeb37
fit!(mach); # fit the machine

# ╔═╡ f63c05c4-eefe-11eb-030e-6f7484849852
fitted_params(mach) # show the result

# ╔═╡ 25eb24a4-5414-481b-88d5-5ad50d75f8c9
md"We can use the fitted machine and the `predict` function to compute ŷ (the predictions)."

# ╔═╡ bd9a7214-f39e-4ae5-995a-6ec02d613fda
predict(mach) # predictions on the training input

# ╔═╡ f0f6ac8f-2e61-4b08-9383-45ff5b02c5b1
predict(mach, DataFrame(x = [0.5, 1.5])) # the second argument is the new input.

# ╔═╡ 7336e7b4-a7c2-4409-8d9d-00285cf633fb
md"Next, we define a function to compute the mean squared error for a given machine `mach` and a data set `data`."

# ╔═╡ f63c05c4-eefe-11eb-2d60-f7bf7beb8118
L(mach, data) = mean((predict(mach, select(data, :x)) .- data.y).^2);

# ╔═╡ f63c05ce-eefe-11eb-1ccd-c9b7714f3bd5
L(mach, training_data) # this is the training data

# ╔═╡ 8fe94709-3673-44cc-8702-83bd7e2cad51
md"We can get the same result by squaring the result obtained with the MLJ function `rmse` (root mean squared error)."

# ╔═╡ 27f44b28-f6d6-4e1f-8dbe-6fe7789b9a18
rmse(predict(mach), training_data.y)^2

# ╔═╡ 43b28121-5d07-4400-8710-287de7b978a4
md"For plotting we will define a new function `fitted_linear_func` that extracts the fitted parameters from a machine and returns a linear function with these parameters."

# ╔═╡ f63c05d8-eefe-11eb-2a11-fd8f954bf059
function fitted_linear_func(mach)
    θ̂ = fitted_params(mach)
    θ̂₀ = θ̂.intercept
    θ̂₁ = θ̂.coefs[1][2]
    x -> θ̂₀ + θ̂₁ * x
end;

# ╔═╡ f63c05d8-eefe-11eb-0a2e-970192b02d61
begin
    scatter(training_data.x, training_data.y,
            xrange = (-1, 3), label = "training data", legend = :topleft)
    plot!(x -> 2x - 1, label = "mean data generator")
	plot!(fitted_linear_func(mach), color = :green, label = "fit")
end

# ╔═╡ f63c05e4-eefe-11eb-012a-9bae1d87f2b5
md"Next we actually generate data from the data generator defined in the blackboard example."

# ╔═╡ f63c0592-eefe-11eb-2a76-15de55eff3ad
function data_generator(x, σ)
	y = 2x .- 1 .+ σ * randn(length(x))
	DataFrame(x = x, y = y)
end;

# ╔═╡ ca8e60e6-1d12-4076-8684-9a7c9e5ac509
md"We can use this generator to create a large test set and compute empirically the test loss of the fitted machine."

# ╔═╡ 4987ad0c-d100-4a5d-a725-12f23bcc7aa3
L(mach, data_generator(fill(1.5, 10^4), .5)) # test loss at 1.5 for σ = 0.5

# ╔═╡ 93f529ab-a386-4746-94de-63a536b5e2aa
L(mach, data_generator(randn(10^4), .5)) # test loss of the joint generator

# ╔═╡ 673a773c-2f6b-4676-8631-49ba08ec28a7
md"
## Training and test loss as a function of $n$ and $\sigma$
We will now look at how the training and the test loss depend on the number of samples $n$ and the noise of the data generator $\sigma$.

For reproducibility we fix the seed of the random number generator with the slider below. To do so we load the package `Random`."

# ╔═╡ f63c059c-eefe-11eb-208c-d5781711deb7
md"With the sliders below you can define the noise level of the data generator, the number n of training points and the seed of the (pseudo-)random number generator.

σ = $(@bind σ Slider(.1:.1:2, show_value = true))

n = $(@bind n Slider(2:5:100, default = 52, show_value = true))

seed = $(@bind seed Slider(1:30, show_value = true))
"

# ╔═╡ f63c05ec-eefe-11eb-0a92-215711928d02
begin
    Random.seed!(seed)
    generated_data = data_generator(randn(n), σ)
	generated_test_data = data_generator(-1:.01:3, σ)
	mach2 = machine(LinearRegressor(),
	            select(generated_data, :x),
	            generated_data.y)
	fit!(mach2)
end;

# ╔═╡ f63c0600-eefe-11eb-3779-df58afdb52ad
begin
    scatter(generated_data.x, generated_data.y, label = "training data", xrange = (-3.3, 3.3), yrange = (-9, 7), legend = :topleft)
    plot_residuals!(generated_data.x, generated_data.y, fitted_linear_func(mach2))
    plot!(fitted_linear_func(mach2), color = :green, label = "fit")
end

# ╔═╡ f63c0600-eefe-11eb-0023-9d5e5ddd987e
(training_loss = L(mach2, generated_data),
 test_loss = L(mach2, generated_test_data))

# ╔═╡ 0e6efce9-093a-4758-a649-6cff525711a5
md"# Distributions

Use the sliders below to get a feeling for how the distributions depend on ``f(x)``.
In the plots, the triangle indicates the current value of ``f(x)``. Note that ``f(x)`` can be an arbitrarily complicated function of ``x``. What matters is that the output of this function is used to shape the probability distribution.
"

# ╔═╡ fbc70eaa-df15-423a-9885-93a5fa27fbc5
begin
	normalpdf(m) = x -> pdf(Normal(m, 1), x)
	logistic(x) = 1 / (1 + exp(-x))
	softmax(x) = exp.(x)/sum(exp.(x))
end;

# ╔═╡ 4b3b509f-7995-46b3-9130-8d881ee1e0ae
md"## Normal

``f(x)`` = $(@bind ŷₙ Slider(-5:.1:5, show_value = true))
"

# ╔═╡ 2adab164-1b4f-4f2f-ade3-dd70544c692e
begin
    plot(normalpdf(ŷₙ), xrange = (-8, 8), label = nothing,
         yrange = (0, .41), xlabel = "y", ylabel = "probability density p(y|x)", size = (350, 250))
    scatter!([ŷₙ], [0], markersize = 8, shape = :utriangle, label = nothing)
end

# ╔═╡ 8060861a-1931-4fae-b6ee-afb94ef2a3d5
md"## Bernoulli

f(x) = $(@bind ŷₛ Slider(-5:.1:5, show_value = true))
"

# ╔═╡ 6124e1eb-0aec-4ac9-aee6-383b25865220
let
	p1 = bar([0, 1], [logistic(ŷₛ), logistic(-ŷₛ)],
               xtick = ([0, 1], ["A", "B"]), ylabel = "probability P(class|x)",
               xlabel = "class", ylim = (0, 1), xlim = (-.5, 1.5), label = nothing)
    p2 = plot(logistic, xlim = (-5, 5), ylim = (0, 1),
               xlabel = "f(x)", ylabel = "probability of A", label = nothing)
    scatter!([ŷₛ], [0], markersize = 8, shape = :utriangle, label = nothing)
    vline!([ŷₛ], linestyle = :dash, color = :black)
    hline!([logistic(ŷₛ)], linestyle = :dash, color = :black)
    plot(p1, p2, size = (700, 350), legend = false)
end

# ╔═╡ 5d3efbe5-8937-45c0-894d-d24ece9b2262
md"## Categorical

f₁(x) = $(@bind ŷ₁ Slider(-5:.1:5, show_value = true))

f₂(x) = $(@bind ŷ₂ Slider(-5:.1:5, show_value = true))

f₃(x) = $(@bind ŷ₃ Slider(-5:.1:5, show_value = true))

f₄(x) = $(@bind ŷ₄ Slider(-5:.1:5, show_value = true))
"

# ╔═╡ 897fc69d-f1e9-49e9-9a61-25eb6849f7ec
bar(1:4, softmax([ŷ₁, ŷ₂, ŷ₃, ŷ₄]), ylabel = "probability P(class|x)",
    ylim = (0, 1), xlim = (.5, 4.5), xlabel = "class", legend = false,
    xtick = (1:4, ["A", "B", "C", "D"]), size = (350, 250))

# ╔═╡ f63c0616-eefe-11eb-268a-271fdc2ddd98
md"# Logistic Regression (Linear Classification)

## Blackboard Example"

# ╔═╡ f63c0616-eefe-11eb-1cea-dfdaa64e6233
classification_data = DataFrame(x = [0., 2., 3.],
	                            y = categorical(["B", "A", "B"],
									            levels = ["B", "A"],
								                ordered = true))

# ╔═╡ 7dfa42e3-c409-43a7-9952-a64fbad63c7f
md"In the cell above we use the function `categorical` to tell the computer that the strings \"A\" and \"B\" indicate membership in different categories."

# ╔═╡ f63c061e-eefe-11eb-095b-8be221b33d49
mach3 = machine(LogisticClassifier(penalty = :none), # model
                select(classification_data, :x),     # input
                classification_data.y);              # output

# ╔═╡ 9a208086-1fcd-4520-b566-ea295e934d75
fit!(mach3);

# ╔═╡ 1bbbb7fe-494b-4692-aaa4-e39599851327
fitted_params(mach3)

# ╔═╡ 38ccf735-5195-4e00-868f-95a895c05985
md"In the following cell we define the loss function to compute the loss of the optimal parameters (found by the `fit!`)."

# ╔═╡ d0c4804f-c66d-4405-a7dc-1e85974e261f
ll(θ) = log(logistic(-θ[1])) +
        log(logistic(θ[1] + 2θ[2])) +
        log(logistic(-θ[1] - 3θ[2]))

# ╔═╡ b8971599-494a-4fb3-9fa0-fb374fb7d79d
ll([-1.28858, .338548])

# ╔═╡ b8b81c1b-0faf-4ce9-b690-2f6cc9542b0f
md"The likelihood of the optimal parameters is approximately -1.9. We could have obtained the same result using the `MLJ` function `log_loss`. This function computes the negative log-likelihood for each individual data point. To get the total likelihood we need to take the negative of the sum of these values."

# ╔═╡ d034a44b-e331-4929-9053-351e7fe9aa94
-sum(log_loss(predict(mach3), classification_data.y))

# ╔═╡ 0e775dfb-0da4-4536-886c-ada8c176a073
md"For the `LogisticClassifier` the `predict` method returns the conditional probabilities of the classes. Click on the little gray triangle below to toggle the display of the full output."

# ╔═╡ f63c061e-eefe-11eb-3b91-7136b4a16616
p̂ = predict(mach3, DataFrame(x = -1:.5:2))

# ╔═╡ 0c90f5b6-8a3b-41d8-9f51-d7d7c6b06ba0
md"If we want to extract the probability of a given response, we can use the `pdf` function."

# ╔═╡ 5224d406-4e02-424d-9502-a22e0614cb96
pdf.(p̂, "A")

# ╔═╡ 34c49e49-a5e7-48ad-807b-c0624a59a367
md"If we want to get as a response the class with the highest probability, we can use the function `predict_mode`."

# ╔═╡ f63c0628-eefe-11eb-3125-077e533456d9
predict_mode(mach3, DataFrame(x = -1:.5:2))

# ╔═╡ f63c0632-eefe-11eb-3b93-8549af7aaed9
md"""# Exercises

## Conceptual
1. We have some training data ``((x_1 = 0, y_1 = -1), (x_2 = 2, y_2 = 4),
   (x_3 = 2, y_3 = 3))`` and a family of probability densities
   ``p(y|x) = \frac1{\sqrt{2\pi\sigma^2}}\exp\left(-\frac{(y - \theta_0 - \theta_1 x)^2}{2\sigma^2}\right)``.
   - Write the log-likelihood function of the parameters ``\theta_0`` and ``\theta_1``
     for this data and model.
   - Find the parameters ``\hat\theta_0`` and ``\hat\theta_1`` that maximize the log-likelihood function and compare your result to the solution we found in our blackboard example of linear regression as a loss minimizing machine.
   - Show for general training data ``((x_1, y_1), (x_2, y_2), \ldots, (x_n, y_n))`` that the log-likelihood function is maximized by the same ``\hat\theta`` that minimizes the loss function of linear regression.
1. Explain whether each scenario is a classification or regression problem, and
   indicate whether we are most interested in inference/interpretation or prediction.
   Finally, provide ``n`` and ``p``.
   - We collect a set of data on the top 500 firms in the US. For each firm we
     record profit, number of employees, industry and the CEO salary. We are
     interested in understanding which factors affect CEO salary.
   - We are considering launching a new product and wish to know
     whether it will be a success or a failure. We collect data on 20
     similar products that were previously launched. For each product
     we have recorded whether it was a success or failure, price
     charged for the product, marketing budget, competition price,
     and ten other variables.
   - We are interested in predicting the % change in the USD/Euro
     exchange rate in relation to the weekly changes in the world
     stock markets. Hence we collect weekly data for all of 2012. For
     each week we record the % change in the USD/Euro, the %
     change in the US market, the % change in the British market,
     and the % change in the German market.
1. Take one of the examples of supervised machine learning applications
   discussed in the introductory lecture of this course. Discuss with a colleague
   - the data generating process for X
   - the data generating process for Y|X
   - where the noise comes from
   - a distribution that could be used to model Y|X

## Applied
1. Change the noise level ``\sigma``, the size of the training data ``n`` and the seed with the sliders of the section "Linear Regression" and observe the training and the test losses. Write down your observations when
    - ``n`` is small.
    - ``n`` is large.
    - Compare training loss and test loss when ``n`` is large for different seed values.
    - Compare for large ``n`` the test error to ``\sigma^2`` for different values of ``\sigma``.
1. Write a data generator function that samples inputs ``x`` normally distributed
   with mean 2 and standard deviation 3. The response ``y\in\{\mbox{true, false}\}``
   should be sampled from a Bernoulli distribution with rate of ``\mbox{true}``
   equal to ``\sigma(0.5x - 2.7)`` where ``\sigma(x) = 1/(1 + e^{-x})`` is the sigmoid (or logistic) function.
   - Create a training set of size ``n = 20``.
   - Fit the data with logistic regression.
   - Look at the fitted parameters.
   - Predict the probability of class `true` on the training input
   - Determine the class with highest predicted probability and compare the result to the labels of the training set.
   - Create a test set of size ``n = 10^4`` where the input is always at ``x = 4``. Estimate the average test error at ``x = 4`` using this test set. Use the negative log-likelihood as error function.
   - Compute the test error at ``x = 4`` using the fitted parameters and compare your result to the previous result. *Hint:* Have a look at the slides for how to compute the test error when the parameters of the generator and the fitted function are known.
   - Rerun your solution with different training sets of size ``n = 20`` and write down your observations.
   - Rerun your solution multiple times with training sets of size ``n = 10^4``. Compare the fitted parameters to the one of the generator and look at the test error. Write down your observations.
"""

# ╔═╡ e3aa458a-486c-4d5b-a330-67fb68e6f516
MLCourse.footer()

# ╔═╡ Cell order:
# ╟─f63c04dc-eefe-11eb-1e24-1d02a686920a
# ╠═f63c04e8-eefe-11eb-39f6-83b31ebe73e7
# ╠═f63c04e8-eefe-11eb-1a14-8305504a6f1c
# ╟─f63c04f2-eefe-11eb-2562-e530426e4300
# ╠═f63c04f2-eefe-11eb-3ceb-ff1e36a2a302
# ╠═f63c04fc-eefe-11eb-35b6-5345dda134e7
# ╟─f63c04fc-eefe-11eb-043c-7fec2fc41913
# ╠═f63c0506-eefe-11eb-1857-e3eaf731c152
# ╟─c82fb649-e60d-43fc-b158-444b27b6e9cf
# ╠═f63c0506-eefe-11eb-399f-455ee06548da
# ╟─f63c050e-eefe-11eb-1709-0ddce2aee0ed
# ╠═f63c051a-eefe-11eb-0db2-c519292f88a2
# ╠═f63c051a-eefe-11eb-16eb-a572fe41aa43
# ╠═f63c0524-eefe-11eb-3abd-63d677b12db9
# ╟─f63c0524-eefe-11eb-3732-258d6989e217
# ╠═f63c052e-eefe-11eb-3a14-e5e8f3d578a8
# ╠═f63c052e-eefe-11eb-0884-7bd433ce0c5e
# ╠═f63c0538-eefe-11eb-2a4b-d19efe6b6689
# ╠═f63c0538-eefe-11eb-2eca-2d8bf12fae95
# ╟─f63c0538-eefe-11eb-2808-0f32a2fa84cf
# ╠═f63c0540-eefe-11eb-01d6-454a05f9182a
# ╠═f63c054c-eefe-11eb-13db-171084b417a9
# ╟─f63c0560-eefe-11eb-0821-73afecf00357
# ╠═f63c0560-eefe-11eb-0251-2d0443448e66
# ╠═f63c056a-eefe-11eb-199e-2b05b1385be4
# ╠═f63c056a-eefe-11eb-218c-afc51fdf1704
# ╟─f63c0574-eefe-11eb-2d7c-bd74b7d83f23
# ╟─f63c057e-eefe-11eb-1672-19db9efd0d32
# ╟─f63c057e-eefe-11eb-0dd2-3333b4e2d220
# ╟─f63c0588-eefe-11eb-0687-a7c44f9177a4
# ╟─f63c0588-eefe-11eb-1b0b-db185fcbdc4e
# ╟─f63c05a6-eefe-11eb-0ef4-df7e458e9ec6
# ╠═f63c05a6-eefe-11eb-2927-b1e9dbbe032d
# ╠═f63c05b2-eefe-11eb-03cd-ef71f9055eb3
# ╟─f63c05b2-eefe-11eb-21dc-1f79ae779316
# ╠═f63c0592-eefe-11eb-0d69-91d4d396c540
# ╠═f63c05ba-eefe-11eb-18b5-7522b326ab65
# ╠═f63c05ba-eefe-11eb-1b7a-21cbd9bbeb37
# ╠═f63c05c4-eefe-11eb-030e-6f7484849852
# ╟─25eb24a4-5414-481b-88d5-5ad50d75f8c9
# ╠═bd9a7214-f39e-4ae5-995a-6ec02d613fda
# ╠═f0f6ac8f-2e61-4b08-9383-45ff5b02c5b1
# ╟─7336e7b4-a7c2-4409-8d9d-00285cf633fb
# ╠═f63c05c4-eefe-11eb-2d60-f7bf7beb8118
# ╠═f63c05ce-eefe-11eb-1ccd-c9b7714f3bd5
# ╟─8fe94709-3673-44cc-8702-83bd7e2cad51
# ╠═27f44b28-f6d6-4e1f-8dbe-6fe7789b9a18
# ╟─43b28121-5d07-4400-8710-287de7b978a4
# ╠═f63c05d8-eefe-11eb-2a11-fd8f954bf059
# ╠═f63c05d8-eefe-11eb-0a2e-970192b02d61
# ╟─f63c05e4-eefe-11eb-012a-9bae1d87f2b5
# ╠═f63c0592-eefe-11eb-2a76-15de55eff3ad
# ╟─ca8e60e6-1d12-4076-8684-9a7c9e5ac509
# ╠═4987ad0c-d100-4a5d-a725-12f23bcc7aa3
# ╠═93f529ab-a386-4746-94de-63a536b5e2aa
# ╟─673a773c-2f6b-4676-8631-49ba08ec28a7
# ╠═f63c05e4-eefe-11eb-0848-9b4eb61d920c
# ╟─f63c059c-eefe-11eb-208c-d5781711deb7
# ╠═f63c05ec-eefe-11eb-0a92-215711928d02
# ╟─f63c0600-eefe-11eb-3779-df58afdb52ad
# ╠═f63c0600-eefe-11eb-0023-9d5e5ddd987e
# ╟─0e6efce9-093a-4758-a649-6cff525711a5
# ╠═36a44a52-8a29-4f92-8136-2c66a81a1f5f
# ╠═fbc70eaa-df15-423a-9885-93a5fa27fbc5
# ╟─4b3b509f-7995-46b3-9130-8d881ee1e0ae
# ╟─2adab164-1b4f-4f2f-ade3-dd70544c692e
# ╟─8060861a-1931-4fae-b6ee-afb94ef2a3d5
# ╟─6124e1eb-0aec-4ac9-aee6-383b25865220
# ╟─5d3efbe5-8937-45c0-894d-d24ece9b2262
# ╟─897fc69d-f1e9-49e9-9a61-25eb6849f7ec
# ╟─f63c0616-eefe-11eb-268a-271fdc2ddd98
# ╠═f63c0616-eefe-11eb-1cea-dfdaa64e6233
# ╟─7dfa42e3-c409-43a7-9952-a64fbad63c7f
# ╠═f63c061e-eefe-11eb-095b-8be221b33d49
# ╠═9a208086-1fcd-4520-b566-ea295e934d75
# ╠═1bbbb7fe-494b-4692-aaa4-e39599851327
# ╟─38ccf735-5195-4e00-868f-95a895c05985
# ╠═d0c4804f-c66d-4405-a7dc-1e85974e261f
# ╠═b8971599-494a-4fb3-9fa0-fb374fb7d79d
# ╟─b8b81c1b-0faf-4ce9-b690-2f6cc9542b0f
# ╠═d034a44b-e331-4929-9053-351e7fe9aa94
# ╟─0e775dfb-0da4-4536-886c-ada8c176a073
# ╠═f63c061e-eefe-11eb-3b91-7136b4a16616
# ╟─0c90f5b6-8a3b-41d8-9f51-d7d7c6b06ba0
# ╠═5224d406-4e02-424d-9502-a22e0614cb96
# ╟─34c49e49-a5e7-48ad-807b-c0624a59a367
# ╠═f63c0628-eefe-11eb-3125-077e533456d9
# ╟─f63c0632-eefe-11eb-3b93-8549af7aaed9
# ╟─0cfbd543-8b9b-406e-b3b4-c6cafbbec212
# ╟─f63c04d4-eefe-11eb-3fda-1729ac6de2cb
# ╟─e3aa458a-486c-4d5b-a330-67fb68e6f516

max_stars_repo_name: baltisberger/MLCourse

max_issues_repo_name: baltisberger/MLCourse

max_forks_repo_name: baltisberger/MLCourse


text:
using SumTrees
if VERSION < v"0.5-"
    using BaseTestNext
else
    using Base.Test
end

@testset "Tiny" begin
    @testset "Default Type" begin
        tree = SumTree(1)
        @test isa(tree, SumTree)
        @test isa(tree, SumTree{Float64, 1})
        @test eltype(tree) === Float64

        @test tree[1] == 0.0
        @test sum(tree) == 0.0

        tree[1] = 3.0
        @test tree[1] == 3.0
        @test sum(tree) == 3.0

        tree[1] = 8
        @test tree[1] == 8.0
        @test sum(tree) == 8.0
        @test isa(tree[1], Float64)

        @test_throws BoundsError tree[0]
        @test_throws BoundsError tree[2]

        @test_throws BoundsError tree[0] = 1.0
        @test_throws BoundsError tree[2] = 1.0

        @test [x for x in tree] == [8.0]
    end

    @testset "Rational" begin
        tree = SumTree(Rational, 1)
        @test isa(tree, SumTree)
        @test isa(tree, SumTree{Rational, 1})

        @test tree[1] == zero(Rational)
        @test sum(tree) == zero(Rational)

        tree[1] = 3//4
        @test tree[1] == 3//4
        @test sum(tree) == 3//4
        @test isa(tree[1], Rational)

        @test [x for x in tree] == [3//4]
    end
end

@testset "Normal" begin
    @testset "Default Type" begin
        tree = SumTree(8)

        @test collect(tree) == zeros(Float64, 8)

        buffer = IOBuffer()
        show(buffer, tree)

        @test takebuf_string(buffer) == "SumTree{Float64, 8} ::\n    (0.0,)\n    (0.0,0.0,)\n    (0.0,0.0,0.0,0.0,)\n    (0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,)\n"

        tree[1] = 3
        tree[4] = 99

        @test tree[1] == 3
        @test tree[4] == 99
        @test tree[2] == 0
        @test sum(tree) == 102

        tree[8] = 0.0001
        @test tree[8] == 0.0001
        @test tree[1] == 3
        @test tree[4] == 99
        @test tree[2] == 0
        @test sum(tree) == 102.0001
    end
end


text:
function create(c::Client, ::Type{Overwrite}, ch::DiscordChannel, r::Role; kwargs...)
    return edit_channel_permissions(c, ch.id, r.id; kwargs..., type=OT_ROLE)
end
function create(c::Client, ::Type{Overwrite}, ch::DiscordChannel, u::User; kwargs...)
    return edit_channel_permissions(c, ch.id, u.id; kwargs..., type=OT_MEMBER)
end

function update(c::Client, o::Overwrite, ch::DiscordChannel; kwargs...)
    return edit_channel_permissions(c, ch.id, o.id; kwargs..., type=o.type)
end

function delete(c::Client, o::Overwrite, ch::DiscordChannel)
    return delete_channel_permission(c, ch.id, o.id)
end


text:
using Plots
using Random
Random.seed!(1)

# Differential Dynamic Programming
include_ddp()
include(joinpath(pwd(), "src/differential_dynamic_programming/multiple_models.jl"))

include_model("acrobot")

function M(model::Acrobot, x, w)
    a = (model.J1 + model.J2 + model.m2 * model.l1 * model.l1
         + 2.0 * model.m2 * model.l1 * model.lc2 * cos(x[2]))

    b = model.J2 + model.m2 * model.l1 * model.lc2 * cos(x[2])

    c = model.J2

    @SMatrix [a b;
              b c]
end

function τ(model::Acrobot, x, w)
    a = (-1.0 * model.m1 * model.g * model.lc1 * sin(x[1])
         - model.m2 * model.g * (model.l1 * sin(x[1])
         + model.lc2 * sin(x[1] + x[2])))

    b = -1.0 * model.m2 * model.g * model.lc2 * sin(x[1] + x[2])

    @SVector [a,
              b]
end

function C(model::Acrobot, x, w)
    a = -2.0 * model.m2 * model.l1 * model.lc2 * sin(x[2]) * x[4]
    b = -1.0 * model.m2 * model.l1 * model.lc2 * sin(x[2]) * x[4]
    c = model.m2 * model.l1 * model.lc2 * sin(x[2]) * x[3]
    d = 0.0

    @SMatrix [a b;
              c d]
end

function B(model::Acrobot, x, w)
    @SMatrix [0.0;
              1.0]
end

function f(model::Acrobot, x, u, w)
    q = view(x, 1:2)
    v = view(x, 3:4)
    qdd = M(model, q, w) \ (-1.0 * C(model, x, w) * v
            + τ(model, q, w) + B(model, q, w) * u[1:1] - [model.b1 + w[1]; model.b2 + w[2]] .* v)
    @SVector [x[3],
              x[4],
              qdd[1],
              qdd[2]]
end
state_output(model::Acrobot, x) = x

n, m, d = 4, 1, 2
model = Acrobot{Midpoint, FixedTime}(n, m, d, 1.0, 0.33, 1.0, 0.5, 1.0, 0.33, 1.0, 0.5, 9.81, 0.1, 0.1)

n = model.n
m = model.m
d = model.d

# Models
N = 100#2 * n + 1

# Time
T = 101
h = 0.05
tf = h * (T - 1)
t = range(0, stop = tf, length = T)

# Initial conditions, controls, disturbances
# x1 = [p_ref[1]; 0.0]
x1 = [0.0, 0.0, 0.0, 0.0]
xT = [π, 0.0, 0.0, 0.0] # goal state
_ū = [1.0e-1 * rand(model.m) for t = 1:T-1]
_ū = u_ref
ū = [_ū for i = 1:N]

W = Distributions.MvNormal(zeros(model.d), Diagonal([1.0e-8, 1.0e-8]))
wi = [vec(rand(W, 1)) for i = 1:N]
w = [[zeros(model.d) for t = 1:T-1], [[wi[i] for t = 1:T-1] for i = 1:N-1]...]
# w = [[vec(rand(W, 1)) for t = 1:T-1] for i = 1:N]
# w = [[zeros(model.d) for t = 1:T-1] for i = 1:N]

# Rollout
x̄ = [rollout(model, x1, ū[1], w[i], h, T) for i = 1:N]

p_idx = (1:2)
plt = plot()
for i = 1:N
	plt = plot!(t, hcat(x̄[i]...)[p_idx, :]',
		label = "")
end
display(plt)

# Objective
Q = [(t < T ? Diagonal(1.0e-3 * ones(model.n))
    : Diagonal(100.0 * ones(model.n))) for t = 1:T]
R = [Diagonal(1.0e-5 * ones(model.m)) for t = 1:T-1]
obj = StageCosts([QuadraticCost(Q[t], nothing,
	t < T ? R[t] : nothing, nothing) for t = 1:T], T)

function g(obj::StageCosts, x, u, t)
    T = obj.T
    if t < T
		Q = obj.cost[t].Q
		R = obj.cost[t].R
        return (x - xT)' * Q * (x - xT) + u' * R * u
    elseif t == T
		Q = obj.cost[t].Q
        return (x - xT)' * Q * (x - xT)
    else
        return 0.0
    end
end

models_data = models(model, obj, deepcopy(x̄), deepcopy(ū), w, h, T)

prob = problem_data(models_data)

function ddp_solve!(prob::ProblemData;
    max_iter = 100,
    grad_tol = 1.0e-5,
    verbose = true)

	println()
    verbose && println("Differential Dynamic Programming")

	# data
	p_data = prob.p_data
	m_data = prob.m_data
	s_data = prob.s_data

    # compute objective
    s_data.obj = objective(m_data, mode = :nominal)

	N = length(prob.m_data)

    for i = 1:max_iter
        # derivatives
        derivatives!(m_data)

        # backward pass
        backward_pass!(p_data, m_data)

		# # set nominal trajectories
		# x_ref = [sum([m.x̄[t] for m in prob.m_data]) ./ N for t = 1:T]
		# u_ref = [sum([m.ū[t] for m in prob.m_data]) ./ N for t = 1:T-1]
		#
		# # disturbances
		# w = [[rand(W, 1) for t = 1:T-1] for i = 1:N]
		#
		# for i = 1:N
		# 	# m_data[i].x̄ .= deepcopy(x_ref)
		# 	# m_data[i].ū .= deepcopy(u_ref)
		# 	m_data[i].w .= deepcopy(w[i])
		# end

		if i > 1
			for i = 2:N
				m_data[i].x̄ .= deepcopy(m_data[1].x̄)
				m_data[i].ū .= deepcopy(m_data[1].ū)
				# m_data[i].w .= deepcopy(w[i])
			end
		end
        # forward pass
        forward_pass!(p_data, m_data, s_data)

		# if i > 1
		# 	for i = 2:N
		# 		m_data[i].x̄ .= deepcopy(m_data[1].x̄)
		# 		m_data[i].ū .= deepcopy(m_data[1].ū)
		# 		# m_data[i].w .= deepcopy(w[i])
		# 	end
		# end
		# # set nominal trajectories
		# x_ref = [sum([m.x̄[t] for m in prob.m_data]) ./ N for t = 1:T]
		# u_ref = [sum([m.ū[t] for m in prob.m_data]) ./ N for t = 1:T-1]

        # check convergence
        grad_norm = norm(s_data.gradient, Inf)
        verbose && println("     iter: $i
             cost: $(s_data.obj)
			 grad_norm: $(grad_norm)")
		grad_norm < grad_tol && break
        !s_data.status && break
    end
end

objective(prob.m_data)
prob.m_data[1].x̄

# Solve
@time ddp_solve!(prob,
    max_iter = 250, verbose = true,
	grad_tol = 1.0e-8)

x = [m.x for m in prob.m_data]
u = [m.u for m in prob.m_data]

x̄ = [m.x̄ for m in prob.m_data]
ū = [m.ū for m in prob.m_data]

x_ref = deepcopy(prob.m_data[1].x̄)#[sum([m.x̄[t] for m in prob.m_data]) ./ N for t = 1:T]
u_ref = deepcopy(prob.m_data[1].ū)#[sum([m.ū[t] for m in prob.m_data]) ./ N for t = 1:T-1]

# Visualize
idx = (1:2)
_plt = plot(hcat([[p_ref[t]; 0.0][idx] for t = 1:T]...)',
    width = 2.0, color = :black, label = "")
for i = 1:N
	_plt = plot!(hcat(x[i]...)[idx, :]', color = :magenta, label = "")
end
_plt = plot!(hcat(x_ref...)[idx, :]', color = :orange, width = 5.0, label = "")
display(_plt)
# plot(hcat(u..., u[end])', linetype = :steppost)

# Simulate policy
include_dpo()
include(joinpath(pwd(), "examples/direct_policy_optimization/simulate.jl"))

# Model
model_sim = model
x1_sim = copy(x1)
T_sim = 10 * (T-1) + 1

# Disturbance distributions

# Time
t_sim = range(0, stop = tf, length = T_sim)
dt_sim = tf / (T_sim - 1)

# Policy
policy = linear_feedback(model.n, model.m)
# K, P = tvlqr(model, x̄, ū, h, [Q for t = 1:T], [R for t = 1:T-1])
K = [-K for K in prob.p_data.K]

# Simulate
N_sim = 100
x_sim = []
u_sim = []
J_sim = []
Random.seed!(1)
W_sim = W#Distributions.MvNormal(zeros(model.d), Diagonal([0.0, 0.0, 1.0e-1]))
w_sim = rand(W_sim, T_sim)
for i = 1:N_sim
	println("sim: $i")
	wi = vec(rand(W_sim, 1))
	w_sim = hcat([wi for t = 1:T_sim]...)

	x_ddp, u_ddp, J_ddp, Jx_ddp, Ju_ddp = simulate(
		model_sim,
		linear_feedback(model.n, model.m), K,
	    x_ref, u_ref,
		Q, R,
		T_sim, h,
		x1_sim,
		w_sim)

	push!(x_sim, x_ddp)
	push!(u_sim, u_ddp)
	push!(J_sim, J_ddp)
end

# Visualize
idx = (1:2)
_plt = plot(t, hcat(x_ref...)[idx, :]',
    width = 2.0, color = :black, label = "")

for (i, xi) in enumerate(x_sim)
	_plt = plot!(t_sim, hcat(xi...)[idx, :]',
    	width = 1.0, color = :magenta, label = i == 1 ? "sim" : "")
end
_plt = plot!(t, hcat(x_ref...)[idx, :]',
    width = 2.0, color = :black, label = "ref")
_plt = plot!(ylims = (-2, π + 2.0); xlabel = "time (s)",
	ylabel = "state",
	title = "N_mc = $N, N_sim = $(N_sim), J_avg = $(round(mean(J_sim), digits = 3))",
	legend = :bottomright)
display(_plt)

# _plt = plot(t, hcat(u_ref..., u_ref[end])[1:1, :]',
#     width = 2.0, color = :black, label = "", linetype = :steppost)
# for (i, ui) in enumerate(u_sim)
# 	_plt = plot!(t_sim, hcat(ui..., ui[end])[1:1, :]',
#     	width = 1.0, color = :magenta, label = i == 1 ? "sim" : "")
# end
# _plt = plot!(t, hcat(u_ref..., u_ref[end])[1:1, :]',
#     width = 2.0, color = :black, label = "ref", linetype = :steppost)
# _plt = plot!(ylims = (-20.0, 20.0); xlabel = "time (s)",
# 	ylabel = "control",
# 	title = "N_mc = $N, N_sim = $(N_sim), J_avg = $(round(mean(J_sim), digits = 2))",
# 	linetype = :steppost)
# display(_plt)
@show mean(J_sim) # 0.00157
round(mean(J_sim), digits = 3)


text:
using DynamicalSystemsBase
using Test

ti = time()

@testset "DynamicalSystemsBase tests" begin
include("dynsys_integrators.jl")
include("dynsys_tangent.jl")
include("dynsys_inference.jl")
include("continuous_systems.jl")
include("discrete_systems.jl")
include("integrators_with_callbacks.jl")
include("norm_tests.jl")
include("projected_integrator_tests.jl")
end

ti = time() - ti
println("\nTest took total time of:")
println(round(ti, digits=3), " seconds or ", round(ti/60, digits=3), " minutes")


text:
# plot kernek density estimation

module MarginalKde200
    using DataFrames, CSV, StatsPlots
    pyplot()

    function main(is_test=false)
        # input
        data_path = joinpath(split(@__FILE__, "src")[1], "data/sensor_data_200.txt")
        df = CSV.read(data_path, DataFrame,
                      header=["date", "time", "ir", "lidar"],
                      delim=' ')

        # plot marginal kde
        marginalkde(df.ir, df.lidar, c=:ice,
                    xlabel="ir", ylabel="liDAR")

        if is_test == false
            save_path = joinpath(split(@__FILE__, "src")[1], "src/prob_stats/multi_dim_gauss_dist/marginal_kde/marginal_kde_200.png")
            savefig(save_path)
        end
    end
end

max_stars_repo_name: ShisatoYano/JuliaAutonomy

max_issues_repo_name: ShisatoYano/JuliaAutonomy

max_forks_repo_name: ShisatoYano/JuliaAutonomy


text:
function perform_conjugate_gradient(A,y)

    """
     perform_conjugate_gradient - perform (bi)-conjugate gradient

         x = perform_conjugate_gradient(A,y);

       Solves for A*x=y.
       Works for vector x,y

       Important: the algorithm assumes that the matrix is symmetric definite
       positive.

       Copyright (c) 2007 Gabriel Peyre
    """

    niter = 100
    epsilon = 1e-5
    is_sdp = 1
    x = zeros(size(A,2),1)
    normb = epsilon
    r  = y - A*x
    p = r
    r0 = sum(vec(r).^2)

    err = [sum(r0)]
    for it in 1:niter
        # auxiliary vector
        w  = A*p

        d = sum(vec(p) .* vec(w));

        if abs(d)<eps()
            d=1
        end
        alpha = repeat( [r0 / d], outer=(size(x,1), 1) );           # found optimal alpha in line search
        x = x + alpha.*p                       # new guess
        r = r - alpha.*w                       # the residual is in fact r=b-A*x

        rprev = r0;                             # save norm of the old residual
        r0 = sum(r.^2);                         # compute norm of new residual

        append!(err, sqrt( sum(r.^2) ))

        if err[end]<normb
            break
        end

        # search direction
        beta = r0./rprev;
        p = r + repeat([beta], outer=(size(x,1), 1)).*p

    end
    return x
end


text:
# This file should contain all recipes for plotting correlators, effective masses, volume plots, fits and channel-spectra (something else?)
# Automatically creating titles depending on settings in theory.jl (use Plots.default() for this)


text:
# This file is part of JuliaFMI.
# Licensed under MIT: https://github.com/AnHeuermann/JuliaFMI/blob/master/LICENSE.txt

# This file contains type definitions and constructores for those types.

using CBinding

"""
Declaration of FMI2 Types
"""

@cenum FMUType {
    modelExchange,
    coSimulation
}

@cenum ModelState {
    modelUninstantiated,
    modelInstantiated,
    modelInitializationMode,
    modelEventMode,
    modelContinuousTimeMode,
    modelStepComplete,
    modelStepProgress,
    modelStepFailed,
    modelStepCancelled,
    modelTerminated,
    modelError,
    modelFatal
}

@cenum FMI2Status {
    fmi2OK,
    fmi2Warning,
    fmi2Discard,
    fmi2Error,
    fmi2Fatal,
    fmi2Pending
}

@ctypedef DynamicPointersCTypedef @cstruct DynamicPointers {
    realValues::Ptr{Nothing} # Vector of Float64
    instanceName::Ptr{Nothing} # Vector of Char
    eventIndicators::Ptr{Nothing} # Vector of Cint
    # A lookup that checks if the inputted index (corresponding to a category)
    # is checked for logging.
    categoriesToLogLookup::Ptr{Nothing} # Vector of Cint
}
DynamicPointers() = DynamicPointers(zero)

@ctypedef CallbackFunctionsCTypedef @cstruct CallbackFunctions {
    logger::Ptr{Nothing}
    allocateMemory::Ptr{Nothing}
    freeMemory::Ptr{Nothing}
    stepFinished::Ptr{Nothing}
    componentEnvironment::Ptr{Nothing}
}

function CallbackFunctions(;stepFinished::Ptr=C_NULL, componentEnvironment::Ptr=C_NULL)

    fmi2CallbacLogger_Cfunc = dlsym(dlopen(@libLogger), :logger)
    fmi2AllocateMemory_funcWrapC = @cfunction(fmi2AllocateMemory, Ptr{Cvoid}, (Csize_t, Csize_t))
    fmi2FreeMemory_funcWrapC = @cfunction(fmi2FreeMemory, Cvoid, (Ptr{Cvoid},))

    return CallbackFunctions(
        fmi2CallbacLogger_Cfunc,
        fmi2AllocateMemory_funcWrapC,
        fmi2FreeMemory_funcWrapC,
        stepFinished,
        componentEnvironment)
end

function CallbackFunctions(a::Ptr, b::Ptr, c::Ptr, d::Ptr, e::Ptr)
    functions = CallbackFunctions(zero)
    functions.logger = a
    functions.allocateMemory = b
    functions.freeMemory = c
    functions.stepFinished = d
    functions.componentEnvironment = e
    return functions
end

@ctypedef EventInfoCTypedef @cstruct EventInfo {
    newDiscreteStatesNeeded::Cint
    terminateSimulation::Cint
    nominalsOfContinuousStatesChanged::Cint
    valuesOfContinuousStatesChanged::Cint
    nextEventTimeDefined::Cint
    nextEventTime::Cdouble
}
function EventInfo()
    ceventinfo = EventInfo(zero)
    ceventinfo.newDiscreteStatesNeeded = 1
    ceventinfo.terminateSimulation = 1
    ceventinfo.nominalsOfContinuousStatesChanged = 1
    ceventinfo.valuesOfContinuousStatesChanged = 1
    ceventinfo.nextEventTimeDefined = 1
    ceventinfo.nextEventTime = -1.0
    return ceventinfo
end

## Forward reference to FMI2Component
## Pointed to by the fmi2Component void pointer.
@ctypedef FMI2ComponentCType @cstruct FMI2Component {

    functions::Cconst{Ptr{CallbackFunctions}}

    (dynPtrs)::@cstruct {
        realValues::Ptr{Cdouble} # Vector of Float64
        instanceName::Ptr{Cchar} # Vector of Char
        eventIndicators::Ptr{Cint} # Vector of Cint
        # A lookup that checks if the inputted index (corresponding to a category)
        # is checked for logging.
        categoriesToLogLookup::Ptr{Cint} # Vector of Cint
    }

    eventInfo::EventInfo

    # Caches the status of the do step routines

    doStepStatus::FMI2Status;

    time::Cdouble
    startTime::Cdouble

    stopTime::Cdouble
    stopTimeDefinedFlag::Cint

    tolerance::Cdouble
    toleranceDefinedFlag::Cint

    computeVarFlag::Cint
    loggingOnFlag::Cint

    state::ModelState
    type::FMUType
}

@cenum NamingConvention {
    flat
    structured
}
function str2NamingConvention(inStr::String)
    if inStr =="flat"
        return flat
    elseif inStr == "structured"
        return structured
    else
        error("Can not convert String \"$in\" to NamingConvention.")
    end
end

max_stars_repo_name: AnHeuermann/JuliaFMU

max_issues_repo_name: AnHeuermann/JuliaFMU

max_forks_repo_name: jwscook/JuliaFMI


text:
# Not in DiffEqDiffTools because `u` -> scalar isn't used anywhere else,
# but could be upstreamed.
mutable struct UGradientWrapper{fType,tType,P} <: Function
  f::fType
  t::tType
  p::P
end

(ff::UGradientWrapper)(uprev) = ff.f(uprev,ff.p,ff.t)

Base.@pure function determine_chunksize(u,alg::DiffEqBase.AbstractSensitivityAlgorithm)
  determine_chunksize(u,get_chunksize(alg))
end

Base.@pure function determine_chunksize(u,CS)
  if CS != 0
    return CS
  else
    return ForwardDiff.pickchunksize(length(u))
  end
end

function jacobian!(J::AbstractMatrix{<:Number}, f, x::AbstractArray{<:Number},
                   fx::AbstractArray{<:Number}, alg::DiffEqBase.AbstractSensitivityAlgorithm, jac_config)
  if alg_autodiff(alg)
    ForwardDiff.jacobian!(J, f, fx, x, jac_config)
  else
    DiffEqDiffTools.finite_difference_jacobian!(J, f, x, jac_config)
  end
  nothing
end

function gradient!(df::AbstractArray{<:Number}, f,
                   x::Union{Number,AbstractArray{<:Number}},
                   alg::DiffEqBase.AbstractSensitivityAlgorithm, grad_config)
    if alg_autodiff(alg)
        ForwardDiff.gradient!(df, f, x, grad_config)
    else
        DiffEqDiffTools.finite_difference_gradient!(df, f, x, grad_config)
    end
    nothing
end

"""
  jacobianvec!(Jv, f, x, v, alg, (buffer, seed)) -> nothing

``Jv <- J(f(x))v``
"""
function jacobianvec!(Jv::AbstractArray{<:Number}, f, x::AbstractArray{<:Number},
                      v, alg::DiffEqBase.AbstractSensitivityAlgorithm, config)
  if alg_autodiff(alg)
    buffer, seed = config
    TD = typeof(first(seed))
    T  = typeof(first(seed).partials)
    @. seed = TD(x, T(tuple(v)))
    f(buffer, seed)
    Jv .= ForwardDiff.partials.(buffer, 1)
  else
    buffer1, buffer2 = config
    f(buffer1,x)
    T = eltype(x)
    # Should it be min? max? mean?
    ϵ = sqrt(eps(real(T))) * max(one(real(T)), abs(norm(x)))
    @. x += ϵ*v
    f(buffer2,x)
    @. x -= ϵ*v
    @. du = (buffer2 - buffer1)/ϵ
  end
  nothing
end

function build_jac_config(alg,uf,u)
  if alg_autodiff(alg)
    jac_config = ForwardDiff.JacobianConfig(uf,u,u,
                 ForwardDiff.Chunk{determine_chunksize(u,alg)}())
  else
    if diff_type(alg) != Val{:complex}
      jac_config = DiffEqDiffTools.JacobianCache(similar(u),similar(u),
                                                 similar(u),diff_type(alg))
    else
      tmp = Complex{eltype(u)}.(u)
      du1 = Complex{eltype(u)}.(du1)
      jac_config = DiffEqDiffTools.JacobianCache(tmp,du1,nothing,diff_type(alg))
    end
  end
  jac_config
end

function build_param_jac_config(alg,uf,u,p)
  if alg_autodiff(alg)
    jac_config = ForwardDiff.JacobianConfig(uf,u,p,
                 ForwardDiff.Chunk{determine_chunksize(p,alg)}())
  else
    if diff_type(alg) != Val{:complex}
      jac_config = DiffEqDiffTools.JacobianCache(similar(p),similar(u),
                                                 similar(u),diff_type(alg))
    else
      tmp = Complex{eltype(p)}.(p)
      du1 = Complex{eltype(u)}.(u)
      jac_config = DiffEqDiffTools.JacobianCache(tmp,du1,nothing,diff_type(alg))
    end
  end
  jac_config
end

function build_grad_config(alg,tf,du1,t)
  if alg_autodiff(alg)
    grad_config = ForwardDiff.GradientConfig(tf,du1,
                                    ForwardDiff.Chunk{determine_chunksize(du1,alg)}())
  else
    grad_config = DiffEqDiffTools.GradientCache(du1,t,diff_type(alg))
  end
  grad_config
end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule NCCL_jll
using Base
using Base: UUID
using LazyArtifacts
import JLLWrappers

JLLWrappers.@generate_main_file_header("NCCL")
JLLWrappers.@generate_main_file("NCCL", UUID("4d6d38e4-5b87-5e63-912a-873ff2d649b7"))
end  # module NCCL_jll


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

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

                     @variable(m, -31.730360167469108 <= q <= 25.84895899492546)

                     add_NL_constraint(m, :((0.0048766473277672695 + 0.6370542772950833*(0.25009028296342795 + -0.2547820908673408*$(x[1]) + -0.7856561704813507*$(x[2]) + -0.058708649178761174*$(x[3]) + -0.7608390940355938*$(x[4]) + -0.05318346349824532*$(x[5]))/(1 + exp(-(0.25009028296342795 + -0.2547820908673408*$(x[1]) + -0.7856561704813507*$(x[2]) + -0.058708649178761174*$(x[3]) + -0.7608390940355938*$(x[4]) + -0.05318346349824532*$(x[5])))) + 0.7379502189935647*(-0.9382943152388461 + 0.5207156183879813*$(x[1]) + -0.2471020507819235*$(x[2]) + -0.7557614646986384*$(x[3]) + -0.33572241425474436*$(x[4]) + -0.6239060401766734*$(x[5]))/(1 + exp(-(-0.9382943152388461 + 0.5207156183879813*$(x[1]) + -0.2471020507819235*$(x[2]) + -0.7557614646986384*$(x[3]) + -0.33572241425474436*$(x[4]) + -0.6239060401766734*$(x[5])))) + -0.017781498731094914*(-0.885143874020959 + 0.9842140185190744*$(x[1]) + 0.788813814030803*$(x[2]) + -0.15763218219342123*$(x[3]) + -0.7069686241723065*$(x[4]) + -0.5626424821908671*$(x[5]))/(1 + exp(-(-0.885143874020959 + 0.9842140185190744*$(x[1]) + 0.788813814030803*$(x[2]) + -0.15763218219342123*$(x[3]) + -0.7069686241723065*$(x[4]) + -0.5626424821908671*$(x[5])))) + -0.30964053180650675*(0.6131135785432997 + 0.10922488949962084*$(x[1]) + -0.12917190851327387*$(x[2]) + -0.1680916348200565*$(x[3]) + -0.8331067755239463*$(x[4]) + 0.6030847772122314*$(x[5]))/(1 + exp(-(0.6131135785432997 + 0.10922488949962084*$(x[1]) + -0.12917190851327387*$(x[2]) + -0.1680916348200565*$(x[3]) + -0.8331067755239463*$(x[4]) + 0.6030847772122314*$(x[5])))) + -0.19693223852539354*(-0.16753936699968897 + -0.5793989304507368*$(x[1]) + 0.5784575652786188*$(x[2]) + 0.7312888898953878*$(x[3]) + 0.9375360522038214*$(x[4]) + -0.8845975912485904*$(x[5]))/(1 + exp(-(-0.16753936699968897 + -0.5793989304507368*$(x[1]) + 0.5784575652786188*$(x[2]) + 0.7312888898953878*$(x[3]) + 0.9375360522038214*$(x[4]) + -0.8845975912485904*$(x[5])))))/(1 + exp(-(0.0048766473277672695 + 0.6370542772950833*(0.25009028296342795 + -0.2547820908673408*$(x[1]) + -0.7856561704813507*$(x[2]) + -0.058708649178761174*$(x[3]) + -0.7608390940355938*$(x[4]) + -0.05318346349824532*$(x[5]))/(1 + exp(-(0.25009028296342795 + -0.2547820908673408*$(x[1]) + -0.7856561704813507*$(x[2]) + -0.058708649178761174*$(x[3]) + -0.7608390940355938*$(x[4]) + -0.05318346349824532*$(x[5])))) + 0.7379502189935647*(-0.9382943152388461 + 0.5207156183879813*$(x[1]) + -0.2471020507819235*$(x[2]) + -0.7557614646986384*$(x[3]) + -0.33572241425474436*$(x[4]) + -0.6239060401766734*$(x[5]))/(1 + exp(-(-0.9382943152388461 + 0.5207156183879813*$(x[1]) + -0.2471020507819235*$(x[2]) + -0.7557614646986384*$(x[3]) + -0.33572241425474436*$(x[4]) + -0.6239060401766734*$(x[5])))) + -0.017781498731094914*(-0.885143874020959 + 0.9842140185190744*$(x[1]) + 0.788813814030803*$(x[2]) + -0.15763218219342123*$(x[3]) + -0.7069686241723065*$(x[4]) + -0.5626424821908671*$(x[5]))/(1 + exp(-(-0.885143874020959 + 0.9842140185190744*$(x[1]) + 0.788813814030803*$(x[2]) + -0.15763218219342123*$(x[3]) + -0.7069686241723065*$(x[4]) + -0.5626424821908671*$(x[5])))) + -0.30964053180650675*(0.6131135785432997 + 0.10922488949962084*$(x[1]) + -0.12917190851327387*$(x[2]) + -0.1680916348200565*$(x[3]) + -0.8331067755239463*$(x[4]) + 0.6030847772122314*$(x[5]))/(1 + exp(-(0.6131135785432997 + 0.10922488949962084*$(x[1]) + -0.12917190851327387*$(x[2]) + -0.1680916348200565*$(x[3]) + -0.8331067755239463*$(x[4]) + 0.6030847772122314*$(x[5])))) + -0.19693223852539354*(-0.16753936699968897 + -0.5793989304507368*$(x[1]) + 0.5784575652786188*$(x[2]) + 0.7312888898953878*$(x[3]) + 0.9375360522038214*$(x[4]) + -0.8845975912485904*$(x[5]))/(1 + exp(-(-0.16753936699968897 + -0.5793989304507368*$(x[1]) + 0.5784575652786188*$(x[2]) + 0.7312888898953878*$(x[3]) + 0.9375360522038214*$(x[4]) + -0.8845975912485904*$(x[5]))))))) + (0.6750746912592471 + -0.6689814871888631*(0.25009028296342795 + -0.2547820908673408*$(x[1]) + -0.7856561704813507*$(x[2]) + -0.058708649178761174*$(x[3]) + -0.7608390940355938*$(x[4]) + -0.05318346349824532*$(x[5]))/(1 + exp(-(0.25009028296342795 + -0.2547820908673408*$(x[1]) + -0.7856561704813507*$(x[2]) + -0.058708649178761174*$(x[3]) + -0.7608390940355938*$(x[4]) + -0.05318346349824532*$(x[5])))) + 0.08174772592079638*(-0.9382943152388461 + 0.5207156183879813*$(x[1]) + -0.2471020507819235*$(x[2]) + -0.7557614646986384*$(x[3]) + -0.33572241425474436*$(x[4]) + -0.6239060401766734*$(x[5]))/(1 + exp(-(-0.9382943152388461 + 0.5207156183879813*$(x[1]) + -0.2471020507819235*$(x[2]) + -0.7557614646986384*$(x[3]) + -0.33572241425474436*$(x[4]) + -0.6239060401766734*$(x[5])))) + 0.43040165086652493*(-0.885143874020959 + 0.9842140185190744*$(x[1]) + 0.788813814030803*$(x[2]) + -0.15763218219342123*$(x[3]) + -0.7069686241723065*$(x[4]) + -0.5626424821908671*$(x[5]))/(1 + exp(-(-0.885143874020959 + 0.9842140185190744*$(x[1]) + 0.788813814030803*$(x[2]) + -0.15763218219342123*$(x[3]) + -0.7069686241723065*$(x[4]) + -0.5626424821908671*$(x[5])))) + 0.6360324545158482*(0.6131135785432997 + 0.10922488949962084*$(x[1]) + -0.12917190851327387*$(x[2]) + -0.1680916348200565*$(x[3]) + -0.8331067755239463*$(x[4]) + 0.6030847772122314*$(x[5]))/(1 + exp(-(0.6131135785432997 + 0.10922488949962084*$(x[1]) + -0.12917190851327387*$(x[2]) + -0.1680916348200565*$(x[3]) + -0.8331067755239463*$(x[4]) + 0.6030847772122314*$(x[5])))) + 0.2590621611090005*(-0.16753936699968897 + -0.5793989304507368*$(x[1]) + 0.5784575652786188*$(x[2]) + 0.7312888898953878*$(x[3]) + 0.9375360522038214*$(x[4]) + -0.8845975912485904*$(x[5]))/(1 + exp(-(-0.16753936699968897 + -0.5793989304507368*$(x[1]) + 0.5784575652786188*$(x[2]) + 0.7312888898953878*$(x[3]) + 0.9375360522038214*$(x[4]) + -0.8845975912485904*$(x[5])))))/(1 + exp(-(0.6750746912592471 + -0.6689814871888631*(0.25009028296342795 + -0.2547820908673408*$(x[1]) + -0.7856561704813507*$(x[2]) + -0.058708649178761174*$(x[3]) + -0.7608390940355938*$(x[4]) + -0.05318346349824532*$(x[5]))/(1 + exp(-(0.25009028296342795 + -0.2547820908673408*$(x[1]) + -0.7856561704813507*$(x[2]) + -0.058708649178761174*$(x[3]) + -0.7608390940355938*$(x[4]) + -0.05318346349824532*$(x[5])))) + 0.08174772592079638*(-0.9382943152388461 + 0.5207156183879813*$(x[1]) + -0.2471020507819235*$(x[2]) + -0.7557614646986384*$(x[3]) + -0.33572241425474436*$(x[4]) + -0.6239060401766734*$(x[5]))/(1 + exp(-(-0.9382943152388461 + 0.5207156183879813*$(x[1]) + -0.2471020507819235*$(x[2]) + -0.7557614646986384*$(x[3]) + -0.33572241425474436*$(x[4]) + -0.6239060401766734*$(x[5])))) + 0.43040165086652493*(-0.885143874020959 + 0.9842140185190744*$(x[1]) + 0.788813814030803*$(x[2]) + -0.15763218219342123*$(x[3]) + -0.7069686241723065*$(x[4]) + -0.5626424821908671*$(x[5]))/(1 + exp(-(-0.885143874020959 + 0.9842140185190744*$(x[1]) + 0.788813814030803*$(x[2]) + -0.15763218219342123*$(x[3]) + -0.7069686241723065*$(x[4]) + -0.5626424821908671*$(x[5])))) + 0.6360324545158482*(0.6131135785432997 + 0.10922488949962084*$(x[1]) + -0.12917190851327387*$(x[2]) + -0.1680916348200565*$(x[3]) + -0.8331067755239463*$(x[4]) + 0.6030847772122314*$(x[5]))/(1 + exp(-(0.6131135785432997 + 0.10922488949962084*$(x[1]) + -0.12917190851327387*$(x[2]) + -0.1680916348200565*$(x[3]) + -0.8331067755239463*$(x[4]) + 0.6030847772122314*$(x[5])))) + 0.2590621611090005*(-0.16753936699968897 + -0.5793989304507368*$(x[1]) + 0.5784575652786188*$(x[2]) + 0.7312888898953878*$(x[3]) + 0.9375360522038214*$(x[4]) + -0.8845975912485904*$(x[5]))/(1 + exp(-(-0.16753936699968897 + -0.5793989304507368*$(x[1]) + 0.5784575652786188*$(x[2]) + 0.7312888898953878*$(x[3]) + 0.9375360522038214*$(x[4]) + -0.8845975912485904*$(x[5]))))))) + (-0.5603004507305789 + -0.6080760485894823*(0.25009028296342795 + -0.2547820908673408*$(x[1]) + -0.7856561704813507*$(x[2]) + -0.058708649178761174*$(x[3]) + -0.7608390940355938*$(x[4]) + -0.05318346349824532*$(x[5]))/(1 + exp(-(0.25009028296342795 + -0.2547820908673408*$(x[1]) + -0.7856561704813507*$(x[2]) + -0.058708649178761174*$(x[3]) + -0.7608390940355938*$(x[4]) + -0.05318346349824532*$(x[5])))) + -0.17827205367530619*(-0.9382943152388461 + 0.5207156183879813*$(x[1]) + -0.2471020507819235*$(x[2]) + -0.7557614646986384*$(x[3]) + -0.33572241425474436*$(x[4]) + -0.6239060401766734*$(x[5]))/(1 + exp(-(-0.9382943152388461 + 0.5207156183879813*$(x[1]) + -0.2471020507819235*$(x[2]) + -0.7557614646986384*$(x[3]) + -0.33572241425474436*$(x[4]) + -0.6239060401766734*$(x[5])))) + -0.3024929990799108*(-0.885143874020959 + 0.9842140185190744*$(x[1]) + 0.788813814030803*$(x[2]) + -0.15763218219342123*$(x[3]) + -0.7069686241723065*$(x[4]) + -0.5626424821908671*$(x[5]))/(1 + exp(-(-0.885143874020959 + 0.9842140185190744*$(x[1]) + 0.788813814030803*$(x[2]) + -0.15763218219342123*$(x[3]) + -0.7069686241723065*$(x[4]) + -0.5626424821908671*$(x[5])))) + 0.2876553916986979*(0.6131135785432997 + 0.10922488949962084*$(x[1]) + -0.12917190851327387*$(x[2]) + -0.1680916348200565*$(x[3]) + -0.8331067755239463*$(x[4]) + 0.6030847772122314*$(x[5]))/(1 + exp(-(0.6131135785432997 + 0.10922488949962084*$(x[1]) + -0.12917190851327387*$(x[2]) + -0.1680916348200565*$(x[3]) + -0.8331067755239463*$(x[4]) + 0.6030847772122314*$(x[5])))) + 0.1394466763243627*(-0.16753936699968897 + -0.5793989304507368*$(x[1]) + 0.5784575652786188*$(x[2]) + 0.7312888898953878*$(x[3]) + 0.9375360522038214*$(x[4]) + -0.8845975912485904*$(x[5]))/(1 + exp(-(-0.16753936699968897 + -0.5793989304507368*$(x[1]) + 0.5784575652786188*$(x[2]) + 0.7312888898953878*$(x[3]) + 0.9375360522038214*$(x[4]) + -0.8845975912485904*$(x[5])))))/(1 + exp(-(-0.5603004507305789 + -0.6080760485894823*(0.25009028296342795 + -0.2547820908673408*$(x[1]) + -0.7856561704813507*$(x[2]) + -0.058708649178761174*$(x[3]) + -0.7608390940355938*$(x[4]) + -0.05318346349824532*$(x[5]))/(1 + exp(-(0.25009028296342795 + -0.2547820908673408*$(x[1]) + -0.7856561704813507*$(x[2]) + -0.058708649178761174*$(x[3]) + -0.7608390940355938*$(x[4]) + -0.05318346349824532*$(x[5])))) + -0.17827205367530619*(-0.9382943152388461 + 0.5207156183879813*$(x[1]) + -0.2471020507819235*$(x[2]) + -0.7557614646986384*$(x[3]) + -0.33572241425474436*$(x[4]) + -0.6239060401766734*$(x[5]))/(1 + exp(-(-0.9382943152388461 + 0.5207156183879813*$(x[1]) + -0.2471020507819235*$(x[2]) + -0.7557614646986384*$(x[3]) + -0.33572241425474436*$(x[4]) + -0.6239060401766734*$(x[5])))) + -0.3024929990799108*(-0.885143874020959 + 0.9842140185190744*$(x[1]) + 0.788813814030803*$(x[2]) + -0.15763218219342123*$(x[3]) + -0.7069686241723065*$(x[4]) + -0.5626424821908671*$(x[5]))/(1 + exp(-(-0.885143874020959 + 0.9842140185190744*$(x[1]) + 0.788813814030803*$(x[2]) + -0.15763218219342123*$(x[3]) + -0.7069686241723065*$(x[4]) + -0.5626424821908671*$(x[5])))) + 0.2876553916986979*(0.6131135785432997 + 0.10922488949962084*$(x[1]) + -0.12917190851327387*$(x[2]) + -0.1680916348200565*$(x[3]) + -0.8331067755239463*$(x[4]) + 0.6030847772122314*$(x[5]))/(1 + exp(-(0.6131135785432997 + 0.10922488949962084*$(x[1]) + -0.12917190851327387*$(x[2]) + -0.1680916348200565*$(x[3]) + -0.8331067755239463*$(x[4]) + 0.6030847772122314*$(x[5])))) + 0.1394466763243627*(-0.16753936699968897 + -0.5793989304507368*$(x[1]) + 0.5784575652786188*$(x[2]) + 0.7312888898953878*$(x[3]) + 0.9375360522038214*$(x[4]) + -0.8845975912485904*$(x[5]))/(1 + exp(-(-0.16753936699968897 + -0.5793989304507368*$(x[1]) + 0.5784575652786188*$(x[2]) + 0.7312888898953878*$(x[3]) + 0.9375360522038214*$(x[4]) + -0.8845975912485904*$(x[5]))))))) + (-0.3104079613918116 + -0.8967757209559282*(0.25009028296342795 + -0.2547820908673408*$(x[1]) + -0.7856561704813507*$(x[2]) + -0.058708649178761174*$(x[3]) + -0.7608390940355938*$(x[4]) + -0.05318346349824532*$(x[5]))/(1 + exp(-(0.25009028296342795 + -0.2547820908673408*$(x[1]) + -0.7856561704813507*$(x[2]) + -0.058708649178761174*$(x[3]) + -0.7608390940355938*$(x[4]) + -0.05318346349824532*$(x[5])))) + 0.9554659562201504*(-0.9382943152388461 + 0.5207156183879813*$(x[1]) + -0.2471020507819235*$(x[2]) + -0.7557614646986384*$(x[3]) + -0.33572241425474436*$(x[4]) + -0.6239060401766734*$(x[5]))/(1 + exp(-(-0.9382943152388461 + 0.5207156183879813*$(x[1]) + -0.2471020507819235*$(x[2]) + -0.7557614646986384*$(x[3]) + -0.33572241425474436*$(x[4]) + -0.6239060401766734*$(x[5])))) + 0.4284452595185697*(-0.885143874020959 + 0.9842140185190744*$(x[1]) + 0.788813814030803*$(x[2]) + -0.15763218219342123*$(x[3]) + -0.7069686241723065*$(x[4]) + -0.5626424821908671*$(x[5]))/(1 + exp(-(-0.885143874020959 + 0.9842140185190744*$(x[1]) + 0.788813814030803*$(x[2]) + -0.15763218219342123*$(x[3]) + -0.7069686241723065*$(x[4]) + -0.5626424821908671*$(x[5])))) + -0.8546417933790886*(0.6131135785432997 + 0.10922488949962084*$(x[1]) + -0.12917190851327387*$(x[2]) + -0.1680916348200565*$(x[3]) + -0.8331067755239463*$(x[4]) + 0.6030847772122314*$(x[5]))/(1 + exp(-(0.6131135785432997 + 0.10922488949962084*$(x[1]) + -0.12917190851327387*$(x[2]) + -0.1680916348200565*$(x[3]) + -0.8331067755239463*$(x[4]) + 0.6030847772122314*$(x[5])))) + -0.9550224035691546*(-0.16753936699968897 + -0.5793989304507368*$(x[1]) + 0.5784575652786188*$(x[2]) + 0.7312888898953878*$(x[3]) + 0.9375360522038214*$(x[4]) + -0.8845975912485904*$(x[5]))/(1 + exp(-(-0.16753936699968897 + -0.5793989304507368*$(x[1]) + 0.5784575652786188*$(x[2]) + 0.7312888898953878*$(x[3]) + 0.9375360522038214*$(x[4]) + -0.8845975912485904*$(x[5])))))/(1 + exp(-(-0.3104079613918116 + -0.8967757209559282*(0.25009028296342795 + -0.2547820908673408*$(x[1]) + -0.7856561704813507*$(x[2]) + -0.058708649178761174*$(x[3]) + -0.7608390940355938*$(x[4]) + -0.05318346349824532*$(x[5]))/(1 + exp(-(0.25009028296342795 + -0.2547820908673408*$(x[1]) + -0.7856561704813507*$(x[2]) + -0.058708649178761174*$(x[3]) + -0.7608390940355938*$(x[4]) + -0.05318346349824532*$(x[5])))) + 0.9554659562201504*(-0.9382943152388461 + 0.5207156183879813*$(x[1]) + -0.2471020507819235*$(x[2]) + -0.7557614646986384*$(x[3]) + -0.33572241425474436*$(x[4]) + -0.6239060401766734*$(x[5]))/(1 + exp(-(-0.9382943152388461 + 0.5207156183879813*$(x[1]) + -0.2471020507819235*$(x[2]) + -0.7557614646986384*$(x[3]) + -0.33572241425474436*$(x[4]) + -0.6239060401766734*$(x[5])))) + 0.4284452595185697*(-0.885143874020959 + 0.9842140185190744*$(x[1]) + 0.788813814030803*$(x[2]) + -0.15763218219342123*$(x[3]) + -0.7069686241723065*$(x[4]) + -0.5626424821908671*$(x[5]))/(1 + exp(-(-0.885143874020959 + 0.9842140185190744*$(x[1]) + 0.788813814030803*$(x[2]) + -0.15763218219342123*$(x[3]) + -0.7069686241723065*$(x[4]) + -0.5626424821908671*$(x[5])))) + -0.8546417933790886*(0.6131135785432997 + 0.10922488949962084*$(x[1]) + -0.12917190851327387*$(x[2]) + -0.1680916348200565*$(x[3]) + -0.8331067755239463*$(x[4]) + 0.6030847772122314*$(x[5]))/(1 + exp(-(0.6131135785432997 + 0.10922488949962084*$(x[1]) + -0.12917190851327387*$(x[2]) + -0.1680916348200565*$(x[3]) + -0.8331067755239463*$(x[4]) + 0.6030847772122314*$(x[5])))) + -0.9550224035691546*(-0.16753936699968897 + -0.5793989304507368*$(x[1]) + 0.5784575652786188*$(x[2]) + 0.7312888898953878*$(x[3]) + 0.9375360522038214*$(x[4]) + -0.8845975912485904*$(x[5]))/(1 + exp(-(-0.16753936699968897 + -0.5793989304507368*$(x[1]) + 0.5784575652786188*$(x[2]) + 0.7312888898953878*$(x[3]) + 0.9375360522038214*$(x[4]) + -0.8845975912485904*$(x[5]))))))) + (0.11524778075403264 + -0.08842997013200282*(0.25009028296342795 + -0.2547820908673408*$(x[1]) + -0.7856561704813507*$(x[2]) + -0.058708649178761174*$(x[3]) + -0.7608390940355938*$(x[4]) + -0.05318346349824532*$(x[5]))/(1 + exp(-(0.25009028296342795 + -0.2547820908673408*$(x[1]) + -0.7856561704813507*$(x[2]) + -0.058708649178761174*$(x[3]) + -0.7608390940355938*$(x[4]) + -0.05318346349824532*$(x[5])))) + -0.2446780712327863*(-0.9382943152388461 + 0.5207156183879813*$(x[1]) + -0.2471020507819235*$(x[2]) + -0.7557614646986384*$(x[3]) + -0.33572241425474436*$(x[4]) + -0.6239060401766734*$(x[5]))/(1 + exp(-(-0.9382943152388461 + 0.5207156183879813*$(x[1]) + -0.2471020507819235*$(x[2]) + -0.7557614646986384*$(x[3]) + -0.33572241425474436*$(x[4]) + -0.6239060401766734*$(x[5])))) + 0.5457935761049009*(-0.885143874020959 + 0.9842140185190744*$(x[1]) + 0.788813814030803*$(x[2]) + -0.15763218219342123*$(x[3]) + -0.7069686241723065*$(x[4]) + -0.5626424821908671*$(x[5]))/(1 + exp(-(-0.885143874020959 + 0.9842140185190744*$(x[1]) + 0.788813814030803*$(x[2]) + -0.15763218219342123*$(x[3]) + -0.7069686241723065*$(x[4]) + -0.5626424821908671*$(x[5])))) + -0.1673283141057702*(0.6131135785432997 + 0.10922488949962084*$(x[1]) + -0.12917190851327387*$(x[2]) + -0.1680916348200565*$(x[3]) + -0.8331067755239463*$(x[4]) + 0.6030847772122314*$(x[5]))/(1 + exp(-(0.6131135785432997 + 0.10922488949962084*$(x[1]) + -0.12917190851327387*$(x[2]) + -0.1680916348200565*$(x[3]) + -0.8331067755239463*$(x[4]) + 0.6030847772122314*$(x[5])))) + 0.6343389483483679*(-0.16753936699968897 + -0.5793989304507368*$(x[1]) + 0.5784575652786188*$(x[2]) + 0.7312888898953878*$(x[3]) + 0.9375360522038214*$(x[4]) + -0.8845975912485904*$(x[5]))/(1 + exp(-(-0.16753936699968897 + -0.5793989304507368*$(x[1]) + 0.5784575652786188*$(x[2]) + 0.7312888898953878*$(x[3]) + 0.9375360522038214*$(x[4]) + -0.8845975912485904*$(x[5])))))/(1 + exp(-(0.11524778075403264 + -0.08842997013200282*(0.25009028296342795 + -0.2547820908673408*$(x[1]) + -0.7856561704813507*$(x[2]) + -0.058708649178761174*$(x[3]) + -0.7608390940355938*$(x[4]) + -0.05318346349824532*$(x[5]))/(1 + exp(-(0.25009028296342795 + -0.2547820908673408*$(x[1]) + -0.7856561704813507*$(x[2]) + -0.058708649178761174*$(x[3]) + -0.7608390940355938*$(x[4]) + -0.05318346349824532*$(x[5])))) + -0.2446780712327863*(-0.9382943152388461 + 0.5207156183879813*$(x[1]) + -0.2471020507819235*$(x[2]) + -0.7557614646986384*$(x[3]) + -0.33572241425474436*$(x[4]) + -0.6239060401766734*$(x[5]))/(1 + exp(-(-0.9382943152388461 + 0.5207156183879813*$(x[1]) + -0.2471020507819235*$(x[2]) + -0.7557614646986384*$(x[3]) + -0.33572241425474436*$(x[4]) + -0.6239060401766734*$(x[5])))) + 0.5457935761049009*(-0.885143874020959 + 0.9842140185190744*$(x[1]) + 0.788813814030803*$(x[2]) + -0.15763218219342123*$(x[3]) + -0.7069686241723065*$(x[4]) + -0.5626424821908671*$(x[5]))/(1 + exp(-(-0.885143874020959 + 0.9842140185190744*$(x[1]) + 0.788813814030803*$(x[2]) + -0.15763218219342123*$(x[3]) + -0.7069686241723065*$(x[4]) + -0.5626424821908671*$(x[5])))) + -0.1673283141057702*(0.6131135785432997 + 0.10922488949962084*$(x[1]) + -0.12917190851327387*$(x[2]) + -0.1680916348200565*$(x[3]) + -0.8331067755239463*$(x[4]) + 0.6030847772122314*$(x[5]))/(1 + exp(-(0.6131135785432997 + 0.10922488949962084*$(x[1]) + -0.12917190851327387*$(x[2]) + -0.1680916348200565*$(x[3]) + -0.8331067755239463*$(x[4]) + 0.6030847772122314*$(x[5])))) + 0.6343389483483679*(-0.16753936699968897 + -0.5793989304507368*$(x[1]) + 0.5784575652786188*$(x[2]) + 0.7312888898953878*$(x[3]) + 0.9375360522038214*$(x[4]) + -0.8845975912485904*$(x[5]))/(1 + exp(-(-0.16753936699968897 + -0.5793989304507368*$(x[1]) + 0.5784575652786188*$(x[2]) + 0.7312888898953878*$(x[3]) + 0.9375360522038214*$(x[4]) + -0.8845975912485904*$(x[5]))))))) - $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:
@testset "Regions" begin
  @testset "Rectangle" begin
    r = RectangleRegion((1.,1.), (2.,3.))
    @test GeoStatsBase.center(r) == [3/2, 4/2]
    @test lowerleft(r) == [1., 1.]
    @test upperright(r) == [2., 3.]
    @test sides(r) == [1., 2.]
    @test GeoStatsBase.diagonal(r) == sqrt(1^2 + 2^2)
    @test volume(r) == 1*2
  end
end


text:
using NPZ
using Distributions
using LinearAlgebra
using Distributed
using StatsBase
using Distributions
using SpecialFunctions
using LinearAlgebra
using Random

function generate_sph_gaussian_data(N::Int64, D::Int64, K::Int64)
	x = randn(D,N)
	tpi = rand(Dirichlet(ones(K)))
	tzn = rand(Multinomial(N,tpi))
	tz = zeros(N)

	tmean = zeros(D,K)
	tcov = zeros(D,D,K)

	ind = 1
	println(tzn)
	for i=1:length(tzn)
		indices = ind:ind+tzn[i]-1
		tz[indices] .= i
		tmean[:,i] .= rand(MvNormal(zeros(D), 100*Matrix{Float64}(I, D, D)))
		tcov[:,:,i] .= rand(InverseGamma((D+2)/2,1))*Matrix{Float64}(I, D, D)
		# T = chol(slice(tcov,:,:,i))
		# x[:,indices] = broadcast(+, T*x[:,indices], tmean[:,i]);
		d = MvNormal(tmean[:,i], tcov[:,:,i])
		for j=indices
			x[:,j] = rand(d)
		end
		ind += tzn[i]
	end
	x, tz, tmean, tcov
end

function generate_gaussian_data(N::Int64, D::Int64, K::Int64)
	x = randn(D,N)
	tpi = rand(Dirichlet(ones(K)))
	tzn = rand(Multinomial(N,tpi))
	tz = zeros(N)

	tmean = zeros(D,K)
	tcov = zeros(D,D,K)

	ind = 1
	println(tzn)
	for i=1:length(tzn)
		indices = ind:ind+tzn[i]-1
		tz[indices] .= i
		tmean[:,i] .= rand(MvNormal(zeros(D), 100*Matrix{Float64}(I, D, D)))
		tcov[:,:,i] .= rand(InverseWishart(D+2, Matrix{Float64}(I, D, D)))
		# T = chol(slice(tcov,:,:,i))
		# x[:,indices] = broadcast(+, T'*x[:,indices], tmean[:,i]);
		d = MvNormal(tmean[:,i], tcov[:,:,i])
		for j=indices
			x[:,j] = rand(d)
		end
		ind += tzn[i]
	end
	x, tz, tmean, tcov
end

function generate_mnmm_data(N::Int64, D::Int64, K::Int64, trials::Int64)
	clusters = zeros(D,K)
	x = zeros(D,N)
	labels = rand(1:K,(N,))
	for i=1:K
		alphas = rand(1:20,(D,))
		alphas[i] = rand(30:100)
		clusters[:,i] = rand(Dirichlet(alphas))
	end
	for i=1:N
		x[:,i] = rand(Multinomial(trials,clusters[:,labels[i]]))
	end
	return x, labels, clusters
end

function generate_grouped_mnm_data(N::Int64, D_global::Int64, D_local, K_global::Int64, K_local::Int64, groups_count::Int64, rand_local::Bool, trials::Int64)
	global_weights = rand(Dirichlet(ones(K_global)))
	pts_dict = Dict()
	labels_dict = Dict()
	clusters_g = zeros(D_global,K_global)
	for i=1:length(global_weights)
		alphas = ones(D_global)*2
		alphas[rand(1:D_global,Int(floor(D_global/40)))] .= rand(trials/2:trials)
		clusters_g[:,i] = rand(Dirichlet(alphas))
	end
	for j=1:groups_count
		x = zeros(D_global + D_local,N)
		x_labels = zeros(size(x,2),2)
		if rand_local
			local_k = rand(1:K_local*2)
		else
			local_k = K_local
		end
		clusters_tzn = sample(1:length(global_weights),ProbabilityWeights(global_weights),local_k)
		local_weights = rand(Dirichlet(ones(local_k)*100))
		clusters_l = zeros(D_global+D_local,local_k)
		ind = 1
		group_ind = 1
		group_tzn = rand(Multinomial(N,local_weights))
		# group_tzn = rand(Multinomial(N,ones(local_k)))
		for i=1:length(local_weights)
			indices = ind:ind+group_tzn[i]-1
			# println(indices)

			alphas = ones(D_local)*2
			alphas[rand(1:D_local,Int(floor(D_local/1.25)))] .= rand(trials/2:trials)
			# local_part = rand(Dirichlet(alphas))
			# clusters_l[:,i] .= cat(clusters_g[:,clusters_tzn[i]], alphas, dims = [1])

			pvector = rand(Dirichlet(alphas))
			for a=indices
				x[1:D_global,a] = rand(Multinomial(trials,clusters_g[:,clusters_tzn[i]]))
				x[D_global+1:end,a] = rand(Multinomial(trials,pvector))
				x_labels[a,1] = clusters_tzn[i]
				x_labels[a,2] = i
			end
			ind += group_tzn[i]
		end
		pts_dict[j] = x
		labels_dict[j] = x_labels
	end
	return pts_dict, labels_dict
end

function generate_grouped_gaussian_data(N::Int64, D_global::Int64, D_local, K_global::Int64, K_local::Int64, groups_count::Int64, rand_local::Bool, var_size, hdp=false)
	global_weights = rand(Dirichlet(ones(K_global)))
	pts_dict = Dict()
	labels_dict = Dict()
	tmean = zeros(D_global,K_global)
	tcov = zeros(D_global,D_global,K_global)
	for i=1:length(global_weights)
		tmean[:,i] .= rand(MvNormal(zeros(D_global), var_size*Matrix{Float64}(I, D_global, D_global)))
		tcov[:,:,i] .= rand(InverseWishart(D_global+2, Matrix{Float64}(I, D_global, D_global)))
	end
	for j=1:groups_count
		x = randn(D_global + D_local,N)
		x_labels = zeros(size(x,2),2)
		if rand_local
			local_k = rand(K_local-2:K_local+2)
		else
			local_k = K_local
		end
		clusters_tzn = sample(1:length(global_weights),ProbabilityWeights(global_weights),local_k)
		local_weights = rand(Dirichlet(ones(local_k)*100))
		group_mean = zeros(D_global+D_local,local_k)
		group_cov = zeros(D_global+D_local,D_global+D_local,local_k)
		ind = 1
		group_ind = 1
		group_tzn = rand(Multinomial(N,local_weights))
		# group_tzn = rand(Multinomial(N,ones(local_k)))
		for i=1:length(local_weights)
			indices = ind:ind+group_tzn[i]-1
			g_mean = rand(MvNormal(zeros(D_local), var_size*Matrix{Float64}(I, D_local, D_local)))
			g_cov = rand(InverseWishart(D_local+2, Matrix{Float64}(I, D_local, D_local)))
			group_mean[:,i] .= cat(tmean[:,clusters_tzn[i]], g_mean, dims = [1])
			group_cov[:,:,i] .= cat(tcov[:,:,clusters_tzn[i]], g_cov, dims = [1,2])
			d = MvNormal(group_mean[:,i], group_cov[:,:,i])
			for j=indices
				x[:,j] = rand(d)
				x_labels[j,1] = clusters_tzn[i]
				x_labels[j,2] = i
			end
			ind += group_tzn[i]
		end
		if hdp
			x[D_global+1:end,:] .= 0
		end
		pts_dict[j] = x
		labels_dict[j] = x_labels
	end
	return pts_dict, labels_dict
end

function create_mnmm_data()
	 x,labels,clusters = generate_mnmm_data(10^6,100,6,1000)
	 npzwrite("data/mnmm/1milD100K6.npy",x')
	 x,labels,clusters = generate_mnmm_data(10^6,100,60,100)
	 npzwrite("data/mnmm/1milD100K60.npy",x')
 end

function create_gaussian_data()
	x,labels,clusters = generate_gaussian_data(10^5,2,20)
	npzwrite("data/2d-1mil/samples_100k1.npy",x')
	# x,labels,clusters = generate_gaussian_data(10^5,100,60,100)
	# npzwrite("data/30d-1mil/samples.npy",x')
end

function create_grouped_gaussian_data()
	samples_count = 10000
	global_dim = 2
	local_dim = 1
	global_clusters_count = 4
	local_clusters_count = 5
	groups_count = 5
	random_local_cluster_count = true
	var_size = 60
	hdp = true

	path_to_save = "data/3d-gaussians/"
	prefix = "gg"

	pts, labels = generate_grouped_gaussian_data(samples_count,
		global_dim,
		local_dim,
		global_clusters_count,
		local_clusters_count,
		groups_count,
		random_local_cluster_count,
		var_size,
		hdp)
	for (k,v) in pts
		npzwrite(path_to_save * prefix * string(k) * ".npy",v')
	end
	for (k,v) in labels
		npzwrite(path_to_save * prefix*"_labels" * string(k) * ".npy",v)
	end
end

function create_grouped_mnmm_data()
	samples_count = 10000
	global_dim = 100
	local_dim = 100
	global_clusters_count = 10
	local_clusters_count = 20
	groups_count = 4
	random_local_cluster_count = false
	trials = 1000

	path_to_save = "data/200d-mnm/"
	prefix = "g"

	pts, labels = generate_grouped_mnm_data(samples_count,
		global_dim,
		local_dim,
		global_clusters_count,
		local_clusters_count,
		groups_count,
		random_local_cluster_count,
		trials)
	for (k,v) in pts
		npzwrite(path_to_save * prefix * string(k) * ".npy",v')
	end
	for (k,v) in labels
		npzwrite(path_to_save * prefix*"_labels" * string(k) * ".npy",v)
	end
end

function single_crp_draw(tables,α)
	tables = push!(tables,α)
	sum_tables = sum(tables)
	table_probs = [x / sum_tables for x in tables]
	return sample(ProbabilityWeights(table_probs))
end

function hdp_prior_crf_draws(N,J,α,γ)
	groups_tables_counts = Dict()
	for j=1:J
		table_counts = []
		for i=1:N
			point_table = single_crp_draw(table_counts[:],α)
			if point_table == length(table_counts)+1
				push!(table_counts,1)
			else
				table_counts[point_table] += 1
			end
		end
		groups_tables_counts[j] = table_counts
	end
	global_tables_count = []
	groups_tables_assignments = Dict([i=>[] for i=1:J])
	cur_group_table = Dict([i=>1 for i=1:J])
	is_done = [false for i=1:J]
	while any(is_done .== false)
		for j=1:J
			group_assignments = groups_tables_assignments[j]
			i = cur_group_table[j]
			if i > length(groups_tables_counts[j])
				is_done[j] = true
				continue
			end
			table_assignment = single_crp_draw(global_tables_count,γ)
			if table_assignment == length(global_tables_count)+1
				push!(global_tables_count,1)
			else
				global_tables_count[table_assignment] += 1
			end
			push!(group_assignments,table_assignment)
			groups_tables_assignments[j] = group_assignments
			cur_group_table[j]+=1
		end
	end
	global_mixture = [x / sum(global_tables_count) for x in global_tables_count]
	groups_mixtures = Dict()
	for j=1:J
		local_mixture = zeros(length(global_mixture))
		for i = 1:length(global_mixture)
			if i in groups_tables_assignments[j]
				local_mixture[i] = sum(groups_tables_counts[j][groups_tables_assignments[j] .== i])
			end
		end
		groups_mixtures[j] = local_mixture
	end
	return global_mixture,groups_mixtures,groups_tables_counts,groups_tables_assignments
end

function generate_grouped_gaussian_from_hdp_group_counts(group_counts, dim, var_size)
	pts_dict = Dict()
	labels_dict = Dict()
	K = length(group_counts[1])
	tmean = zeros(dim,K)
	tcov = zeros(dim,dim,K)
	J = length(group_counts)
	components = []
	for i=1:length(group_counts[1])
		tmean[:,i] .= rand(MvNormal(zeros(dim), var_size*Matrix{Float64}(I, dim, dim)))
		tcov[:,:,i] .= rand(InverseWishart(dim+2, Matrix{Float64}(I, dim, dim)))
		push!(components,MvNormal(tmean[:,i], tcov[:,:,i]))
	end
	for j=1:J
		points = reduce(hcat,[rand(components[i],Int(group_counts[j][i])) for i=1:length(group_counts[j])])
		labels = reduce(vcat,[Int.(ones(Int(group_counts[j][i])))*i for i=1:length(group_counts[j])])
		pts_dict[j] = points
		labels_dict[j] = labels
	end
	return pts_dict, labels_dict
end