text:
isfrozenvalue(x) = isfrozenvalue(typeof(x))
isfrozenindex(x) = isfrozenindex(typeof(x))
isfrozen(x) = isfrozen(typeof(x))

isfrozenvalue(::Type{T}) where {T} = isbitstype(T)
isfrozenindex(::Type{T}) where {T} = isbitstype(T)
isfrozen(::Type{T}) where {T} = isfrozenvalue(T) && isfrozenindex(T)

isimmutablevalue(x) = isimmutablevalue(typeof(x))
isimmutableindex(x) = isimmutableindex(typeof(x))

isimmutablevalue(::Type{T}) where {T} = isfrozenvalue(T)
isimmutableindex(::Type{T}) where {T} = isfrozenindex(T)

"""
    readonly(x) -> z

Create a read-only view of `x`.  Mutations on `x` are reflected to `z`.
"""
readonly

"""
    melt(z) -> x
    meltvalue(z) -> x
    meltindex(z) -> x
    freeze(x) -> z
    freezevalue(x) -> z
    freezeindex(x) -> z

`melt` and `meltvalue` create a mutable *copy* of `x`.  `freeze`,
`freezevalue`, and `freezeindex` create an immutable *copy* of `x`.

The result of `melt(::AbstractVector)` is appendable. The result of
`meltvalue(::AbstractVector)` may not be appendable.  `meltindex`
undoes `freezeindex`.

`freezevalue` only freezes the existing values and it is still
possible to append the new items to `z`.  `freezeindex` only freezes
the indices and it is still possible to mutate the values.

Use, e.g., `freeze(move!(x))` and `melt(move!(z))` to freeze or melt the
values without creating a copy (see also [`move!`](@ref)).

[`readonly`](@ref) can also be used to create a read-only view without
creating a copy and without asserting strict absence of ownership.

# Examples
```jldoctest
julia> using Mutabilities

julia> z = freeze([1, 2, 3])
3-element freeze(::Array{Int64,1}) with eltype Int64:
 1
 2
 3

julia> z[1] = 111
ERROR: setindex! not defined for Mutabilities.ImmutableArray{Int64,1,Array{Int64,1}}

julia> y = melt(z)
3-element Array{Int64,1}:
 1
 2
 3
```
"""
(melt, meltvalue, meltindex, freeze, freezevalue, freezeindex)
# TODO: `melt(::Matrix)` can make a matrix appendable. What to do with it?
# Maybe only allow maximum mutability supported by its abstract type
# interface? (current implementation)

readonly(x) = isimmutablevalue(x) && isimmutableindex(x) ? x : readonly_impl(x)
# maybe just?: readonly(x) = freeze(move!(x))

melt(x) = meltvalue(move!(meltindex(x)))
meltvalue(x) = isimmutablevalue(x) ? ownmutablevalue(x) : own(x)
meltindex(x) = isimmutableindex(x) ? ownmutableindex(x) : own(x)

freeze(x) = isfrozen(x) ? unmoved(x) : freeze!(own(x))
freezevalue(x) = isfrozenvalue(x) ? unmoved(x) : freezevalue!(own(x))
freezeindex(x) = isfrozenindex(x) ? unmoved(x) : freezeindex!(own(x))

function freeze! end
function freezevalue! end
function freezeindex! end

meltvalue!(x) = x
meltindex!(x) = x

# `f!(x)` is the version of `f(x)` where the caller asserts that `x`
# is not going to be used by anyone else. i.e., the caller "allows"
# any future mutations of `x`.

own(x) = copy(x)
# own(x) = deepcopy(x)

ownmutablevalue(x) = Base.copymutable(x)
ownmutableindex(x) = Base.copymutable(x)

struct Moved{T}
    value::T
end

own(x::Moved) = x.value
ownmutablevalue(x::Moved) = meltvalue!(x.value)
ownmutableindex(x::Moved) = meltindex!(x.value)

unmoved(x::Moved) = x.value
unmoved(x) = x

isimmutablevalue(::Type{Moved{T}}) where {T} = isimmutablevalue(T)
isimmutableindex(::Type{Moved{T}}) where {T} = isimmutableindex(T)
isfrozenvalue(::Type{Moved{T}}) where {T} = isfrozenvalue(T)
isfrozenindex(::Type{Moved{T}}) where {T} = isfrozenindex(T)

readonly(x::Moved) = readonly(x.value)  # ignore "move!"

"""
    move!(x)

Manually declare that the object `x` has no other owners and the
object `x` is not going to be used by the caller.

# Examples
```jldoctest
julia> using Mutabilities

julia> x = [];

julia> melt(move!(freeze(move!(x)))) === x
true

julia> melt(move!(freeze(x))) === x
false

julia> melt(freeze(move!(x))) === x
false
```

!!! note

    Above examples intentionally violate the rule for using `move!` to
    show how it works.  When `x` is passed to `move!` on the left hand
    side, it is not allowed to use `x` on the right hand side of `===`.
"""
move!(x) = Moved(x)
move!(x::Moved) = x


text:
@testset "pairwise" begin
    rng = MersenneTwister(123456)
    d = SqEuclidean()
    Ns = (4, 5)
    D = 3
    x = [randn(rng, D) for _ in 1:Ns[1]]
    y = [randn(rng, D) for _ in 1:Ns[2]]
    X = hcat(x...)
    Y = hcat(y...)
    K = zeros(Ns)

    @test KernelFunctions.pairwise(d, x, y) ≈ pairwise(d, X, Y, dims=2)
    @test KernelFunctions.pairwise(d, x) ≈ pairwise(d, X, dims=2)
    KernelFunctions.pairwise!(K, d, x, y)
    @test K ≈ pairwise(d, X, Y, dims=2)
    K = zeros(Ns[1], Ns[1])
    KernelFunctions.pairwise!(K, d, x)
    @test K ≈ pairwise(d, X, dims=2)

    x = randn(rng, 10)
    X = reshape(x, :, 1)
    y = randn(rng, 11)
    Y = reshape(y, :, 1)
    K = zeros(10, 11)
    @test KernelFunctions.pairwise(d, x, y) ≈ pairwise(d, X, Y; dims=1)
    @test KernelFunctions.pairwise(d, x) ≈ pairwise(d, X; dims=1)
    KernelFunctions.pairwise!(K, d, x, y)
    @test K ≈ pairwise(d, X, Y, dims=1)
end

ext: jl


text:
####################################################################
# In this example we demonstrate using the Trimming solver.
#
# We test the l2 penalty (with trimming) here. As such, we use
# an example problem where the outlier noise is limited to a
# few specific sensors (this corresponds to mode=2 in the genDMD
# routine that creates the synthetic data)
#
# We also use a prox operator that forces the real part of the
# eigenvalues alpha to be less than or equal
# to one. It is generally advisable to put some bound
# in order to keep e^(Re(alpha)*tmax) from being too large,
# e.g. overflow may occur without such a choice

using RobustDMD

T = Float32

#--------------------------------------------------------------------
# Generate DMD Synthetic Data
#--------------------------------------------------------------------
println("generate synthetic data...")
# dimensions
m = 100;    # temporal dimension
n = 100;    # spatial dimension
h = 80;     # trimming level
k = 3;      # number of modes
# generate data
sigma = T(1e-4); # size of background noise
mu = T(1); # size of spikes
p = T(0.1); # frequency of spikes (as mode=2 this determines the percentage
# of columns that are corrupted)
X, Xclean, t, alphat, betat = genDMD(m, n, k, sigma, mu; seed=123456, p=p, mode=2);

#--------------------------------------------------------------------
# Initial Guess
#--------------------------------------------------------------------
println("initialize alpha and B...")
# create initial guess using trap. rule with exact dmd
alpha0, B0 = dmdexactestimate(m,n,k,X,t);

#--------------------------------------------------------------------
# Trim l2
#--------------------------------------------------------------------
println("Trim l2 loss experiment...")
# loss functions
lossf = (z) -> l2_func(z);
lossg = (z) -> l2_grad!(z);
l2Params = DMDParams(k, X, t, lossf, lossg);
l2DMD = DMDVariables(alpha0, B0, l2Params);
# wrapper: sets up inner solve specialized for l2
l2Svars = DMDVPSolverVariablesLSQ(l2Params);
# prox operator
function prox_mr(αr)
    k = length(αr) >> 1;
    for i = 1:k
        αr[i<<1-1] = min(T(1.0), αr[i<<1-1]);
    end
end

# apply solver
options = Trim_options(h,300,T(1e-5),T(1.2),true,true,10,prox=prox_mr);
Trim_solve_DMD!(l2DMD, l2Params, l2Svars, options);

#--------------------------------------------------------------------
# Compare Result
#--------------------------------------------------------------------
println("compare result...")
err0 = besterrperm(alpha0, alphat);
err1 = besterrperm(l2DMD.alpha, alphat);
println("err exact DMD ..... ", err0);
println("err w/ trimming ... ", err1);


text:
export BrickletSolidStateRelayV2Monoflop
struct BrickletSolidStateRelayV2Monoflop
    state::Bool
    time::Integer
    time_remaining::Integer
end

export BrickletSolidStateRelayV2SPITFPErrorCount
struct BrickletSolidStateRelayV2SPITFPErrorCount
    error_count_ack_checksum::Integer
    error_count_message_checksum::Integer
    error_count_frame::Integer
    error_count_overflow::Integer
end

export BrickletSolidStateRelayV2Identity
struct BrickletSolidStateRelayV2Identity
    uid::String
    connected_uid::String
    position::Char
    hardware_version::Vector{Integer}
    firmware_version::Vector{Integer}
    device_identifier::Integer
end

export BrickletSolidStateRelayV2
"""
Controls AC and DC Solid State Relays
"""
mutable struct BrickletSolidStateRelayV2 <: TinkerforgeDevice
    ipcon::IPConnection
    deviceInternal::PyObject

    """
    Creates an object with the unique device ID *uid* and adds it to
    the IP Connection *ipcon*.
    """
    function BrickletSolidStateRelayV2(uid::String, ipcon::IPConnection)
        package = pyimport("tinkerforge.bricklet_solid_state_relay_v2")
        deviceInternal = package.BrickletSolidStateRelayV2(uid, ipcon.ipconInternal)

        return new(ipcon, deviceInternal)
    end
end

export set_state
"""
    $(SIGNATURES)

Sets the state of the relays *true* means on and *false* means off.

A running monoflop timer will be aborted if this function is called.
"""
function set_state(device::BrickletSolidStateRelayV2, state)
    device.deviceInternal.set_state(state)
end

export get_state
"""
    $(SIGNATURES)

Returns the state of the relay, *true* means on and *false* means off.
"""
function get_state(device::BrickletSolidStateRelayV2)
    return device.deviceInternal.get_state()
end

export set_monoflop
"""
    $(SIGNATURES)

The first parameter  is the desired state of the relay (*true* means on
and *false* means off). The second parameter indicates the time that
the relay should hold the state.

If this function is called with the parameters (true, 1500):
The relay will turn on and in 1.5s it will turn off again.

A monoflop can be used as a failsafe mechanism. For example: Lets assume you
have a RS485 bus and a Solid State Relay Bricklet connected to one of the slave
stacks. You can now call this function every second, with a time parameter
of two seconds. The relay will be on all the time. If now the RS485
connection is lost, the relay will turn off in at most two seconds.
"""
function set_monoflop(device::BrickletSolidStateRelayV2, state, time)
    device.deviceInternal.set_monoflop(state, time)
end

export get_monoflop
"""
    $(SIGNATURES)

Returns the current state and the time as set by
:func:`Set Monoflop` as well as the remaining time until the state flips.

If the timer is not running currently, the remaining time will be returned
as 0.
"""
function get_monoflop(device::BrickletSolidStateRelayV2)
    return device.deviceInternal.get_monoflop()
end

export get_spitfp_error_count
"""
    $(SIGNATURES)

Returns the error count for the communication between Brick and Bricklet.

The errors are divided into

* ACK checksum errors,
* message checksum errors,
* framing errors and
* overflow errors.

The errors counts are for errors that occur on the Bricklet side. All
Bricks have a similar function that returns the errors on the Brick side.
"""
function get_spitfp_error_count(device::BrickletSolidStateRelayV2)
    return device.deviceInternal.get_spitfp_error_count()
end

export set_bootloader_mode
"""
    $(SIGNATURES)

Sets the bootloader mode and returns the status after the requested
mode change was instigated.

You can change from bootloader mode to firmware mode and vice versa. A change
from bootloader mode to firmware mode will only take place if the entry function,
device identifier and CRC are present and correct.

This function is used by Brick Viewer during flashing. It should not be
necessary to call it in a normal user program.
"""
function set_bootloader_mode(device::BrickletSolidStateRelayV2, mode)
    return device.deviceInternal.set_bootloader_mode(mode)
end

export get_bootloader_mode
"""
    $(SIGNATURES)

Returns the current bootloader mode, see :func:`Set Bootloader Mode`.
"""
function get_bootloader_mode(device::BrickletSolidStateRelayV2)
    return device.deviceInternal.get_bootloader_mode()
end

export set_write_firmware_pointer
"""
    $(SIGNATURES)

Sets the firmware pointer for :func:`Write Firmware`. The pointer has
to be increased by chunks of size 64. The data is written to flash
every 4 chunks (which equals to one page of size 256).

This function is used by Brick Viewer during flashing. It should not be
necessary to call it in a normal user program.
"""
function set_write_firmware_pointer(device::BrickletSolidStateRelayV2, pointer)
    device.deviceInternal.set_write_firmware_pointer(pointer)
end

export write_firmware
"""
    $(SIGNATURES)

Writes 64 Bytes of firmware at the position as written by
:func:`Set Write Firmware Pointer` before. The firmware is written
to flash every 4 chunks.

You can only write firmware in bootloader mode.

This function is used by Brick Viewer during flashing. It should not be
necessary to call it in a normal user program.
"""
function write_firmware(device::BrickletSolidStateRelayV2, data)
    return device.deviceInternal.write_firmware(data)
end

export set_status_led_config
"""
    $(SIGNATURES)

Sets the status LED configuration. By default the LED shows
communication traffic between Brick and Bricklet, it flickers once
for every 10 received data packets.

You can also turn the LED permanently on/off or show a heartbeat.

If the Bricklet is in bootloader mode, the LED is will show heartbeat by default.
"""
function set_status_led_config(device::BrickletSolidStateRelayV2, config)
    device.deviceInternal.set_status_led_config(config)
end

export get_status_led_config
"""
    $(SIGNATURES)

Returns the configuration as set by :func:`Set Status LED Config`
"""
function get_status_led_config(device::BrickletSolidStateRelayV2)
    return device.deviceInternal.get_status_led_config()
end

export get_chip_temperature
"""
    $(SIGNATURES)

Returns the temperature as measured inside the microcontroller. The
value returned is not the ambient temperature!

The temperature is only proportional to the real temperature and it has bad
accuracy. Practically it is only useful as an indicator for
temperature changes.
"""
function get_chip_temperature(device::BrickletSolidStateRelayV2)
    return device.deviceInternal.get_chip_temperature()
end

export reset
"""
    $(SIGNATURES)

Calling this function will reset the Bricklet. All configurations
will be lost.

After a reset you have to create new device objects,
calling functions on the existing ones will result in
undefined behavior!
"""
function reset(device::BrickletSolidStateRelayV2)
    device.deviceInternal.reset()
end

export write_uid
"""
    $(SIGNATURES)

Writes a new UID into flash. If you want to set a new UID
you have to decode the Base58 encoded UID string into an
integer first.

We recommend that you use Brick Viewer to change the UID.
"""
function write_uid(device::BrickletSolidStateRelayV2, uid)
    device.deviceInternal.write_uid(uid)
end

export read_uid
"""
    $(SIGNATURES)

Returns the current UID as an integer. Encode as
Base58 to get the usual string version.
"""
function read_uid(device::BrickletSolidStateRelayV2)
    return device.deviceInternal.read_uid()
end

export get_identity
"""
    $(SIGNATURES)

Returns the UID, the UID where the Bricklet is connected to,
the position, the hardware and firmware version as well as the
device identifier.

The position can be 'a', 'b', 'c', 'd', 'e', 'f', 'g' or 'h' (Bricklet Port).
A Bricklet connected to an :ref:`Isolator Bricklet <isolator_bricklet>` is always at
position 'z'.

The device identifier numbers can be found :ref:`here <device_identifier>`.
|device_identifier_constant|
"""
function get_identity(device::BrickletSolidStateRelayV2)
    return device.deviceInternal.get_identity()
end

export register_callback
"""
Registers the given *function* with the given *callback_id*.
"""
function register_callback(device::BrickletSolidStateRelayV2, callback_id, function_)
    if isnothing(function_)
        device.registered_callbacks.pop(callback_id, None)
    else
        device.registered_callbacks[callback_id] = function_
    end
end


text:
export SimulatedAmplifier, SimulatedAmplifierParams

Base.@kwdef struct SimulatedAmplifierParams <: DeviceParams
  mode::AmplifierMode = AMP_VOLTAGE_MODE # This should be the safe default
	voltageMode::AmplifierVoltageMode = AMP_LOW_VOLTAGE_MODE # This should be the safe default
	matchingNetwork::Integer = 1
end
SimulatedAmplifierParams(dict::Dict) = params_from_dict(SimulatedAmplifierParams, dict)

Base.@kwdef mutable struct SimulatedAmplifier <: Amplifier
  "Unique device ID for this device as defined in the configuration."
  deviceID::String
  "Parameter struct for this devices read from the configuration."
  params::SimulatedAmplifierParams
  "Flag if the device is optional."
	optional::Bool = false
  "Flag if the device is present."
  present::Bool = false
  "Vector of dependencies for this device."
  dependencies::Dict{String, Union{Device, Missing}}

  state::Bool = false
  mode::AmplifierMode = AMP_VOLTAGE_MODE
  voltageMode::AmplifierVoltageMode = AMP_LOW_VOLTAGE_MODE
  network::Integer = 1
end

function init(amp::SimulatedAmplifier)
  @debug "Initializing simulated amplifier unit with ID `$(amp.deviceID)`."

  # Set values given by configuration
	mode(amp, amp.params.mode)
	voltageMode(amp, amp.params.voltageMode)
	matchingNetwork(amp, amp.params.matchingNetwork)

  amp.present = true
end

checkDependencies(amp::SimulatedAmplifier) = true

Base.close(amp::SimulatedAmplifier) = nothing

state(amp::SimulatedAmplifier) = amp.state
turnOn(amp::SimulatedAmplifier) = amp.state = true
turnOff(amp::SimulatedAmplifier) = amp.state = false
mode(amp::SimulatedAmplifier) = amp.mode
mode(amp::SimulatedAmplifier, mode::AmplifierMode) = amp.mode = mode
voltageMode(amp::SimulatedAmplifier) = amp.voltageMode
voltageMode(amp::SimulatedAmplifier, mode::AmplifierVoltageMode) = amp.voltageMode = mode
matchingNetwork(amp::SimulatedAmplifier) = amp.network
matchingNetwork(amp::SimulatedAmplifier, network::Integer) = amp.network = network
temperature(amp::SimulatedAmplifier) = 25.0u"°C"


text:
# randSeisEvent_well-behaved
function rse_wb(n::Int64)
  Ev = randSeisEvent(n, s=1.0)

  # Prevent problematic :source fields
  Ev.source.eid   = Ev.hdr.id
  Ev.source.npol  = 0
  Ev.source.notes = String[]
  Ev.source.misc  = Dict{String,Any}(
        "pax_desc"    => "azimuth, plunge, length",
        "mt_id"       => "smi:SeisIO/moment_tensor;fmid="*Ev.source.id,
        "planes_desc" => "strike, dip, rake")
  note!(Ev.source, "+origin ¦ " * Ev.source.src)

  # Prevent read/write of problematic :hdr fields
  Ev.hdr.int        = (0x00, "")
  Ev.hdr.loc.datum  = ""
  Ev.hdr.loc.typ    = ""
  Ev.hdr.loc.rms    = 0.0
  Ev.hdr.mag.src    = Ev.hdr.loc.src * ","
  Ev.hdr.notes      = String[]
  Ev.hdr.misc       = Dict{String,Any}()
  note!(Ev.hdr, "+origin ¦ " * Ev.hdr.src)

  # Ensure flags will re-read accurately
  flags = bitstring(Ev.hdr.loc.flags)
  if flags[1] == '1' || flags[2] == '1'
    flags             = "11" * flags[3:8]
    Ev.hdr.loc.flags  = parse(UInt8, flags, base=2)
  end

  # Prevent very low fs, bad locations, rubbish in :misc, :notes
  for j in 1:Ev.data.n
    Ev.data.fs[j]     = max(Ev.data.fs[j], 1.0)
    Ev.data.misc[j]   = Dict{String,Any}()
    Ev.data.notes[j]  = String[]
    Ev.data.loc[j]    = RandSeis.rand_loc(false)

    # Use time structures that formerly broke check_for_gap! / t_extend
    if j < 3
      Ev.data.t[j] = breaking_tstruct(Ev.data.t[j][1,2],
                                      length(Ev.data.x[j]),
                                      Ev.data.fs[j])
    end
  end
  return Ev
end


text:
"""
$(TYPEDEF)

A tag for inplace unitary function. `func(cache, t)` is the actual inplace update function.

# Fields

$(FIELDS)
"""
struct InplaceUnitary
    "inplace update function"
    func::Any
end

isinplace(::Any) = false
isinplace(::InplaceUnitary) = true

"""
$(SIGNATURES)

Calculate the Lamb shift of spectrum `γ` at angular frequency `ω`. All keyword arguments of `quadgk` function is supported.
"""
function lambshift(ω, γ; kwargs...)
    integrand = (x)->(γ(ω+x) - γ(ω-x))/x
    integral, = quadgk(integrand, 0, Inf; kwargs...)
    # TODO: do something with the error information
    -integral / 2 / π
end


text:
module SpringCollab2020AnimatedJumpthruPlatform

using ..Ahorn, Maple

@mapdef Entity "SpringCollab2020/AnimatedJumpthruPlatform" AnimatedJumpthruPlatform(x::Integer, y::Integer, width::Integer=8,
    animationPath::String="SpringCollab2020Example/rainbowwood", animationDelay::Number=0.1, surfaceIndex::Int=-1)

const placements = Ahorn.PlacementDict(
    "Jump Through (Animated) (Spring Collab 2020)" => Ahorn.EntityPlacement(
        AnimatedJumpthruPlatform,
        "rectangle"
    )
)

Ahorn.editingOptions(entity::AnimatedJumpthruPlatform) = Dict{String, Any}(
    "surfaceIndex" => Maple.tileset_sound_ids
)

Ahorn.minimumSize(entity::AnimatedJumpthruPlatform) = 8, 0
Ahorn.resizable(entity::AnimatedJumpthruPlatform) = true, false

function Ahorn.selection(entity::AnimatedJumpthruPlatform)
    x, y = Ahorn.position(entity)
    width = Int(get(entity.data, "width", 8))

    return Ahorn.Rectangle(x, y, width, 8)
end

function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::AnimatedJumpthruPlatform, room::Maple.Room)
    texture = get(entity.data, "animationPath", "SpringCollab2020Example/rainbowwood")

    # Values need to be system specific integer
    x = Int(get(entity.data, "x", 0))
    width = Int(get(entity.data, "width", 8))

    for i in 0:div(width, 8)-1
        Ahorn.drawImage(ctx, "objects/jumpthru/$(texture)00", 8 * i, 0)
    end
end

end


text:
import Pkg
Pkg.add("Documenter")
using Documenter
push!(LOAD_PATH,"../src/")
makedocs(
    #modules=[Home],
    sitename="LaQuiMoBio",
    pages = [
        "Home" => "index.md",
        "Motivation" => "motivation.md",
        "People" => "people.md",
        "Publications" => "publications.md",
        "Patents" => "patents.md",
        "Videos" => "videos.md",
        "Contact" => "contact.md",
    ],
    format = Documenter.HTML(canonical="https://carezende.github.io/home/stable",
                             edit_link=nothing,),

)
deploydocs(
    repo = "github.com/carezende/home.git",
    target = "build",
    branch = "gh-pages",
    versions = ["stable" => "v^", "v#.#" ],
)

max_stars_repo_name: carezende/homepage

max_issues_repo_name: carezende/homepage

max_forks_repo_name: carezende/homepage


text:
#=
By starting at the top of the triangle below and moving to adjacent numbers on the row below, the maximum total from top to bottom is 23.

3
7 4
2 4 6
8 5 9 3

That is, 3 + 7 + 4 + 9 = 23.

Find the maximum total from top to bottom in triangle.txt (right click and 'Save Link/Target As...'), a 15K text file containing a triangle with one-hundred rows.

NOTE: This is a much more difficult version of Problem 18. It is not possible to try every route to solve this problem, as there are 2^99 altogether! If you could check one trillion (1012) routes every second it would take over twenty billion years to check them all. There is an efficient algorithm to solve it. ;o)
=#
#a =

function calc()
  lines = readlines(open(dirname(@__FILE__()) * "/../Resources/p67.txt"))
  arr = map(x -> [parse(Int64, s) for s = split(strip(x), ',')], lines)
  result = Array[]
  push!(result, arr[length(arr)])
  for row in reverse(1:length(arr)-1)
    nextrow = Int64[]
    for column in 1:length(arr[row])
      lowleft = last(result)[column]
      lowright = last(result)[column+1]
      center = arr[row][column]
      push!(nextrow, max(center+lowleft, center+lowright))
    end
    push!(result, nextrow)
  end
  last(result)[1]
end
@time println(calc())

max_stars_repo_name: daniel-beard/JuliaProjectEuler

max_issues_repo_name: daniel-beard/JuliaProjectEuler

max_forks_repo_name: daniel-beard/JuliaProjectEuler


text:
using ExamplePackage_rosa
using Test, Polynomials, Random

@testset "ExamplePackage_rosa.jl" begin
    testfn(x)= (x-1)^3
    testfn_prime(x)=3*(x-1)^2
    testfn1(x)= exp(x)-1
    testfn_prime1(x)=exp(x)
    testfn2(x)= atan(x)
    testfn_prime2(x)= 1/((x)^2+1)

    Random.seed!(1234)

    #test for functions with derivative

    root1, conv, iter  = ExamplePackage_rosa.newtonroot(testfn,testfn_prime,x_0=0.0)
    root2, conv, iter  = ExamplePackage_rosa.newtonroot(testfn1,testfn_prime1,x_0=1.0)
    root3, conv, iter  = ExamplePackage_rosa.newtonroot(testfn1,testfn_prime1,x_0=-0.3)

    @test abs(root1-1.0)<1e-6
    @test abs(root2 - 0.0)<1e-6
    @test abs(root3 - 0.0)<1e-6

    #do the same for functions without derivative

    root1, conv, iter  = ExamplePackage_rosa.newtonroot(testfn,x_0=0.0)
    root2, conv, iter  = ExamplePackage_rosa.newtonroot(testfn1,x_0=1.0)
    root3, conv, iter  = ExamplePackage_rosa.newtonroot(testfn1,x_0=-.4)

    @test abs(root1-1.0)<1e-6
    @test abs(root2 - 0.0)<1e-6
    @test abs(root3 - 0.0)<1e-6

    #check maxiter works
    root1, conv, iter  = ExamplePackage_rosa.newtonroot(testfn,testfn_prime,x_0=0.5,  maxiter = 3)
    # if the function converged, return number of iteration, else return nothing
    @test iter == nothing ||iter <= 3
    #check for tolerance, increasing it to 1e-8
    root1, conv, iter  = ExamplePackage_rosa.newtonroot(testfn,testfn_prime,x_0=0.5,  tol = 1e-8)
    @test abs(root1-1.0)<1e-7

    #check lack iof convergence

    test_nc(x) = x^2+1

    rootnc, conv, iter  = ExamplePackage_rosa.newtonroot(test_nc,x_0=1.0)
    @test rootnc==nothing

    #check bigfloot
    x0 = BigFloat(0.53)
    root1, conv, iter  = ExamplePackage_rosa.newtonroot(testfn,x_0=x0)

    @test abs(root1-1.0)<1e-6

    #check for the polynomials roots

    p_1 = Polynomial([2, -5, 2])
    p_2 = Polynomial(rand(4))
    p_3 = convert(Polynomial,ChebyshevT(rand(4)))

    root1, conv, iter  = ExamplePackage_rosa.newtonroot(p_1,x_0=1.0)
    root2, conv, iter  = ExamplePackage_rosa.newtonroot(p_2,x_0=1.0)
    root3, conv, iter  = ExamplePackage_rosa.newtonroot(p_3,x_0=1.0)

    @test minimum(abs.(roots(p_1).-root1))<=1e-6
    @test minimum(abs.(roots(p_2).-root2))<=1e-6
    @test minimum(abs.(roots(p_3).-root3))<=1e-6

end


text:
function generaVol3d_v4(V_predio, V_bruto, rasante, dcn, dcp)

    numLadosPredio = size(V_predio, 1);
    conjuntoLados = 1:numLadosPredio
    conjuntoLadosCalle = dcp.ladosConCalle;
    numCalles = length(conjuntoLadosCalle)
    conjuntoLadosVecinos = setdiff(conjuntoLados, conjuntoLadosCalle);
    alturaMax = dcn.ALTURAMAX
    sepVecinos = dcn.SEPMIN
    antejardin = dcn.ANTEJARDIN
    anchoEspacioPublico = dcp.ANCHOESPACIOPUBLICO

    # Genera vector de separaciones entre el comienzo de la rasante (separacion predial o eje esp. público)
    # y la línea de edificación
    vecSeparacion = zeros(numLadosPredio)
    for i in conjuntoLadosVecinos
        vecSeparacion[i] = sepVecinos
    end
    numCalles = length(conjuntoLadosCalle)
    for i= 1:numCalles
        vecSeparacion[conjuntoLadosCalle[i]] = anchoEspacioPublico[i]/2 + antejardin
    end
    conjSepDist = 1:.1:alturaMax/rasante
    numSepDistintas = length(conjSepDist)

    # Genera Volumen Teórico a partir de múltiples cortes a distintas alturas
    vec_psVolteor = Array{PolyShape, 1}(undef, numSepDistintas)
    vec_altVolteor = zeros(numSepDistintas,1)
    id = 1
    for j = 1:numSepDistintas # Para cada una de las separaciones distintas
        vecDelta_j = max.(0,vecSeparacion .- conjSepDist[j])
        ps_corte = PolyShape([copy(V_bruto)],1)
        ps_corte = polyShape.polyExpandSides_v2(ps_corte, -vecDelta_j, conjuntoLados)
        ps_corte = polyShape.polyExpand(ps_corte, -conjSepDist[j]);
        if ps_corte.NumRegions >= 1
            vec_psVolteor[id] = ps_corte
            vec_altVolteor[id] = conjSepDist[j]*rasante
            id += 1
        end
    end
    vec_psVolteor = vec_psVolteor[1:id-1]
    vec_psVolteor = [vec_psVolteor[1]; vec_psVolteor]
    vec_altVolteor = vec_altVolteor[1:id-1]
    vec_altVolteor = [0; vec_altVolteor]

    # Determina las alturas en las cuales hay cambios
    vec_psVolteor_ = copy(vec_psVolteor)
    vec_altVolteor_ = copy(vec_altVolteor)
    flagAlt = zeros(length(vec_altVolteor),1)
    flagAlt[1] = 1
    for i = 3:length(vec_altVolteor_)
        n_i = size(vec_psVolteor_[i].Vertices[1],1)
        n_i1 = size(vec_psVolteor_[i-1].Vertices[1],1)
        n_i2 = size(vec_psVolteor_[i-2].Vertices[1],1)
        if n_i < n_i1
            flagAlt[i-1] = 1
            flagAlt[i] = 1
        elseif n_i == n_i1 && n_i1 == n_i2
            pos_i = vec_psVolteor_[i].Vertices[1]
            pos_i1 = vec_psVolteor_[i-1].Vertices[1]
            pos_i2 = vec_psVolteor_[i-2].Vertices[1]
            dist_01 = sqrt.(sum((pos_i .- pos_i1).^2, dims=2))
            dist_12 = sqrt.(sum((pos_i1 .- pos_i2).^2, dims=2))
            for j = 1:length(dist_01)
                if dist_01[j] > 0.01 && dist_12[j] < 0.01
                    flagAlt[i] = 1
                end
            end
        end
    end
    flagAlt[end] = 1
    vec_psVolteor_ = vec_psVolteor_[flagAlt[:] .== 1]
    vec_altVolteor_ = vec_altVolteor_[flagAlt[:] .== 1]

    # Genera arreglo de ids de los vertices de cada nivel
    vecVertices = []
    cont = 0
    for i = 1:size(vec_psVolteor_,1)
        vecAux = zeros(size(vec_psVolteor_[i].Vertices[1],1),1)
        for j = 1:size(vec_psVolteor_[i].Vertices[1],1)
            cont += 1
            vecAux[j] = Int(cont)
        end
        vecVertices = push!(vecVertices, vecAux)
    end

    # Genera matriz de conexiones entre vertices
    matConexionVertices = zeros(size(vec_psVolteor_,1) * size(vec_psVolteor_[1].Vertices[1],1), 2)
    cont = 0
    for i = 1:length(vec_psVolteor_)-1
        V_il = vec_psVolteor_[i].Vertices[1] # matriz Vertices capa low
        V_iu = vec_psVolteor_[i+1].Vertices[1] # matriz Vertices capa up
        for j = 1:size(V_il,1)
            min_dist = 10000
            cont += 1
            if j == 1 || size(V_il,1) != size(V_iu,1) # Para el j inicial y para los casos en que las capas tienen distinto n° de vertices
                for k = 1:size(V_iu,1) # Busca el vertice k (capa up) más cercano al vertice j (capa low)
                    p_l_j = V_il[j,:] # coordenadas vertice j de la capa low
                    p_u_k = V_iu[k,:] # coordenadas vertice k de la capa up
                    dist_jk = sqrt.(sum((p_l_j .- p_u_k).^2))
                    id_l_j = vecVertices[i][j] # id vertice j de la capa low
                    id_u_k = vecVertices[i+1][k] # id vertice k de la capa up
                    if dist_jk < min_dist - 1
                        matConexionVertices[cont,:] = [id_l_j id_u_k]
                        min_dist = dist_jk
                    end
                end
            else #para los casos en que las capas tienen igual n° de vertices
                id_u_a = matConexionVertices[cont-1,2] # id vertice capa up del vertice anterior
                pos_u_a = findall(x->x == id_u_a, vecVertices[i+1])[1][1] # pos capa up del vertice id_u_a
                id_l_j = vecVertices[i][j] # id vertice j de la capa low
                if pos_u_a + 1 <= size(V_iu,1)
                    id_u = vecVertices[i+1][pos_u_a + 1] # id siguiente vertice de la capa up
                else
                    id_u = vecVertices[i+1][1] # id siguiente vertice de la capa up
                end
                matConexionVertices[cont,:] = [id_l_j id_u]
            end
        end
    end
    matConexionVertices = matConexionVertices[matConexionVertices[:,1] .>= 1,:]

    matVolteor = [0 0 0]
    for i = 1:size(vec_psVolteor_,1)
        matVolteor = [matVolteor; vec_psVolteor_[i].Vertices[1] vec_altVolteor_[i]*ones(size(vec_psVolteor_[i].Vertices[1],1),1)]
    end
    matVolteor = matVolteor[2:end, :]
    ps_volteor = PolyShape([matVolteor],1)

    return matConexionVertices, vecVertices, ps_volteor

end

# vec_psVolteor, vec_altVolteor, vec_psVolteor_, vec_altVolteor_


text:
@testset "Sudoku" begin

    @testset "Sudoku MOI/Integers/Int8" begin
        grid = Int8[
            0 2 1 0 7 9 0 8 5
            0 4 5 3 1 0 0 0 9
            0 7 0 0 4 0 0 1 0
            0 0 0 1 0 8 0 3 6
            0 6 0 0 0 0 2 0 8
            0 0 0 0 0 3 0 0 4
            6 0 8 0 0 0 0 0 0
            0 9 4 0 0 7 8 0 0
            2 0 0 5 0 0 0 4 0
        ]
        grid .-= 6
        grid[grid .== 3] .= 10
        grid[grid .== -2] .= 7

        m = CS.Optimizer()
        # don't use 1-9 here but some other integers to test offset and alldifferent without all numbers
        x = [
            [
                MOI.add_constrained_variable(
                    m,
                    CS.Integers([-5, -4, -3, 7, -1, 0, 1, 2, 10]),
                ) for i in 1:9
            ] for j in 1:9
        ]

        # set variables
        for r in 1:9, c in 1:9
            if grid[r, c] != -6
                sat = [MOI.ScalarAffineTerm(Int8(1), x[r][c][1])]
                MOI.add_constraint(
                    m,
                    MOI.ScalarAffineFunction{Int8}(sat, 0),
                    MOI.EqualTo(grid[r, c]),
                )
            end
        end

        # sudoku constraints
        moi_add_sudoku_constr!(m, x)
        MOI.set(m, MOI.TimeLimitSec(), 10.0)

        MOI.optimize!(m)
        @test MOI.get(m, MOI.TerminationStatus()) == MOI.OPTIMAL
        solution = zeros(Int, 9, 9)
        for r in 1:9
            solution[r, :] = [MOI.get(m, MOI.VariablePrimal(), x[r][c][1]) for c in 1:9]
        end
        @test jump_fulfills_sudoku_constr(solution)
        @test CS.get_inner_model(m).options.time_limit == 10.0
    end

    @testset "Hard sudoku with table constraint" begin
        grid = zeros(Int, (9, 9))

        grid[1, :] = [3 8 0 6 0 0 0 0 0]
        grid[2, :] = [0 0 9 0 0 0 0 0 0]
        grid[3, :] = [0 2 0 0 3 0 5 1 0]
        grid[4, :] = [0 0 0 0 0 5 0 0 0]
        grid[5, :] = [0 3 0 0 1 0 0 6 0]
        grid[6, :] = [0 0 0 4 0 0 0 0 0]
        grid[7, :] = [0 1 7 0 5 0 0 8 0]
        grid[8, :] = [0 0 0 0 0 0 9 0 0]
        grid[9, :] = [0 0 0 0 0 7 0 3 2]

        offset = -2
        grid .+= offset

        m = Model(optimizer_with_attributes(CS.Optimizer))
        @variable(m, 1 + offset <= x[1:9, 1:9] <= 9 + offset, Int)
        # set variables
        nvars_set = 0
        for r in 1:9, c in 1:9
            if grid[r, c] != offset
                @constraint(m, x[r, c] == grid[r, c])
                nvars_set += 1
            end
        end

        table = Array{Int64}(undef, (factorial(9), 9))
        i = 1
        for row in permutations(1:9)
            table[i, :] = row .+ offset
            i += 1
        end

        # sudoku constraints
        for rc in 1:9
            @constraint(m, x[rc, :] in CS.TableSet(table))
            @constraint(m, x[:, rc] in CS.TableSet(table))
        end

        for br in 0:2
            for bc in 0:2
                @constraint(
                    m,
                    vec(x[(br * 3 + 1):((br + 1) * 3), (bc * 3 + 1):((bc + 1) * 3)]) in
                    CS.TableSet(table)
                )
            end
        end

        optimize!(m)
        @test JuMP.termination_status(m) == MOI.OPTIMAL
        @test jump_fulfills_sudoku_constr(JuMP.value.(x))
    end

    @testset "Hard sudoku infeasible" begin
        grid = [
            0 0 0 5 4 6 0 0 9
            0 2 0 0 0 0 0 0 7
            0 0 3 9 0 0 0 0 4
            9 0 5 0 0 0 0 7 3
            7 0 0 0 0 0 0 2 0
            0 0 0 0 9 3 0 0 0
            0 5 6 0 0 8 0 0 0
            0 1 0 0 3 9 0 0 0
            0 0 0 0 0 0 8 0 6
        ]

        m = Model(CSJuMPTestOptimizer(; branch_strategy = :ABS))
        @variable(m, 1 <= x[1:9, 1:9] <= 9, Int)
        # set variables
        for r in 1:9, c in 1:9
            if grid[r, c] != 0
                @constraint(m, x[r, c] == grid[r, c])
            end
        end
        # sudoku constraints
        jump_add_sudoku_constr!(m, x)

        optimize!(m)
        @test JuMP.termination_status(m) == MOI.INFEASIBLE
    end

    @testset "Hard fsudoku repo" begin
        grid = zeros(Int, (9, 9))
        grid[1, :] = [0 0 0 0 0 0 0 0 0]
        grid[2, :] = [0 1 0 6 2 0 0 9 0]
        grid[3, :] = [0 0 2 0 0 9 3 1 0]
        grid[4, :] = [0 0 4 0 0 6 0 8 0]
        grid[5, :] = [0 0 8 7 0 2 1 0 0]
        grid[6, :] = [0 3 0 8 0 0 5 0 0]
        grid[7, :] = [0 6 9 1 0 0 4 0 0]
        grid[8, :] = [0 8 0 0 7 3 0 5 0]
        grid[9, :] = [0 0 0 0 0 0 0 0 0]

        m = Model(CSJuMPTestOptimizer())
        @variable(m, 1 <= x[1:9, 1:9] <= 9, Int)
        # set variables
        nvars_set = 0
        for r in 1:9, c in 1:9
            if grid[r, c] != 0
                @constraint(m, x[r, c] == grid[r, c])
                nvars_set += 1
            end
        end

        @test nvars_set == length(filter(n -> n != 0, grid))

        # sudoku constraints
        jump_add_sudoku_constr!(m, x)

        optimize!(m)

        @test JuMP.termination_status(m) == MOI.OPTIMAL
        com = CS.get_inner_model(m)
        @test_reference "refs/hard_fsudoku" test_string([
            constraint.indices for constraint in com.constraints
        ])

        # check that it actually solves the given sudoku
        for r in 1:9, c in 1:9
            if grid[r, c] != 0
                @test JuMP.value(x[r, c]) == grid[r, c]
            end
        end

        @test jump_fulfills_sudoku_constr(JuMP.value.(x))
    end

    @testset "Infeasible sudoku at start" begin
        grid = zeros(Int, (9, 9))
        grid[1, :] = [0 0 0 0 0 0 0 0 0]
        grid[2, :] = [0 1 0 6 2 0 0 9 0]
        grid[3, :] = [0 0 2 0 0 9 3 1 0]
        grid[4, :] = [0 0 4 0 0 6 0 8 0]
        grid[5, :] = [0 0 8 7 0 2 1 0 0]
        grid[6, :] = [0 0 0 0 0 0 5 0 0]
        grid[7, :] = [0 6 9 1 0 0 2 1 0]
        grid[8, :] = [0 0 0 0 7 0 0 5 0]
        grid[9, :] = [0 0 1 0 0 0 0 0 0]

        m = Model(CSJuMPTestOptimizer())
        @variable(m, 1 <= x[1:9, 1:9] <= 9, Int)
        # set variables
        nvars_set = 0
        for r in 1:9, c in 1:9
            if grid[r, c] != 0
                @constraint(m, x[r, c] == grid[r, c])
                nvars_set += 1
            end
        end

        # sudoku constraints
        jump_add_sudoku_constr!(m, x)

        optimize!(m)

        @test JuMP.termination_status(m) == MOI.INFEASIBLE
    end

    @testset "top95 some use backtracking" begin
        grids = sudokus_from_file("data/top95")
        c = 0
        for grid in grids
            m = Model(optimizer_with_attributes(
                CS.Optimizer,
                "solution_type" => Int8,
                "logging" => [],
            ))

            @variable(m, 1 <= x[1:9, 1:9] <= 9, Int)
            # set variables
            for r in 1:9, c in 1:9
                if grid[r, c] != 0
                    @constraint(m, x[r, c] == grid[r, c])
                end
            end

            # sudoku constraints
            jump_add_sudoku_constr!(m, x)
            @objective(m, Min, x[1, 1])

            optimize!(m)
            com = CS.get_inner_model(m)

            @test typeof(com.best_sol) == Int8
            @test JuMP.objective_value(m) == JuMP.value(x[1, 1]) == com.best_sol
            @test JuMP.termination_status(m) == MOI.OPTIMAL
            @test is_solved(com)
            c += 1
        end
        # check that actually all 95 problems were tested
        @test c == 95
    end

    @testset "Number 7 in top95.txt w/o backtracking" begin
        m = Model(optimizer_with_attributes(
            CS.Optimizer,
            "backtrack" => false,
            "logging" => [],
        ))

        grid = Int[
            6 0 2 0 5 0 0 0 0
            0 0 0 0 0 3 0 4 0
            0 0 0 0 0 0 0 0 0
            4 3 0 0 0 8 0 0 0
            0 1 0 0 0 0 2 0 0
            0 0 0 0 0 0 7 0 0
            5 0 0 2 7 0 0 0 0
            0 0 0 0 0 0 0 8 1
            0 0 0 6 0 0 0 0 0
        ]
        grid = transpose(reshape(grid, (9, 9)))

        @variable(m, 1 <= x[1:9, 1:9] <= 9, Int)
        # set variables
        for r in 1:9, c in 1:9
            if grid[r, c] != 0
                @constraint(m, x[r, c] == grid[r, c])
            end
        end

        # sudoku constraints
        jump_add_sudoku_constr!(m, x)

        optimize!(m)

        com = CS.get_inner_model(m)
        @test JuMP.termination_status(m) == MOI.OTHER_LIMIT

        @test !com.info.backtracked
        @test com.info.backtrack_fixes == 0
        @test com.info.in_backtrack_calls == 0
        @show com.info
    end

    @testset "Two solutions" begin
        m = Model(optimizer_with_attributes(
            CS.Optimizer,
            "all_solutions" => true,
            "logging" => [],
        ))

        grid = Int[
            9 0 6 0 7 0 4 0 3
            0 0 0 4 0 0 2 0 0
            0 7 0 0 2 3 0 1 0
            5 0 0 0 0 0 1 0 0
            0 4 0 2 0 8 0 6 0
            0 0 3 0 0 0 0 0 5
            0 3 0 7 0 0 0 5 0
            0 0 7 0 0 5 0 0 0
            4 0 5 0 1 0 7 0 8
        ]

        @variable(m, 1 <= x[1:9, 1:9] <= 9, Int)
        # set variables
        for r in 1:9, c in 1:9
            if grid[r, c] != 0
                @constraint(m, x[r, c] == grid[r, c])
            end
        end

        # sudoku constraints
        jump_add_sudoku_constr!(m, x)

        optimize!(m)

        com = CS.get_inner_model(m)

        @test JuMP.result_count(m) == 2
        @test JuMP.value.(x) != JuMP.value.(x, result = 2)
        if JuMP.value.(x[6, 5:6]) == [9, 4]
            @test JuMP.value.(x[7, 5:6]) == [4, 9]
        else
            @test JuMP.value.(x[6, 5:6]) == [4, 9]
            @test JuMP.value.(x[7, 5:6]) == [9, 4]
        end
    end

    @testset "Two solutions but optimize on one value" begin
        # all solutions should give all solutions even if they are not optimal
        m = Model(optimizer_with_attributes(
            CS.Optimizer,
            "all_solutions" => true,
            "logging" => [],
        ))

        grid = Int[
            9 0 6 0 7 0 4 0 3
            0 0 0 4 0 0 2 0 0
            0 7 0 0 2 3 0 1 0
            5 0 0 0 0 0 1 0 0
            0 4 0 2 0 8 0 6 0
            0 0 3 0 0 0 0 0 5
            0 3 0 7 0 0 0 5 0
            0 0 7 0 0 5 0 0 0
            4 0 5 0 1 0 7 0 8
        ]

        @variable(m, 1 <= x[1:9, 1:9] <= 9, Int)
        # set variables
        for r in 1:9, c in 1:9
            if grid[r, c] != 0
                @constraint(m, x[r, c] == grid[r, c])
            end
        end

        # sudoku constraints
        jump_add_sudoku_constr!(m, x)
        @objective(m, Max, x[6, 5])

        optimize!(m)

        com = CS.get_inner_model(m)

        @test JuMP.result_count(m) == 2
        @test JuMP.value.(x) != JuMP.value.(x, result = 2)
        # the better one should be solution 1 of course
        @test JuMP.objective_value(m) == 9
        @test JuMP.value.(x[6, 5:6]) == [9, 4]
        @test JuMP.value.(x[7, 5:6]) == [4, 9]
        # second worse solution
        @test JuMP.objective_value(m, result = 2) == 4
        @test JuMP.value.(x[6, 5:6], result = 2) == [4, 9]
        @test JuMP.value.(x[7, 5:6], result = 2) == [9, 4]
    end

    @testset "Two solutions but only one optimal" begin
        # all solutions should give all solutions even if they are not optimal
        m = Model(optimizer_with_attributes(
            CS.Optimizer,
            "all_optimal_solutions" => true,
            "logging" => [],
        ))

        grid = Int[
            9 0 6 0 7 0 4 0 3
            0 0 0 4 0 0 2 0 0
            0 7 0 0 2 3 0 1 0
            5 0 0 0 0 0 1 0 0
            0 4 0 2 0 8 0 6 0
            0 0 3 0 0 0 0 0 5
            0 3 0 7 0 0 0 5 0
            0 0 7 0 0 5 0 0 0
            4 0 5 0 1 0 7 0 8
        ]

        @variable(m, 1 <= x[1:9, 1:9] <= 9, Int)
        # set variables
        for r in 1:9, c in 1:9
            if grid[r, c] != 0
                @constraint(m, x[r, c] == grid[r, c])
            end
        end

        # sudoku constraints
        jump_add_sudoku_constr!(m, x)
        @objective(m, Max, x[6, 5])

        optimize!(m)

        com = CS.get_inner_model(m)

        @test JuMP.result_count(m) == 1
        # the better one should be solution 1 of course
        @test JuMP.objective_value(m) == 9
        @test JuMP.value.(x[6, 5:6]) == [9, 4]
        @test JuMP.value.(x[7, 5:6]) == [4, 9]
    end

    @testset "Two solutions but both optimal" begin
        # all solutions should give all solutions even if they are not optimal
        m = Model(optimizer_with_attributes(
            CS.Optimizer,
            "all_optimal_solutions" => true,
            "logging" => [],
        ))

        grid = Int[
            9 0 6 0 7 0 4 0 3
            0 0 0 4 0 0 2 0 0
            0 7 0 0 2 3 0 1 0
            5 0 0 0 0 0 1 0 0
            0 4 0 2 0 8 0 6 0
            0 0 3 0 0 0 0 0 5
            0 3 0 7 0 0 0 5 0
            0 0 7 0 0 5 0 0 0
            4 0 5 0 1 0 7 0 8
        ]

        @variable(m, x[1:9, 1:9], CS.Integers(1:9))
        # set variables
        for r in 1:9, c in 1:9
            if grid[r, c] != 0
                @constraint(m, x[r, c] == grid[r, c])
            end
        end

        # sudoku constraints
        jump_add_sudoku_constr!(m, x)
        @objective(m, Max, x[1, 1])

        optimize!(m)

        com = CS.get_inner_model(m)
        @test com.info.n_constraint_types.alldifferent == 27
        @test length(com.constraints) == 27
        @test length(com.search_space) == 81

        # the better one should be solution 1 of course
        @test JuMP.objective_value(m) == 9
        @test JuMP.result_count(m) == 2
        @test JuMP.value.(x) != JuMP.value.(x, result = 2)
        if JuMP.value.(x[6, 5:6]) == [9, 4]
            @test JuMP.value.(x[7, 5:6]) == [4, 9]
        else
            @test JuMP.value.(x[6, 5:6]) == [4, 9]
            @test JuMP.value.(x[7, 5:6]) == [9, 4]
        end
    end

    @testset "All optimal solutions" begin
        n = 9
        g = 3

        grid = zeros(Int, (n, n))

        m = Model(optimizer_with_attributes(
            CS.Optimizer,
            "all_optimal_solutions" => true,
            "time_limit" => 1.0,
            "keep_logs" => true,
            "logging" => [],
        ))

        @variable(m, 1 <= x[1:n, 1:n] <= n, Int)
        # set variables
        for r in 1:n, c in 1:n
            if grid[r, c] != 0
                @constraint(m, x[r, c] == grid[r, c])
            end
        end
        for rc in 1:n
            @constraint(m, x[rc, :] in CS.AllDifferentSet())
            @constraint(m, x[:, rc] in CS.AllDifferentSet())
        end
        for br in 0:(g - 1)
            for bc in 0:(g - 1)
                @constraint(
                    m,
                    vec(x[(br * g + 1):((br + 1) * g), (bc * g + 1):((bc + 1) * g)]) in
                    CS.AllDifferentSet()
                )
            end
        end

        optimize!(m)
        @test MOI.get(m, MOI.SolveTime()) >= 1.0
        # at least more than 1 but in that time frame it should find a lot ;)
        @test MOI.get(m, MOI.ResultCount()) >= 10
        @test JuMP.termination_status(m) == MOI.TIME_LIMIT
        com = CS.get_inner_model(m)
        general_tree_test(com)
        logs = CS.get_logs(com)
        @test CS.sanity_check_log(logs[:tree])
    end

    @testset "All solutions" begin
        n = 9
        g = 3

        grid = zeros(Int, (n, n))

        m = Model(optimizer_with_attributes(
            CS.Optimizer,
            "all_solutions" => true,
            "time_limit" => 1.0,
            "logging" => [],
        ))

        @variable(m, 1 <= x[1:n, 1:n] <= n, Int)
        # set variables
        for r in 1:n, c in 1:n
            if grid[r, c] != 0
                @constraint(m, x[r, c] == grid[r, c])
            end
        end
        for rc in 1:n
            @constraint(m, x[rc, :] in CS.AllDifferentSet())
            @constraint(m, x[:, rc] in CS.AllDifferentSet())
        end
        for br in 0:(g - 1)
            for bc in 0:(g - 1)
                @constraint(
                    m,
                    vec(x[(br * g + 1):((br + 1) * g), (bc * g + 1):((bc + 1) * g)]) in
                    CS.AllDifferentSet()
                )
            end
        end

        optimize!(m)
        @test MOI.get(m, MOI.SolveTime()) >= 1.0
        # at least more than 1 but in that time frame it should find a lot ;)
        @test MOI.get(m, MOI.ResultCount()) >= 10
        @test JuMP.termination_status(m) == MOI.TIME_LIMIT
    end
end


text:
@testset "GeoMean" begin
    mock.optimize! = (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, [ones(4); 2; √2; √2])
    bridged_mock = MOIB.GeoMean{Float64}(mock)
    MOIT.geomean1vtest(bridged_mock, config)
    MOIT.geomean1ftest(bridged_mock, config)
    # Dual is not yet implemented for GeoMean bridge
    ci = first(MOI.get(bridged_mock, MOI.ListOfConstraintIndices{MOI.VectorAffineFunction{Float64}, MOI.GeometricMeanCone}()))
    test_delete_bridge(bridged_mock, ci, 4, ((MOI.VectorAffineFunction{Float64}, MOI.RotatedSecondOrderCone, 0),
                                            (MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64},      1)))
end


text:
"save_attributes: Saves attributes of Stochastic Solutions to HDF5 file."
function Solutions.save_attributes(sol::StochasticSolution, h5::HDF5.File)
    attributes(h5)["ntime"] = ntime(sol)
    attributes(h5)["nsave"] = nsave(sol)
    attributes(h5)["conv"]  = string(conv(sol))
    attributes(h5)["nd"] = sol.nd
    attributes(h5)["nm"] = sol.nm
    attributes(h5)["ns"] = sol.ns
    attributes(h5)["K"]  = sol.K
end

function Solutions.init_solution_in_hdf5(solution::SolutionSDE{DT,TT,2}, h5::HDF5.File) where {DT,TT}
    q = create_dataset(h5, "q", DT, ((solution.nd, solution.nt+1), (solution.nd, -1)), chunk=(solution.nd,1))
    q[:,1] = solution.q[:,0]
end

function Solutions.init_solution_in_hdf5(solution::SolutionSDE{DT,TT,3}, h5::HDF5.File) where {DT,TT}
    q = create_dataset(h5, "q", DT, ((solution.nd, solution.nt+1, solution.ns),(solution.nd, -1, solution.ns)), chunk=(solution.nd,1,1))
    q[:,1,:] = solution.q[:,0,:]
end

function Solutions.init_solution_in_hdf5(solution::SolutionPSDE{DT,TT,2}, h5::HDF5.File) where {DT,TT}
    q = create_dataset(h5, "q", DT, ((solution.nd, solution.nt+1), (solution.nd, -1)), chunk=(solution.nd,1))
    p = create_dataset(h5, "p", DT, ((solution.nd, solution.nt+1), (solution.nd, -1)), chunk=(solution.nd,1))
    q[:, 1] = solution.q[:, 0]
    p[:, 1] = solution.p[:, 0]
end

function Solutions.init_solution_in_hdf5(solution::SolutionPSDE{DT,TT,3}, h5::HDF5.File) where {DT,TT}
    q = create_dataset(h5, "q", DT, ((solution.nd, solution.nt+1, solution.ns),(solution.nd, -1, solution.ns)), chunk=(solution.nd,1,1))
    p = create_dataset(h5, "p", DT, ((solution.nd, solution.nt+1, solution.ns),(solution.nd, -1, solution.ns)), chunk=(solution.nd,1,1))
    q[:, 1, :] = solution.q[:, 0, :]
    p[:, 1, :] = solution.p[:, 0, :]
end

function init_increments_in_hdf5(solution::StochasticSolution{DT,TT,NQ,2}, h5::HDF5.File) where {DT,TT,NQ}
    dW = create_dataset(h5, "ΔW", DT, ((solution.nm, solution.ntime),(solution.nm, -1)), chunk=(solution.nm,1))
    dZ = create_dataset(h5, "ΔZ", DT, ((solution.nm, solution.ntime),(solution.nm, -1)), chunk=(solution.nm,1))
end

function init_increments_in_hdf5(solution::StochasticSolution{DT,TT,NQ,3}, h5::HDF5.File) where {DT,TT,NQ}
    dW = create_dataset(h5, "ΔW", DT, ((solution.nm, solution.ntime, solution.ns),(solution.nm, -1, solution.ns)), chunk=(solution.nm,1,1))
    dZ = create_dataset(h5, "ΔZ", DT, ((solution.nm, solution.ntime, solution.ns),(solution.nm, -1, solution.ns)), chunk=(solution.nm,1,1))
end

"Creates HDF5 file and initialises datasets for stochastic solution object."
function Solutions.create_hdf5(solution::StochasticSolution, file::AbstractString; save_W=true)
    # create HDF5 file and save attributes and common parameters
    h5 = createHDF5(solution, file)

    # save attributes
    save_attributes(solution, h5)

    # create dataset
    init_timeteps_in_hdf5(solution, h5)
    init_solution_in_hdf5(solution, h5)

    if save_W
        init_increments_in_hdf5(solution, h5)
    end

    return h5
end

function Solutions.copy_solution_to_hdf5(solution::SolutionSDE{DT,TT,1}, h5::HDF5.File, j1, j2, n1, n2) where {DT <: Number, TT}
    if size(h5["q"])[end] < j2
        HDF5.set_extent_dims(h5["q"], (j2,))
    end
    h5["q"][j1:j2] = solution.q[n1:n2]
end

function Solutions.copy_solution_to_hdf5(solution::SolutionSDE{AT,TT,1}, h5::HDF5.File, j1, j2, n1, n2) where {DT, AT <: Array{DT}, TT}
    elaxes = axes(solution.q[begin])
    if size(h5["q"])[end] < j2
        HDF5.set_extent_dims(h5["q"], (size(h5["q"])[1:end-1]..., j2))
    end
    for i in eachindex(j1:j2, n1:n2)
        j = (j1:j2)[i]
        n = (n1:n2)[i]
        h5["q"][elaxes..., j] = solution.q[n]
    end
end

function Solutions.copy_solution_to_hdf5(solution::SolutionSDE{DT,TT,2}, h5::HDF5.File, j1, j2, n1, n2) where {DT <: Number, TT}
    if size(h5["q"],1) < j2
        HDF5.set_extent_dims(h5["q"], (j2, size(h5["q"],2)))
    end
    h5["q"][j1:j2, :] = solution.q[n1:n2, :]
end

function Solutions.copy_solution_to_hdf5(solution::SolutionSDE{AT,TT,2}, h5::HDF5.File, j1, j2, n1, n2) where {DT, AT <: Array{DT}, TT}
    elaxes = axes(solution.q[begin,begin])
    if size(h5["q"])[end-1] < j2
        HDF5.set_extent_dims(h5["q"], (size(h5["q"])[begin:end-2]..., j2, size(h5["q"])[end]))
    end
    for k in 1:nsamples(solution.q)
        for i in eachindex(j1:j2, n1:n2)
            j = (j1:j2)[i]
            n = (n1:n2)[i]
            h5["q"][elaxes..., j, k] = solution.q[n, k]
        end
    end
end

function Solutions.copy_solution_to_hdf5(solution::SolutionPSDE{DT,TT,1}, h5::HDF5.File, j1, j2, n1, n2) where {DT <: Number, TT}
    if size(h5["q"])[end] < j2
        HDF5.set_extent_dims(h5["q"], (j2,))
    end
    if size(h5["p"])[end] < j2
        HDF5.set_extent_dims(h5["p"], (j2,))
    end
    h5["q"][j1:j2] = solution.q[n1:n2]
    h5["p"][j1:j2] = solution.p[n1:n2]
end

function Solutions.copy_solution_to_hdf5(solution::SolutionPSDE{AT,TT,1}, h5::HDF5.File, j1, j2, n1, n2) where {DT, AT <: Array{DT}, TT}
    elaxes = axes(solution.q[begin])
    if size(h5["q"])[end] < j2
        HDF5.set_extent_dims(h5["q"], (size(h5["q"])[1:end-1]..., j2))
    end
    for i in eachindex(j1:j2, n1:n2)
        j = (j1:j2)[i]
        n = (n1:n2)[i]
        h5["q"][elaxes..., j] = solution.q[n]
    end

    elaxes = axes(solution.p[begin])
    if size(h5["p"])[end] < j2
        HDF5.set_extent_dims(h5["p"], (size(h5["p"])[1:end-1]..., j2))
    end
    for i in eachindex(j1:j2, n1:n2)
        j = (j1:j2)[i]
        n = (n1:n2)[i]
        h5["p"][elaxes..., j] = solution.p[n]
    end
end

function Solutions.copy_solution_to_hdf5(solution::SolutionPSDE{DT,TT,2}, h5::HDF5.File, j1, j2, n1, n2) where {DT <: Number, TT}
    if size(h5["p"],1) < j2
        HDF5.set_extent_dims(h5["p"], (j2, size(h5["p"],2)))
    end
    if size(h5["p"],1) < j2
        HDF5.set_extent_dims(h5["p"], (j2, size(h5["p"],2)))
    end
    h5["q"][j1:j2, :] = solution.q[n1:n2, :]
    h5["p"][j1:j2, :] = solution.p[n1:n2, :]
end

function Solutions.copy_solution_to_hdf5(solution::SolutionPSDE{AT,TT,2}, h5::HDF5.File, j1, j2, n1, n2) where {DT, AT <: Array{DT}, TT}
    elaxes = axes(solution.q[begin,begin])
    if size(h5["q"])[end-1] < j2
        HDF5.set_extent_dims(h5["q"], (size(h5["q"])[begin:end-2]..., j2, size(h5["q"])[end]))
    end
    for k in 1:nsamples(solution.q)
        for i in eachindex(j1:j2, n1:n2)
            j = (j1:j2)[i]
            n = (n1:n2)[i]
            h5["q"][elaxes..., j, k] = solution.q[n, k]
        end
    end

    elaxes = axes(solution.p[begin,begin])
    if size(h5["p"])[end-1] < j2
        HDF5.set_extent_dims(h5["p"], (size(h5["p"])[begin:end-2]..., j2, size(h5["p"])[end]))
    end
    for k in 1:nsamples(solution.p)
        for i in eachindex(j1:j2, n1:n2)
            j = (j1:j2)[i]
            n = (n1:n2)[i]
            h5["p"][elaxes..., j, k] = solution.p[n, k]
        end
    end
end

function copy_increments_to_hdf5(solution::StochasticSolution{DT,TT,NQ,2}, h5::HDF5.File, j1, j2, n1, n2) where {DT,TT,NQ}
    if size(h5["ΔW"],2) < j2
        HDF5.set_extent_dims(h5["ΔW"], (size(h5["ΔW"],1), j2))
    end
    if size(h5["ΔZ"],2) < j2
        HDF5.set_extent_dims(h5["ΔZ"], (size(h5["ΔZ"],1), j2))
    end
    h5["ΔW"][:,j1:j2] = solution.W.ΔW[:,n1:n2]
    h5["ΔZ"][:,j1:j2] = solution.W.ΔZ[:,n1:n2]
end

function copy_increments_to_hdf5(solution::StochasticSolution{DT,TT,NQ,3}, h5::HDF5.File, j1, j2, n1, n2) where {DT,TT,NQ}
    if size(h5["ΔW"],2) < j2
        HDF5.set_extent_dims(h5["ΔW"], (size(h5["ΔW"],1), j2, size(h5["ΔW"],3)))
    end
    if size(h5["ΔZ"],2) < j2
        HDF5.set_extent_dims(h5["ΔZ"], (size(h5["ΔZ"],1), j2, size(h5["ΔZ"],3)))
    end
    h5["ΔW"][:,j1:j2,:] = solution.W.ΔW[:,n1:n2,:]
    h5["ΔZ"][:,j1:j2,:] = solution.W.ΔZ[:,n1:n2,:]
end

"""
Append solution to HDF5 file.
  soffset - start writing the solution q at the position soffset+2
  woffset - start writing the increments ΔW, ΔZ at the position woffset+1
"""
function Common.write_to_hdf5(solution::StochasticSolution, h5::HDF5.File=hdf5(solution), soffset=offset(solution), woffset=ioffset(solution))
    # set convenience variables and compute ranges
    js1 = soffset+2
    js2 = soffset+1+solution.nt
    jw1 = woffset+1
    jw2 = woffset+solution.nwrite

    # copy data from solution to HDF5 dataset
    copy_timeteps_to_hdf5(solution, h5, js1, js2, 1, solution.nt)
    copy_solution_to_hdf5(solution, h5, js1, js2, 1, solution.nt)

    if exists(h5, "ΔW") && exists(h5, "ΔZ")
        # copy the Wiener process increments from solution to HDF5 dataset
        copy_increments_to_hdf5(solution, h5, jw1, jw2, 1, solution.nwrite)
    end

    return nothing
end


text:
function linedist(xold::Vector{Float64},fold::Float64,g::Vector{Float64},p::Vector{Float64},
    stpmax::Float64,f::Function,J::Function,maxiter::Int16,V::Float64,Q::Float64,
    Vs::Float64,vs::Float64,ts::Float64,V0::Float64,rho::Float64,beta::Float64,C::Float64,
    W1end::Float64,c0::Float64,A0::Float64,h::Float64)
    alpha = 1e-4;
    tolx = 1e-7;
    check = false;

    # scale if attempted step too large
    np = norm(p);
    # println(p)
    if np > stpmax
        p = p*stpmax/np;
    end
    # println(p)

    slope = dot(g,p);
    if slope >= 0
        print(g)
        print(p)
        error("Roundoff problem in line search.")
    end

    # compute minimum step fraction
    test = 0.;
    for i = 1:length(xold)
        temp = abs(p[i])/max(abs(xold[i]),1.);
        if temp > test
            test = temp;
        end
    end
    alamin = tolx/test;

    omega = 1;
    alam = omega; # start w/ desired fraction of Newton step
    N = 1;

    while N < maxiter
        x = xold+alam*p;
        JJ = J(x,Vs,vs,ts,C,W1end,c0,rho,beta,h);
        # println(JJ)
        D = diagm(maximum!(zeros(length(x)),abs.(JJ)).^-1);
        fvec = D*f(x,V,Q,Vs,vs,ts,V0,rho,beta,C,W1end,c0,A0,h);
        # println(D)
        # println(fvec)
        fn = 0.5*dot(fvec,fvec);
        if alam < alamin
            x = xold;
            # println(alam)
            # println(alamin)
            # println(fn)
            # println(fold+alpha*alam*slope)
            # println("Δx converged in distal line search. Verify in Newton loop.")
            check = true;
            return fn,x,check
        elseif fn < fold+alpha*alam*slope
            # println(fn)
            # println(fold)
            # println(alpha*alam*slope)
            return fn,x,check
        else
            if alam == omega
                tmplam = -slope/(2*(fn-fold-slope));
            else
                rhs1 = fn-fold-alam*slope
                rhs2 = fn2-fold-alam2*slope
                a = (rhs1/alam^2-rhs2/alam2^2)/(alam-alam2);
                b = (-alam2*rhs1/alam^2+alam*rhs2/alam2^2)/(alam-alam2);
                if a == 0
                    tmplam = -slope/(2*b);
                else
                    disc = b^2-3*a*slope;
                    if disc < 0
                        tmplam = 0.5*alam;
                    elseif b <= 0
                        tmplam = (-b+sqrt(disc)/(3*a));
                    else
                        tmplam = -slope/(b+sqrt(disc));
                    end
                end
                if tmplam > 0.5*alam
                    tmplam = 0.5*alam;
                end
            end
        end
        alam2 = alam;
        fn2 = fn;
        alam = max(tmplam,0.1*alam);
        N+=1;
        if N == maxiter
            println(xn)
            println(f(xn,V,Q,Vs,vs,ts,V0,rho,beta,C,W1end,c0,A0,h))
            error("Distal line search iteration failed to converge.");
        end
    end
end

max_stars_repo_name: dcanuto/CardioModeling

max_issues_repo_name: dcanuto/CardioModeling

max_forks_repo_name: dcanuto/CardioModeling


text:
abstract type AbstractEvent end

_undefined_error(meth, T, args) = error("Method`$(meth)(e::$T, $args)` not defined!")

## ------------------------------------------------------------------
old_state(e::AbstractEvent, args...) = error("Method `old_state(e::$(typeof(e)), args...)` not defined!")
new_state(e::AbstractEvent, args...) = error("Method `new_state(e::$(typeof(e)), args...)` not defined!")

"""
This should ensure that `old_state` returns an object that do not `trigger` the handler over the return of `new_state`.
That is, it configure the handler to react to a 'incoming event'.
"""
update!

update!(e::AbstractEvent, args...) = error("Method `update!(e::$(typeof(e)), args...)` not defined!")

reset!(e::AbstractEvent, args...) = error("Method `reset!(e::$(typeof(e)), args...)` not defined!")
trigger(e::AbstractEvent, old_state, new_state) = error("Method `trigger(e::$(typeof(e)), old_state, new_state)` not defined!")
states(e::AbstractEvent) = error("Method `states(e::$(typeof(e)))::AbstractDict` not defined!")
istraking(e::AbstractEvent, key) = haskey(states(e), key)

## ------------------------------------------------------------------
Base.getindex(e::AbstractEvent, key) = getindex(states(e), key)
Base.setindex!(e::AbstractEvent, value, key) = setindex!(states(e), value, key)
Base.get!(e::AbstractEvent, key, default) = get!(states(e), key, default)
Base.get!(f::Function, e::AbstractEvent, key) = get!(f, states(e), key)
Base.get(e::AbstractEvent, key, default) = get(states(e), key, default)
Base.get(f::Function, e::AbstractEvent, key) = get(f, states(e), key)

## ------------------------------------------------------------------
"""
Returns true if the event happened.
It do not update the handler state (see bang version).
"""
has_event

has_event(e::AbstractEvent, args...) = trigger(e, old_state(e, args...), new_state(e, args...))

"""
Returns true if the event happened.
Additionally, it will call `update!` over the handler (only if `has_event` triggered).
"""
has_event!

function has_event!(e::AbstractEvent, args...)
    flag = has_event(e, args...)
    flag && update!(e, args...)
    return flag
end

"""
Executes `f()` if `has_event` returns true.
Returns the results of `f()` or nothing.
"""
on_event

on_event(f::Function, e::AbstractEvent, args...) = (has_event(e, args...) ? f() : nothing)

"""
Executes `f()` if `has_event` returns true.
Additionally, it will call `update!` over the handler (only if `has_event` triggered).
Returns the results of `f()` or nothing.
"""
on_event!

function on_event!(f::Function, e::AbstractEvent, args...)
    flag = has_event!(e, args...)
    return flag ? f() : nothing
end


text:
using WordleSolver
using Test

@testset "WordleSolver.jl" begin
    # Write your tests here.
end


text:
using Random
using Plots
pyplot(size = (600,1100), legend = true)
"""

Arguments:
    - `f`: Function to be optimized
    - `g`: Gradient function for `f`
    - `c`: Constraints function. Evaluates all the constraint equations and return a Vector of values
    - `x0`: (Vector) Initial position to start from
    - `n`: (Int) Number of evaluations allowed. Remember `g` costs twice of `f`
    - `prob`: (String) Name of the problem. So you can use a different strategy for each problem
"""
function optimize(f, g, c, x0, n, prob)
    p = squared_sum_penalty_function(c)
    x_best, xlist,fpvalues = penalty_method(f, p, x0, n)
    return x_best
end

function squared_sum_penalty_function(c)
    f = function(x)
        ceval = c(x)
        sum(ceval[ceval .> 0].^2)
    end
    return f
end

function penalty_method(f, p, x, n; ρ=1, γ=2)
    numEvalsLeft = n
    hj_max_eval = 10;
    xlist = Array{Float64}[]
    fpvalues = Array{Float64}[]
    push!(xlist, x)
    while numEvalsLeft >= hj_max_eval
        x, evalsUsed, newFeasFound, fpvalues = hooke_jeeves(f, ρ, p, x, hj_max_eval, fpvalues)
        push!(xlist, x)
        numEvalsLeft = numEvalsLeft - evalsUsed
        # if !newFeasFound
        #     ρ /= γ
        # else
            ρ *= γ
        # end
    end
    return x,xlist,fpvalues
end

function last_feasible_x(currXFeas, x, feasFound, last_x_feas)
    if currXFeas
        return x
    else
        if !feasFound
            return x
        else
            return last_x_feas
        end
    end
end

basis(i, n) = [k == i ? 1.0 : 0.0 for k in 1 : n]

function hooke_jeeves(f, ρ, p, x, maxEval, fpvalues; α=0.2, γ=0.1)
    xf = f(x)
    xp = p(x)
    push!(fpvalues,[xf,xp])
    y, m = xf + ρ*xp, length(x)
    evalsUsed = 1
    last_x_feas = x
    newFeasFound = false
    if xp==0
        feasFound = true
        currXFeas = true
    else
        feasFound = false
        currXFeas = false
    end
    while evalsUsed < maxEval
        improved = false
        for i in 1 : m
            x′ = x + α*basis(i, m)
            if evalsUsed >= maxEval
                return last_feasible_x(currXFeas, x, feasFound, last_x_feas), evalsUsed, newFeasFound, fpvalues
            end
            xf′ = f(x′)
            xp′ = p(x′)
            push!(fpvalues,[xf′,xp′])
            y′ = xf′ + ρ*xp′
            evalsUsed = evalsUsed + 1
            if y′ < y
                x, y, improved = x′, y′, true
                if xp′==0
                    last_x_feas = x
                    newFeasFound = true
                    feasFound = true
                    currXFeas = true
                else
                    currXFeas = false
                end
            else
                x′ = x - α*basis(i, m)
                if evalsUsed >= maxEval
                    return last_feasible_x(currXFeas, x, feasFound, last_x_feas), evalsUsed, newFeasFound, fpvalues
                end
                xf′ = f(x′)
                xp′ = p(x′)
                push!(fpvalues,[xf′,xp′])
                y′ = xf′ + ρ*xp′
                evalsUsed = evalsUsed + 1
                if y′ < y
                    x, y, improved = x′, y′, true
                    if xp′==0
                        last_x_feas = x
                        newFeasFound = true
                        feasFound = true
                        currXFeas = true
                    else
                        currXFeas = false
                    end
                end
            end
        end
        if !improved
            α *= γ
        end
    end
    return last_x_feas, evalsUsed, newFeasFound, fpvalues
end

function grad(x::Vector)
    return 0
end

function simple1(x::Vector)
    return -x[1] * x[2]
end

function simple1_constraints(x::Vector)
    return [x[1] + x[2]^2 - 1,
            -x[1] - x[2]]
end

function simple1_init()
    return rand(2) * 2.0
end

function simple2(x::Vector)
    return (1.0 - x[1])^2 + 100.0 * (x[2] - x[1]^2)^2
end

function simple2_constraints(x::Vector)
    return [(x[1]-1)^3 - x[2] + 1,
            x[1] + x[2] - 2]
end

function simple2_init()
    return rand(2) .* 2.0 .- 1.0
end

function simple3(x::Vector)
    return x[1] - 2*x[2] + x[3]
end

function simple3_constraints(x::Vector)
    return [x[1]^2 + x[2]^2 + x[3]^2 - 1]
end

function simple3_init()
    b = 2.0 .* [1.0, -1.0, 0.0]
    a = -2.0 .* [1.0, -1.0, 0.0]
    return rand(3) .* (b-a) + a
end

function get_full_list(f, g, c, x0, n, prob)
    p = squared_sum_penalty_function(c)
    x_best, xlist, fpvalues = penalty_method(f, p, x0, n)
    return x_best, xlist, fpvalues
end

function create_path_plot(problemNumber, fn, constr, numConstr, lvls, initial_points, true_min, xMin, xMax, yMin, yMax, lineColors, contourColors)
    pyplot(size = (600,1000), legend = true)
    # Plot function contours
    x = xMin:0.1:xMax
    y = yMin:0.1:yMax
    f(x, y) = begin
            fn([x,y])
        end
    p1 = contour(x, y, f, fill=false, seriescolor=:grays,w=3, label="function")

    # Plot constraint contours
    for i in 1:numConstr
        c(x, y) = begin
            simple1_constraints([x,y])[i]
        end
        scolor = contourColors[i]
        contour!(x, y, c, seriescolor=scolor, w=2,  levels=lvls)
    end

    # Plot true minimum
    scatter!([true_min[1]],[true_min[2]],markercolor=:green,markersize=8,label="true minimum")

    for i in 1:size(initial_points,1)
        x0 = initial_points[i]
        lineColor = lineColors[i]

        # Plot initial point
        scatter!([x0[1]],[x0[2]],marker=:d,markersize=5,markercolor=lineColor,label="initial point")

        # Plot trajectory and final point
        x_best, xlist, _  = get_full_list(simple1, grad, simple1_constraints, x0, 1000, "simple_1")
        xnlist=hcat(xlist...)'
        plot!(xnlist[:,1],xnlist[:,2], color=lineColor, marker=:dot, label="path of feasible points", w=2)
        scatter!([x_best[1]],[x_best[2]],marker=:d,markersize=5,markercolor=:red, label="best found")
    end
    title!("Simple $problemNumber path taken")
    return p1
end

problemNumber = 1
fn = simple1
constr = simple1_constraints
numConstr = 2
lvls = [-2, -1, 0]
initial_points1 = [[1.51183, 1.63151],
 [0.447248, 0.233549],
 [1.09013, 0.345318]]
# for i in 1:3
#     initial_points[i] = simple1_init()
# end
true_min = [2.0/3.0; 1.0/sqrt(3)]
xMin = -2
xMax = 2
yMin = -2
yMax = 2
lineColors = [:indigo, :orchid, :lightpink]
contourColors = [:blues, :reds, :greens]

p1 = create_path_plot(problemNumber, fn, constr, numConstr, lvls, initial_points1, true_min, xMin, xMax, yMin, yMax, lineColors, contourColors)

problemNumber = 2
fn = simple2
constr = simple2_constraints
lvls = [-3, -2, -1, 0]
initial_points2 = [[-0.303123, 0.733811],
 [-0.807351, -0.710301],
 [-0.450702, 0.209268]]
# for i in 1:3
#     initial_points[i] = simple2_init()
# end
true_min = [1;1]
p2 = create_path_plot(problemNumber, fn, constr, numConstr, lvls, initial_points2, true_min, xMin, xMax, yMin, yMax, lineColors, contourColors)

plt1 = plot(p1,p2,layout=(2,1))
gui(plt1)

function create_convergence_plot(fpvalues, problemNumber)
    pyplot(size = (1000,1000), legend = true)
    fp=hcat(fpvalues...)'
    colors = [:green, :red]
    color = fp[2,2] <= 0 ? colors[1] : colors[2]
    firstRed = false
    firstGreen = false
    p1 = plot(fp[:,1], label= "", color = :black)
    # p2=plot([1],[fp[1,1]],marker=:circle, markersize=5, color=(fp[1,2] <= 0 ? colors[1] : colors[2]), label = "")
    for i in 1:size(fp,1)
        lbl = ""
        if !firstRed && fp[i,2] > 0
            lbl = "infeasible"
            firstRed = true
        elseif !firstGreen && fp[i,2] <= 0
            lbl = "feasible"
            firstGreen = true
        end

        plot!([i],[fp[i,1]],marker=:circle, markersize=5, color=(fp[i,2] <= 0 ? colors[1] : colors[2]), label = lbl)
    end
    # plot(p2)
    xlabel!("Iteration number")
    ylabel!("Value of f(x)")
    title!("Simple $problemNumber convergence")
    return p1
end

# xbest, _, fpvalues  = get_full_list(simple1, grad, simple1_constraints, initial_points1[1], 1000, "simple_1")
# problemNumber = 1
# p1 = create_convergence_plot(fpvalues, problemNumber)

# xbest, _, fpvalues  = get_full_list(simple2, grad, simple2_constraints, initial_points2[1], 1000, "simple_2")
# problemNumber = 2
# p2 = create_convergence_plot(fpvalues, problemNumber)

# init3 = [-1.9022234854723372,
#  -1.3978361756095108,
#   0.0]
# xbest, _, fpvalues  = get_full_list(simple3, grad, simple1_constraints, init3, 1000, "simple_3")
# problemNumber = 3
# p3 = create_convergence_plot(fpvalues, problemNumber)

# plt1 = plot(p1,p2,p3,layout=(3,1))
# gui(plt1)


text:
"""
    ⊂(X::LazySet{N}, Y::LazySet, [witness]::Bool=false) where {N}

Strict inclusion check.

### Input

- `X`       -- first set
- `Y`       -- second set
- `witness` -- (optional, default: `false`) compute a witness if activated

### Output

* If `witness` option is deactivated: `true` iff ``X ⊂ Y``
* If `witness` option is activated:
  * `(true, v)` iff ``X ⊂ Y`` and ``v ∈ Y \\setminus X``
  * `(false, [])` iff not ``X ⊂ Y``

### Algorithm

We check inclusion of `X` in `Y` and then check inclusion of `Y` in `X`:

```math
X ⊂ Y \\Leftrightarrow X ⊆ Y \\land ¬ (Y ⊆ X)
```
"""
function ⊂(X::LazySet{N}, Y::LazySet, witness::Bool=false) where {N}
    if witness
        res, w = ⊆(X, Y, witness)
        if res
            res, w = ⊆(Y, X, witness)
            if res
                return (false, N[])
            else
                return (true, w)
            end
        else
            return (res, N[])
        end
    end
    return (X ⊆ Y) && !(Y ⊆ X)
end


text:
using JLD2
using Statistics

num_bins = 12

mutable struct OneDStateVisitation <: LoggerKeyData
    state_visitations::Array{Int32,3}
    temp_state_visitations::Array{Int32, 2}
    log_interval::Int

    function OneDStateVisitation(logger_init_info)
        num_logged_steps = fld(logger_init_info[LoggerInitKey.TOTAL_STEPS], logger_init_info[LoggerInitKey.INTERVAL]) + 1
        new(zeros(num_bins, num_bins, num_logged_steps), zeros(num_bins, num_bins), logger_init_info[LoggerInitKey.INTERVAL])
    end
end

function lg_start!(self::OneDStateVisitation, env, agent)
end

function lg_step!(self::OneDStateVisitation, env, agent, s, a, s_next, r, is_terminal, cur_step_in_episode, cur_step_total)
    x = ceil(Int32, s_next[1] * num_bins)
    y = ceil(Int32, s_next[2] * num_bins)
    if (x == 0)
        x = 1
    end
    if (y == 0)
        y = 1
    end
    if (s_next[1] == 0.49999999999999994)
        println()
        println("found state")
        println("state: $(s[1:2])")
        println("next state: $(s_next[1:2])")
        print("================")
    end
    if (y == 3 && x != 6)
        println()
        println("BAD STATE")
        println("state: $(s[1:2])")
        println("next state: $(s_next[1:2])")
    end

    # println("state: $(s[1]), $(s[2])")
    # println("$(x), $(y)")

    # println(s)
    # print(x, y)

    self.temp_state_visitations[x, y] += 1
    if (rem(cur_step_total, self.log_interval) == 0)
        ind = fld(cur_step_total, self.log_interval) + 1
        self.state_visitations[:, :, ind] = self.temp_state_visitations

        # Reset temp state Visitation
        # println(sum(self.state_visitations))
        self.temp_state_visitations .*= 0
    end
end

function lg_episode_end!(self::OneDStateVisitation, cur_step_in_episode, cur_step_total)
end

function save_log(self::OneDStateVisitation, save_dict::Dict)
    save_dict[:oned_tmaze_state_visitation] = self.state_visitations
end

max_stars_repo_name: MatthewMcLeod/curiosity

max_issues_repo_name: MatthewMcLeod/curiosity

max_forks_repo_name: MatthewMcLeod/curiosity


text:
"Parametric representation of expressions"
module ParamExpression

"""
A parametric representation of an expression

This encodes a function from some known primitives types to an expression.

```
θ::Bool -> Θ ? (+ 1 1) : (- 1 1)
```
```
"""
struct ParamExpr
end

end


text:
# License is MIT: LICENSE.md

include("setup.jl")
include("strtest.jl")
include("bench.jl")

@time books = load_books()
@time save_books(books)
@time tst = checktests(test=true)


text:
# Test for gausshermite()

tol = 1e-14

n = 18;
x,w = gausshermite( n )
@test isa(x,Vector{Float64})
@test isa(w,Vector{Float64})
@test (length(x) == n && length(w) == n)
@test (dot(w,x) < tol && abs(dot(w,x.^2) - sqrt(pi)/2) < tol)

# Test a small n:
n = 42
x,w = gausshermite( n )
@test isa(x,Vector{Float64})
@test isa(w,Vector{Float64})
@test (length(x) == n && length(w) == n)
@test (dot(w,x) < tol && abs(dot(w,x.^2) - sqrt(pi)/2) < tol)
@test abs(x[37] - 5.660357581283058) < tol
@test abs(w[17] - 0.032202101288908) < tol

# Test a larger n:
n = 251
x,w = gausshermite( n )
@test isa(x,Vector{Float64})
@test isa(w,Vector{Float64})
@test (length(x) == n && length(w) == n)
@test (dot(w,x) < tol && abs(dot(w,x.^2) - sqrt(pi)/2) < 300*tol)
@test abs(x[37] - -13.292221459334638) < tol
@test abs(w[123] - 0.117419270715955) < 2*tol


text:
@testset "Gibbs Sampling" begin
    seed!(42)
    L = 3
    D = 10
    N = 20
    nSamples = 10
    b = 5
    i = GibbsSampling(; nBurnin=0)
    @test repr(i) == "Gibbs Sampler"

    @test AGP.ρ(i) == 1
    @test AGP.is_stochastic(i) == false
end


text:
"""
Package: Forecast

    arsim(Φ,Φ0,x0,n; Σ,E,fix)
    arsim(AR,n;fix)

Simulate a multivariate autoregressive series model.

The simulated series follows the model:

```math
Xt = \\Phi_0 + \\sum_{i=1}^p \\Phi_i \\cdot X_{t-i} + E
```

# Arguments
- `AR`:           AR struct coming from an `ar` model.
- `Φ`:            Array with dimensions (m,m,p) for the parameters in the AR model.
- `Φ0`:           Vector size `m` for the constant in the AR model. Default value is 0.
- `x0`:           Array with dimensions (m,p) for the initial value in the AR model. Default value is a random value from zero to one.
- `n`:            Number of simulations.

- `fix`:          Matrix{Union{Missing,Float64}} containing values to be fixed in th simulation.
- `Σ2`:           Variance Covariance matrix for the AR model with a MvNormal distribution for the noise. Default value is an identity Matrix.
- `E`:            Distribution for the error.

# Returns
A multivariate series simulating an AR model each column containing a dimension and ordered by time ascending rows.

# Examples
```julia-repl
julia> arsim(1,1,1,10)
10-element Vector{Float64,1}:
[...]
"""
function arsim(xar::AR, n::Integer,
               fix::Union{Nothing,AbstractVector{Union{Missing,T}}} = nothing) where T<:Real
    Φ = compact(xar.Φ)
    Φ0 = compact(xar.Φ0)
    x0 = compact(zeros(size(Φ,2),size(Φ,3)))
    Σ2 = compact(xar.Σ2)
    arsim(Φ,Φ0,x0,n;Σ2,fix)
end

function arsim(Φ::T, Φ0::T, x0::T, n::Integer;
               Σ2::T = T(1), E::Union{Nothing,Distribution} = Normal(T(0),T(1)),
               fix::Union{Nothing,AbstractVector{Union{Missing,T}}} = nothing) where T<:Real

    x = Array{T}(undef,n)

    if Σ2 == T(0) && isnothing(fix)
        for i in 1:n
            x[i] = x0 = Φ * x0 + Φ0
        end
        return x
    end

    if Σ2 == T(0) && !isnothing(fix)
        @assert size(x) == size(fix)
        for i in 1:n
            x[i] = x0 = ismissing(fix[i]) ? Φ * x0 + Φ0 : fix[i]
        end
        return x
    end

    E =  Σ2 != T(1) ? Normal(0,sqrt(Σ2)) : E
    e = rand(E,n)

    for i in 1:n
        x[i] = x0 = Φ * x0 + Φ0 + e[i]
    end

    return x

end

function arsim(Φ::AbstractVector{T}, Φ0::T, x0::AbstractVector{T}, n::Integer;
               Σ2::T = T(1), E::Union{Nothing,Distribution} = Normal(T(0),T(1)),
               fix::Union{Nothing,AbstractVector{Union{Missing,T}}} = nothing) where T<:Real

    np = length(Φ)

    @assert length(Φ) == length(x0)

    x = Array{T}(undef,n)

    if Σ2 == T(0) && isnothing(fix)
        for i in 1:n
            x[i] = sum(Φ .* x0) + Φ0
            x0 = vcat(x[i],x0)[1:np]
        end
        return x
    end

    if Σ2 == T(0) && !isnothing(fix)
        @assert size(x) == size(fix)
        for i in 1:n
            x[i] = ismissing(fix[i]) ? sum(Φ .* x0) + Φ0 : fix[i]
            x0 = vcat(x[i],x0)[1:np]
        end
        return x
    end

    E =  Σ2 != T(1) ? Normal(0,sqrt(Σ2)) : E
    e = rand(E,n)

    for i in 1:n
        x[i] = sum(Φ .* x0) + Φ0 + e[i]
        x0 = vcat(x[i],x0)[1:np]
    end

    return x

end

function arsim(Φ::AbstractArray{T}, Φ0::AbstractVector{T}, x0::AbstractArray{T}, n::Integer;
               Σ2::Matrix{T} = collect(T(1)*I(length(Φ0))),
               E::Union{Nothing,Distribution} = MvNormal(length(Φ0),1),
               fix::Union{Nothing,AbstractMatrix{Union{Missing,T}}} = nothing) where T<:Real

    m = size(Φ,1)
    np = size(Φ,3)

    @assert length(Φ0) == m
    @assert size(x0,1) == m
    @assert size(x0,2) == np

    mnp = m*np
    x0 = reshape(x0,mnp)
    Φ = reshape(Φ,m,mnp)

    x = Array{T}(undef,n,m)

    if Σ2 == zeros(m,m) && isnothing(fix)
        for i in 1:n
            x[i,:] = Φ * x0 + Φ0
            x0 = vcat(x[i,:],x0)[1:mnp]
        end

        return x
    end

    if Σ2 == zeros(m,m) && !isnothing(fix)
        @assert size(x) == size(fix)
        for i in 1:n
            x[i,:] =  Φ * x0 + Φ0
            fixi = findall((!ismissing).(fix[i,:]))
            x[i,fixi] = fix[i,fixi]
            x0 = vcat(x[i,:],x0)[1:mnp]
        end
        return x
    end

    E =  Σ2 != collect(I(length(Φ0))) ? MvNormal(Σ2) : E
    e = rand(E,n)

    for i in 1:n
        x[i,:] = Φ * x0 + Φ0 + e[:,i]
        x0 = vcat(x[i,:],x0)[1:mnp]
    end

    return x

end


text:
const DEFAULT_BUF_SIZE = UInt64(1024 * 1024 * 30)
const DEFAULT_CHUNK_SIZE = UInt64(1024 * 1024 * 4)

mutable struct SequenceBuffer{T<:ReadDatastore}
    data_store::T
    bufferdata::Vector{UInt8}
    chunk_size::Int
    buffer_position::Int
end

function SequenceBuffer(ds::T) where {T<:ReadDatastore}
    return SequenceBuffer{T}(ds, Vector{UInt8}(undef, DEFAULT_BUF_SIZE), DEFAULT_CHUNK_SIZE, typemax(Int))
end

@inline chunksize(sb::SequenceBuffer) = sb.chunk_size
@inline buffersize(sb::SequenceBuffer) = length(bufferdata(sb))
@inline bufferpos(sb::SequenceBuffer) = sb.buffer_position
@inline bufferpos!(sb::SequenceBuffer, pos) = sb.buffer_position = pos
@inline datastore(sb::SequenceBuffer) = sb.data_store
@inline stream(sb::SequenceBuffer) = stream(datastore(sb))
@inline bufferdata(sb::SequenceBuffer) = sb.bufferdata

"""
    _load_sequence_data!(seq::LongSequence, buffer::Vector{UInt8}, offset::Integer)

This internal method reads the packed data of a sequence from a blob of bytes.
"""
function _load_sequence_data!(seq::LongSequence, buffer::Vector{UInt8}, offset::Integer)
    seqdata = BioSequences.encoded_data(seq)
    for i in eachindex(seqdata)
        seqdata[i] = unsafe_load(convert(Ptr{UInt64}, pointer(buffer, offset + 1)))
        offset = offset + sizeof(UInt64)
    end
    return seq
end

function _load_sequence(data::Vector{UInt8}, offset::Integer)
    sequence_length = unsafe_load(convert(Ptr{UInt64}, pointer(data, offset + 1)))
    offset = offset + sizeof(UInt64)
    seq = LongSequence{DNAAlphabet{4}}(sequence_length)
    return _load_sequence_data!(seq, data, offset)
end

Base.eltype(sb::SequenceBuffer) = Base.eltype(datastore(sb))
Base.firstindex(sb::SequenceBuffer) = firstindex(datastore(sb))
Base.lastindex(sb::SequenceBuffer) = lastindex(datastore(sb))
Base.length(sb::SequenceBuffer) = length(datastore(sb))
Base.eachindex(sb::SequenceBuffer) = eachindex(datastore(sb))
Base.IteratorSize(sb::SequenceBuffer) = Base.IteratorSize(datastore(sb))
Base.IteratorEltype(sb::SequenceBuffer) = Base.IteratorEltype(datastore(sb))

@inline function Base.getindex(sb::SequenceBuffer, idx::Integer)
    @boundscheck checkbounds(datastore(sb), idx)

    # Get the position in the io stream that the read data begins.
    roif = unsafe_read_offset_in_file(datastore(sb), idx)

    # Knowing the position in the iostream, check if the buffer will need
    # refreshing.
    if roif < bufferpos(sb) || (roif + chunksize(sb)) > bufferpos(sb) + buffersize(sb)
        bufferpos!(sb, roif)
        seek(stream(sb), roif)
        readbytes!(stream(sb), bufferdata(sb), buffersize(sb))
    end

    return _load_sequence(bufferdata(sb), roif - bufferpos(sb))
end

@inline function Base.iterate(sb::SequenceBuffer, state = 1)
    @inbounds if firstindex(sb) ≤ state ≤ lastindex(sb)
        return sb[state], state + 1
    else
        return nothing
    end
end


text:
abstract type AbstractHydroFormulation end

abstract type HydroCurtailmentForm <: AbstractHydroFormulation end

include("hydro_generation/hydro_variables.jl")
include("hydro_generation/curtailment_constraints.jl")


text:
function construct(start::Vector{FP.Polynomial{T}}, target::Vector{FP.Polynomial{T}}
    ) where {T<:Number}
    T, start, target
end

function construct(::Type{T}, start::Vector{FP.Polynomial{T}}, target::Vector{FP.Polynomial{T}}
    ) where {T<:Number}
    start, target
end

function construct(::Type{T}, start::Vector{<:FP.Polynomial}, target::Vector{<:FP.Polynomial}
    ) where {T<:Number}
    convert(Vector{FP.Polynomial{T}}, start), convert(Vector{FP.Polynomial{T}}, target)
end

function construct(start::Vector{FP.Polynomial{U}}, target::Vector{FP.Polynomial{V}}
    ) where {U<:Number, V<:Number}
    P = promote_type(U, V, Float64)
    P, convert(Vector{FP.Polynomial{P}}, start), convert(Vector{FP.Polynomial{P}}, target)
end

function construct(::Type{T}, start::Vector{<:MP.AbstractPolynomial}, target::Vector{<:MP.AbstractPolynomial}
    ) where {T<:Number}
    convert(Vector{FP.Polynomial{T}}, start), convert(Vector{FP.Polynomial{T}}, target)
end

function construct(start::Vector{<:MP.AbstractPolynomial{U}}, target::Vector{<:MP.AbstractPolynomial{V}}
    ) where {U<:Number, V<:Number}
    P = promote_type(U, V, Float64)
    P, convert(Vector{FP.Polynomial{P}}, start), convert(Vector{FP.Polynomial{P}}, target)
end

function construct(::Type{T}, start::MP.AbstractPolynomial, target::MP.AbstractPolynomial
    ) where {T<:Number}
    s, t = convert(Vector{FP.Polynomial{T}}, [start, target])
    [s], [t]
end

function construct(start::MP.AbstractPolynomial{U}, target::MP.AbstractPolynomial{V}
    ) where {U<:Number, V<:Number}
    P = promote_type(U, V, Float64)
    s, t = convert(Vector{FP.Polynomial{P}}, [start, target])
    P, [s], [t]
end


text:
using JuMP
using CSV
using DataFrames
#using AmplNLWriter

data_mass = dropmissing(CSV.read("plecak.csv"))
data_plan = CSV.read("plan.csv")

plan_summary = by(data_plan, [:Przedmiot], plan_summary -> size(plan_summary, 1))
mass_summary = by(data_mass, [:Przedmiot], mass_summary -> sum(mass_summary[:Waga]))

przybory_mass = mass_summary[mass_summary[:Przedmiot] .== "Przybory", :][:x1] / 1000

mass_summary = mass_summary[mass_summary[:Przedmiot] .!= "Przybory", :]

dict_plan = Dict(zip(plan_summary[:Przedmiot], plan_summary[:x1]))
dict_mass = Dict(zip(mass_summary[:Przedmiot], mass_summary[:x1]))

data_in= DataFrame(Przedmiot = String[],
    Sztuk = Int64[],
    Waga = Float64[])

for row in 1:size(plan_summary, 1)
    course = plan_summary[:Przedmiot][row]
    push!(data_in,[course, get(dict_plan, course, 0), get(dict_mass, course, 0) / 1000])
end
println(data_in)

course_no = size(data_in,1)
max_course_occurrences = maximum(data_in[:Sztuk])
days = 5
max_per_day = 8
min_per_day = 3

avg_mass = sum(data_in[:Sztuk].*data_in[:Waga]) / days + przybory_mass[1]

course_map = Dict(zip(1:course_no, data_in[:Przedmiot]))
mass_map = Dict(zip(1:course_no, data_in[:Waga]))

max_occurences_per_day  = 2
if max_occurences_per_day < max_course_occurrences
    max_course_occurrences = max_occurences_per_day
end

print("Masa przyborów: ")
println(przybory_mass[1])


text:
export BinningVectorObservable, push!, tau, reset!
export extrapolate_tau, extrapolate_stderror
export BinningVectorObservableSet

mutable struct BinningVectorObservable <: VectorObservable
    ## index of these vectors denotes the level of bins (each bin stores the mean of 2^(i-1) values)
    raw_ts :: Vector{Vector{Float64}}        ## Time series of raw data
    bins :: Vector{Vector{Float64}}          ## Time series of bins
    sum :: Vector{Vector{Float64}}           ## Summation of stored values
    sum2 :: Vector{Vector{Float64}}          ## Summation of square of bin's mean
    entries :: Vector{Int}           ## Number of bins
    binsize :: Int
    lastbin :: Int
    minbinnum :: Int
    maxlevel :: Int
end

function BinningVectorObservable(minbinnum::Int = 128)
    BinningVectorObservable( Vector{Float64}[], Vector{Float64}[], Vector{Float64}[], Vector{Float64}[],
                             zeros(Int,1), 1, 0, minbinnum, 1)
end

function reset!(b::BinningVectorObservable)
    b.raw_ts = Vector{Float64}[]
    b.bins = Vector{Float64}[]
    b.sum = Vector{Float64}[]
    b.sum2 = Vector{Float64}[]
    b.entries = zeros(Int,1)
    b.binsize = 1
    b.lastbin = 0
    b.maxlevel = 1

    return b
end

maxlevel(b::BinningVectorObservable) = b.maxlevel
count(b::BinningVectorObservable, level::Int = 1) = b.entries[level]
sum(b::BinningVectorObservable, level::Int = 1) = b.sum[level]
sum2(b::BinningVectorObservable, level::Int = 1) = b.sum2[level]

function push!(b::BinningVectorObservable, x::Vector)
    ## time series
    push!(b.raw_ts, deepcopy(x))
    if b.lastbin == 0
        push!(b.bins, deepcopy(x))
        b.lastbin += 1
    else
        b.bins[end] += x
        b.lastbin += 1
    end
    if b.lastbin == b.binsize
        b.bins[end] /= b.binsize
        b.lastbin = 0
        if length(b.bins) == 2(b.minbinnum)
            new_bins = Vector{Vector{Float64}}(b.minbinnum)
            for i in 1:b.minbinnum
                new_bins[i] = 0.5(b.bins[2i-1]+b.bins[2i])
            end
            b.bins = new_bins
            b.binsize <<= 1
            b.maxlevel += 1
        end
    end

    ## binning
    if b.entries[1] == 0
        push!(b.sum, deepcopy(x))
        push!(b.sum2, squared(x))
    else
        b.sum[1] += x
        b.sum2[1] += squared(x)
    end
    i = b.entries[1]
    b.entries[1] += 1
    level = 2
    bsize = 2
    while i & 1 > 0
        if level > length(b.sum)
            push!(b.sum, zeros(length(x)))
            push!(b.sum2,zeros(length(x)))
            push!(b.entries, 0)
        end
        lastbin = b.sum[1]/bsize - b.sum[level]
        if b.entries[level] == 0
            b.sum[level] = lastbin
            b.sum2[level] = squared(lastbin)
            b.entries[level] += 1
        else
            b.sum[level] += lastbin
            b.sum2[level] += squared(lastbin)
            b.entries[level] += 1
        end

        level += 1
        bsize <<= 1
        i >>= 1
    end

    return b
end

function mean(b::BinningVectorObservable, level::Int = 1)
    return sum(b, level)./count(b, level)
end

function var(b::BinningVectorObservable, level::Int = 1)
    n = count(b, level)
    s = sum(b, level)
    s2 = sum2(b, level)
    if n > 1
        v2 = s2 .- (s.^2)./n
        return map(maxzero, v2/(n-1))
    else
        return fill(NaN, length(b.sum[1]))
    end
end

function stderror(b::BinningVectorObservable, level::Int = maxlevel(b))
    return sqrt.(var(b,level)./count(b,level))
end

function confidence_interval(b::BinningVectorObservable, confidence_rate::Real, level::Int = maxlevel(b))
    q = 0.5+0.5*confidence_rate
    correction = quantile( TDist(count(b,level)), q)
    serr = stderror(b, level)
    return correction * serr
end
function confidence_interval(b::BinningVectorObservable, confidence_rate_symbol::Symbol = :sigma1, level::Int = maxlevel(b))
    n = parsesigma(confidence_rate_symbol)
    return confidence_interval(b, erf(0.5n*sqrt(2.0)), level)
end

function tau(b::BinningVectorObservable, level::Int = maxlevel(b))
    binsize = 1<<(level-1)
    return 0.5.*( (binsize.*var(b,level))./var(b) .- 1.0)
end

#=
linearmodel(x::Float64, p::Vector{Float64}) = p[1] + x*p[2]
linearmodel(xs::Vector{Float64}, p::Vector{Float64}) = map(x->linearmodel(x,p),xs)
=#

function extrapolate_detail(op :: Function, b::BinningVectorObservable, point::Int)
    ml = maxlevel(b)
    ll = max(ml-point+1, 1)
    levels = ll:ml
    ns = map( level->1<<(level-1), levels)
    ninvs = ns .\ 1.0
    ys = map( level->op(b, level), levels)
    fit = curve_fit(linearmodel, ninvs, ys, [ys[end], 0.0])
    return fit.param[1], estimate_errors(fit)[1]
end
extrapolate_tau(b::BinningVectorObservable, point::Int = 5) = extrapolate_detail(tau, b, point)
extrapolate_stderror(b::BinningVectorObservable, point::Int = 5) = extrapolate_detail(stderror, b, point)

function show(io::IO, obs::BinningVectorObservable)
    if count(obs) > 0
        print(io, mean(obs), " +/- ", stderror(obs), "; tau = ", tau(obs))
    else
        print(io, "No entries")
    end
end

const BinningVectorObservableSet = MCObservableSet{BinningVectorObservable}


text:
export
Line,
    set!

mutable struct Line{T <: AbstractFloat}
    p1::Point{T}
    p2::Point{T}
    # Transformed vertices
    bucket::Array{Point{Float64},1}

    function Line{T}() where {T <: AbstractFloat}
        new(Point{T}(), Point{T}(), [])
    end

    function Line{T}(x1::T, y1::T, x2::T, y2::T) where {T <: AbstractFloat}
        new(Point{T}(x1, y1), Point{T}(x2, y2))
    end
end

# copy specializations
copy(p::Line) = Line{Float64}(p.x1, p.y1, p.x2, p.y2)

# setters/getters set_point
function set!(l::Line{T}, x1::T, y1::T, x2::T, y2::T) where {T <: AbstractFloat}
    set!(l.p1, x1, y1)
    set!(l.p2, x2, y2);
end

function set!(l::Line{T}, from::Line{T}) where {T <: AbstractFloat}
    set!(l.p1, from.x1, from.y1)
    set!(l.p2, from.x2, from.y2);
end

max_stars_repo_name: wdevore/Ranger-Julia-SDL

max_issues_repo_name: wdevore/Ranger-Julia-SDL

max_forks_repo_name: wdevore/Ranger-Julia-SDL


text:
using GeometricFlux, Flux, JLD2, SparseArrays, DiffEqFlux, DifferentialEquations
using Flux: onehotbatch, onecold, logitcrossentropy, throttle
using Flux: @epochs
using Statistics: mean
using CUDA

# Load the dataset
@load "data/cora_features.jld2" features
@load "data/cora_labels.jld2" labels
@load "data/cora_graph.jld2" g

# Model and Data Configuration
num_nodes = 2708
num_features = 1433
hidden = 16
target_catg = 7
epochs = 40

# Preprocess the data and compute adjacency matrix
train_X = Matrix{Float32}(features) |> gpu  # dim: num_features * num_nodes
train_y = Float32.(labels) |> gpu  # dim: target_catg * num_nodes
fg = FeaturedGraph(g) |> gpu

# Define the Neural GDE
diffeqarray_to_array(x) = reshape(gpu(x), size(x)[1:2])

node = NeuralODE(
    GCNConv(fg, hidden=>hidden) |> gpu,
    (0.f0, 1.f0), Tsit5(), save_everystep = false,
    reltol = 1e-3, abstol = 1e-3, save_start = false
)

model = Chain(GCNConv(fg, num_features=>hidden, relu),
              Dropout(0.5),
              node,
              diffeqarray_to_array,
              GCNConv(fg, hidden=>target_catg),
              softmax) |> gpu

# Loss
loss(x, y) = logitcrossentropy(model(x), y)
accuracy(x, y) = mean(onecold(model(x)) .== onecold(y))

# Training
## Model Parameters
ps = Flux.params(model, node.p);
## Training Data
train_data = [(train_X, train_y)]
## Optimizer
opt = ADAM(0.05)
## Callback Function for printing accuracies
evalcb() = @show(accuracy(train_X, train_y))

## Training Loop
@epochs epochs Flux.train!(loss, ps, train_data, opt, cb=throttle(evalcb, 10))