text:
module TMatrices

export LS,
    ODE

module LS
    export SetupMatrix,
        Quadrature

    include("quadrature.jl")
    include("LS.jl")
end

module ODE
    include("ODE.jl")
end

"""
SquareWell(R ; width=1., height=0.1)

Sample potential, given by V = `height` for R < `width` and V = 0 otherwise.
"""
SquareWell(R ; width=1., height=0.1) = (R < width ? height : 0.)

"""
Gaussian(R ; width=1., height=0.1)

Sample potential, given by a Gaussian centred at R=0 with a width of `width` and
scaled by `height`.
"""
Gaussian(R ; width=1., height=0.1) = height*exp(-1/2 * (R/width)^2)

"""
TestCase(R)

A sample potential which includes a repulsive barrier and an attractive part, to
mimic an atomic system.

The scattering length of this potential seems to be around 4.65 a0.
"""
TestCase(R) = SquareWell(R, width=0.5, height=10.0) + Gaussian(R, width=2., height=-0.5)

"""
TestCase2(R)

A sample potential which includes a repulsive barrier and an attractive part, to
mimic an atomic system.

The scattering length of this potential seems to be around -34 a0.
"""
TestCase2(R) = SquareWell(R, width=1.0, height=100.0) + Gaussian(R, width=2., height=-2.0)

end # module


text:
"""
ДАНО: Робот - в произвольной клетке ограниченного прямоугольного поля

РЕЗУЛЬТАТ: Робот - в исходном положении, и все клетки поля промакированы
"""

using HorizonSideRobots
originRobot = Robot(11,11,animate = false)

# moving_in_angle! передвинет робота в Юго-Западный угол поля и вернёт в переменную back_bath ОБРАТНЫЙ путь робота
function moving_in_angle!(robot::Robot)
    back_path = []
    while (!(isborder(robot, West)) || !(isborder(robot, Sud)))
        for side in (Sud, West)
            steps = 0
            while (!(isborder(robot,side)))
                move!(robot,side)
                steps += 1
            end
            push!(back_path, steps)
        end
    end
    return reverse!(back_path)
end

# moving_back_to_start! вернёт робота в начальное положение через обратный путь робота
function moving_back_to_start!(robot::Robot, back_path::Vector)::Nothing
    for (index,value) in enumerate(back_path)
        for step in 1:value
            move!(robot, HorizonSide((((index + 1)% 2) + 3) % 4))
        end
    end
end
# get_lines промаркирует всё поле змейкой
function get_lines!(robot::Robot)
    flag = true
    while (flag)
        for side in [Nord,Sud]
            putmarker!(robot)
            while (!(isborder(robot,side)))
                move!(robot, side)
                putmarker!(robot)
            end
            if (!(isborder(robot,Ost)))
                move!(robot,Ost)
            else
                flag = false
                break
            end
        end
    end
end

# Point3_master! главная функция программы
function master3!(robot::Robot)
    # moving_in_angle! передвинет робота в Юго-Западный угол поля и вернёт в переменную back_bath ОБРАТНЫЙ путь робота
    back_path = moving_in_angle!(robot)

    putmarker!(robot)
    # get_lines промаркирует всё поле змейкой
    get_lines!(robot)

    for side in [Sud,West]
        while (!(isborder(robot,side)))
            move!(robot,side)
        end
    end

    # moving_back_to_start! вернёт робота в начальное положение через обратный путь робота
    moving_back_to_start!(robot, back_path)

end

master3!(originRobot)
show!(originRobot)

ext: jl

max_stars_repo_name: Igetsuga/AlgorithmicsHomework

max_issues_repo_name: Igetsuga/AlgorithmicsHomework

max_forks_repo_name: Igetsuga/AlgorithmicsHomework


text:
#

#
# Horner initially stores a list of coefficients of a polynomial, a evaluation point,
# and optionally a corresponding list of absolute errors of the coefficients.
# On demand it modifies the values in order to obtain the values of the selected derivatives
# at the given point and the corresponding error value.
#

mutable struct Horner{S<:Number,T<:Number,U<:AbstractFloat,V<:AbstractFloat}
    x::S							# evaluation spot
    coeff::AbstractVector{T}		# polynomial coefficients
    devx::U							# absolute error estimatin for x
    devcoeff::AbstractVector{V}		# absolute error estimatins for coefficients
    valid::Int						# evaluation status
end

# Create Horner table for evaluation at point x.
# Provide polynomial coefficients and optionally absolute errors for x and the coefficients.
function Horner(x::S, coeff::AbstractVector{T}, devx::U, devcoeff::AbstractVector{V}) where {S<:Number,T<:Number,U<:AbstractFloat,V<:AbstractFloat}

    if devcoeff == []
        devcoeff = abs(coeff) * eps(promote_type(T,Float64))
    else
        devcoeff = copy(devcoeff)
    end
    if devx < 0
        devx = abs(x) * eps(real(S))
    end
    if length(coeff) != length(devcoeff) throw(ArgumentError("devcoeff length wrong")) end
    R = promote_type(S,T)
	Horner(x, map(R,coeff), devx, devcoeff, -1)
end

# Evaluate Horner table up to and including m deviations.
# return the m'th deviation of the polynomial and the error estimation
# If m >= degree, all deviations ot the polynomial have been calculated at x.
# Further calls just return the cached results.
function call(h::Horner{S,T,U,V}, m::Int) where {S,T,U,V}
#function (h::Horner{S,T,U,V})(m::Int)
    if m < 0 throw(ArgumentError("index $m is not >= 0")) end
	n = length(h.coeff) - 1
    if m > h.valid
        if h.x != 0
            ep = eps(real(T))
            b, e = h.coeff, h.devcoeff
            x, tol = h.x, h.devx
            ax = abs(x)
            for k = max(h.valid+1,0):min(m, n-1)
                for j = 1:n-k
                    bjp = b[j+1] + b[j] * x
                    b[j+1] = bjp
                    e[j+1] = hypot(hypot(abs(e[j+1]), abs(b[j])*tol +ax*e[j]), abs(bjp)*ep)
                end
            end
        end
        h.valid = m
    end
    m <= n ? (h.coeff[n+1-m], h.devcoeff[n+1-m]) : (zero(T), zero(V))
end

#
# evaluate polynomial f and all derivatives at given approximate root-values z
#
function horner_analysis(z::AbstractVector{S}, f::AbstractVector{T}) where {S<:Number,T<:Number}
    thresh = 1000.0
    m = length(z)
    n = length(z)-1
    if m > length(f)-1
        throw(ArgumentError("polynomial of degree $n with $m zeros"))
    end
	h = Array{Horner}(m)
    for i = 1:m
        h[i] = Horner(z[i], f, eps(abs(z[i])) * 100, map(eps, f) * 100)
    end
	for i = 1:m
        println("root $i: $(z[i])")
        mult = 0
        search = true
		for k = 0:n
            hv, he = h[i](k)
            deviation = abs(hv / (he + eps(zero(real(S)))))
            @printf("%d %d %6.2g %8.8g %8.8g\n", i, k, deviation, abs(hv), he)
            if search && deviation < thresh
                mult = k + 1
            else
                search = false
            end
        end
        println("root $i: $(z[i]) estimated multiplicity: $mult")
    end
end


text:
# calculate varianve of sensor data

module CalcVariance
    using DataFrames, CSV, Statistics

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

        # calculate from definition
        zs = df_200_mm.lidar # observation
        mean_def = sum(zs) / length(zs) # mean
        diff_square = [(z - mean_def)^2 for z in zs]

        sampling_var = sum(diff_square) / length(zs) # sampling variance
        unbiased_var = sum(diff_square) / (length(zs) - 1) # unbiased variance

        println("Sampling Variance = $sampling_var")
        println("Unbiased Variance = $unbiased_var")

        # calculate by Statistics
        stats_sampling_var = Statistics.var(zs, corrected=false)
        stats_unbiased_var = Statistics.var(zs, corrected=true)

        println("Sampling Variance by Statistics = $stats_sampling_var")
        println("Unbiased Variance by Statistics = $stats_unbiased_var")
    end
end

max_stars_repo_name: ShisatoYano/JuliaAutonomy

max_issues_repo_name: ShisatoYano/JuliaAutonomy

max_forks_repo_name: ShisatoYano/JuliaAutonomy


text:
@testset "quadratic.jl" begin

    @testset "scaled" begin

        p(x)=1+2*x+3*x*x

        map = NLS_Fit.Map_Quadratic(1=>p(1),3=>p(3),10=>p(10))

        @test parameter_size(map) == 3

        θ = Float32[1, 1, 1]
        X_hat = Int[1,2,3,4,10]
        X = eval_map(map,X_hat,θ)

        @test eltype(X) == Float32
        @test X ≈ p.(X_hat)

    end

    @testset "unscaled" begin

        p(x) = 3+2*x+1*x*x

        map = NLS_Fit.Map_Quadratic(1,3,10)

        @test parameter_size(map) == 3

        θ = Float32[p(1), p(3), p(10)]
        X_hat = Int[1,2,3,4,10]
        X = eval_map(map,X_hat,θ)

        @test eltype(X) == Float32
        @test X ≈ p.(X_hat)

    end
end


text:
"""
    dot(x1,y1,x2,y2) -> Real

Calculate the dot product between two vectors.
"""
dot(x1::Real, y1::Real, x2::Real, y2::Real) = x1*x2 + y1*y2

"""
    dot(z1,z2) -> Real

Calculate the dot product between two complex number (like 2D vectors).
"""
dot(z1::Number, z2::Number) = real(z1)*real(z2) + imag(z1)*imag(z2)

"""
    perp(x,y) -> Tuple{Real,Real}

Return the pependicular vector to \$(x,y)\$.
"""
perp(x::Real, y::Real) = -y, x

"""
    perp(z) -> Complex

Return the pependicular complex number to \$z\$.
"""
perp(z::Number) = -imag(z) + im*real(z)

"""
"""
tripledet(x1::Real, y1::Real, x2::Real, y2::Real, x3::Real, y3::Real) =
  x1*(y2-y3)-y1*(x2-x3)+x2*y3-y2*y3

"""
"""
tripledet(z1::Number, z2::Number, z3::Number) =
  tripledet(real(z1), imag(z1), real(z2), imag(z2), real(z3), imag(z3))

#=
ToDo!
circumcenter
circumradius
=#


text:
using MultiAgentPathFinding
using Documenter

DocMeta.setdocmeta!(MultiAgentPathFinding, :DocTestSetup, :(using MultiAgentPathFinding); recursive=true)

makedocs(;
    modules=[MultiAgentPathFinding],
    authors="Guillaume Dalle <22795598+gdalle@users.noreply.github.com> and contributors",
    repo="https://github.com/gdalle/MultiAgentPathFinding.jl/blob/{commit}{path}#{line}",
    sitename="MultiAgentPathFinding.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://gdalle.github.io/MultiAgentPathFinding.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        "API reference" => "api.md",
    ],
)

deploydocs(;
    repo="github.com/gdalle/MultiAgentPathFinding.jl",
    devbranch="main",
)


text:
module Sorting_Test

using Intan, SpikeSorting, Gtk.ShortNames, Cairo, MAT, JLD

if VERSION > v"0.7-"
    using Test
else
    using Base.Test
end

myamp=RHD2164("PortA1")
d=Debug(string(dirname(Base.source_path()),"/data/qq.mat"),"qq")
myt=Task_NoTask()
myfpga=FPGA(1,myamp)
(myrhd,ss,myfpgas)=makeRHD([myfpga],debug=d)

handles = makegui(myrhd,ss,myt,myfpgas)

sleep(1.0)

#Initialization
Intan.init_cb(handles.init.handle,(handles,myrhd,myt,myfpgas))

#Run
setproperty!(handles.run,:active,true)
sleep(1.0)
Intan.run_cb(handles.run.handle,(handles,myrhd,ss,myt,myfpgas))

#Calibration
sleep(1.0)

myctx2=Gtk.getgc(handles.sc.c2)

sleep(1.0)

#Add Unit

@testset "Add Unit" begin
SpikeSorting.b2_cb_template(handles.run.handle,(handles.sc,))

press=Gtk.GdkEventButton(Gtk.GdkEventType.BUTTON_PRESS, Gtk.gdk_window(handles.sc.c2),Int8(0),UInt32(0),142.0,252.0,convert(Ptr{Float64},C_NULL),UInt32(0),UInt32(1),C_NULL,0.0,0.0)
signal_emit(handles.sc.c2,"button-press-event",Bool,press)

press=Gtk.GdkEventButton(Gtk.GdkEventType.BUTTON_PRESS, Gtk.gdk_window(handles.sc.c2),Int8(0),UInt32(0),110.0,253.0,convert(Ptr{Float64},C_NULL),UInt32(0),UInt32(1),C_NULL,0.0,0.0)
signal_emit(handles.sc.c2,"button-release-event",Bool,press)

sleep(1.0)

    @test handles.total_clus[1] == 1
    @test handles.total_clus[2] == 0

end
#Delete Unit

@testset "Delete Unit" begin
SpikeSorting.b1_cb_template(handles.run.handle,(handles.sc,))

    @test handles.total_clus[1] == 0
    @test handles.total_clus[2] == 0

end
sleep(1.0)

#Add Unit
SpikeSorting.b2_cb_template(handles.run.handle,(handles.sc,))

press=Gtk.GdkEventButton(Gtk.GdkEventType.BUTTON_PRESS, Gtk.gdk_window(handles.sc.c2),Int8(0),UInt32(0),142.0,252.0,convert(Ptr{Float64},C_NULL),UInt32(0),UInt32(1),C_NULL,0.0,0.0)
signal_emit(handles.sc.c2,"button-press-event",Bool,press)

press=Gtk.GdkEventButton(Gtk.GdkEventType.BUTTON_PRESS, Gtk.gdk_window(handles.sc.c2),Int8(0),UInt32(0),110.0,253.0,convert(Ptr{Float64},C_NULL),UInt32(0),UInt32(1),C_NULL,0.0,0.0)
signal_emit(handles.sc.c2,"button-release-event",Bool,press)

sleep(1.0)

#Add Second Unit
SpikeSorting.b2_cb_template(handles.run.handle,(handles.sc,))

press=Gtk.GdkEventButton(Gtk.GdkEventType.BUTTON_PRESS, Gtk.gdk_window(handles.sc.c2),Int8(0),UInt32(0),186.0,287.0,convert(Ptr{Float64},C_NULL),UInt32(0),UInt32(1),C_NULL,0.0,0.0)
signal_emit(handles.sc.c2,"button-press-event",Bool,press)

press=Gtk.GdkEventButton(Gtk.GdkEventType.BUTTON_PRESS, Gtk.gdk_window(handles.sc.c2),Int8(0),UInt32(0),156.0,293.0,convert(Ptr{Float64},C_NULL),UInt32(0),UInt32(1),C_NULL,0.0,0.0)
signal_emit(handles.sc.c2,"button-release-event",Bool,press)

sleep(1.0)

    @test handles.total_clus[1] == 2
    @test handles.total_clus[2] == 0

#Change second unit

press=Gtk.GdkEventButton(Gtk.GdkEventType.BUTTON_PRESS, Gtk.gdk_window(handles.sc.c2),Int8(0),UInt32(0),186.0,287.0,convert(Ptr{Float64},C_NULL),UInt32(0),UInt32(1),C_NULL,0.0,0.0)
signal_emit(handles.sc.c2,"button-press-event",Bool,press)

press=Gtk.GdkEventButton(Gtk.GdkEventType.BUTTON_PRESS, Gtk.gdk_window(handles.sc.c2),Int8(0),UInt32(0),156.0,293.0,convert(Ptr{Float64},C_NULL),UInt32(0),UInt32(1),C_NULL,0.0,0.0)
signal_emit(handles.sc.c2,"button-release-event",Bool,press)

sleep(1.0)

    @test handles.total_clus[1] == 2
    @test handles.total_clus[2] == 0

#Show Templates
SpikeSorting.b4_cb_template(handles.run.handle,(handles.sc,))

sleep(1.0)

#Set draw template checkbox
setproperty!(handles.sort_widgets.check,:active,true)

sleep(1.0)

#Collect Templates
press=Gtk.GdkEventButton(Gtk.GdkEventType.BUTTON_PRESS, Gtk.gdk_window(handles.sort_widgets.b3),Int8(0),UInt32(0),0.0,0.0,convert(Ptr{Float64},C_NULL),UInt32(0),UInt32(1),C_NULL,0.0,0.0)
signal_emit(handles.sort_widgets.b3,"clicked",Bool,press)

sleep(1.0)

press=Gtk.GdkEventButton(Gtk.GdkEventType.BUTTON_PRESS, Gtk.gdk_window(handles.sort_widgets.b3),Int8(0),UInt32(0),0.0,0.0,convert(Ptr{Float64},C_NULL),UInt32(0),UInt32(1),C_NULL,0.0,0.0)
signal_emit(handles.sort_widgets.b3,"clicked",Bool,press)

sleep(1.0)

press=Gtk.GdkEventButton(Gtk.GdkEventType.BUTTON_PRESS, Gtk.gdk_window(handles.sc.c2),Int8(0),UInt32(0),156.0,293.0,convert(Ptr{Float64},C_NULL),UInt32(0),UInt32(1),C_NULL,0.0,0.0)
signal_emit(handles.sc.c2,"button-release-event",Bool,press)

sleep(1.0)

press=Gtk.GdkEventButton(Gtk.GdkEventType.BUTTON_PRESS, Gtk.gdk_window(handles.sc.pause_button),Int8(0),UInt32(0),0.0,0.0,convert(Ptr{Float64},C_NULL),UInt32(0),UInt32(1),C_NULL,0.0,0.0)
signal_emit(handles.sc.pause_button,"clicked",Bool,press)

sleep(1.0)

    @test handles.total_clus[1] == 2
    @test handles.total_clus[2] == 0

end


text:
"""
    MedianVarianceBinning([minsize::Int = 10, maxbins::Int = typemax(Int)])

Dynamic binning scheme of the probability simplex with at most `maxbins` bins that each
contain at least `minsize` samples.

The data set is split recursively as long as it is possible to split the bins while
satisfying these conditions. In each step, the bin with the maximum variance of predicted
probabilities for any component is selected and split at the median of the predicted
probability of the component with the largest variance.
"""
struct MedianVarianceBinning <: AbstractBinningAlgorithm
    minsize::Int
    maxbins::Int

    function MedianVarianceBinning(minsize, maxbins)
        minsize ≥ 1 || error("minimum number of samples must be positive")
        maxbins ≥ 1 || error("maximum number of bins must be positive")

        return new(minsize, maxbins)
    end
end

MedianVarianceBinning(minsize::Int=10) = MedianVarianceBinning(minsize, typemax(Int))

function perform(
    alg::MedianVarianceBinning,
    predictions::AbstractVector{<:AbstractVector{<:Real}},
    targets::AbstractVector{<:Integer},
)
    @unpack minsize, maxbins = alg

    # check if binning is not possible
    nsamples = length(predictions)
    nsamples < minsize && error("at least $minsize samples are required")

    # check if only trivial binning is possible
    minsplit = 2 * minsize
    (nsamples < minsplit || maxbins == 1) && return [Bin(predictions, targets)]

    # find dimension with maximum variance
    idxs_predictions = collect(1:nsamples)
    GC.@preserve idxs_predictions begin
        max_var_predictions, argmax_var_predictions = max_argmax_var(
            predictions, idxs_predictions
        )

        # create priority queue and empty set of bins
        queue = PriorityQueue(
            (idxs_predictions, argmax_var_predictions) => max_var_predictions,
            Base.Order.Reverse,
        )
        bins = Vector{typeof(Bin(predictions, targets))}(undef, 0)

        nbins = 1
        while nbins < maxbins && !isempty(queue)
            # pick the set with the largest variance
            idxs, argmax_var = dequeue!(queue)

            # compute indices of the two subsets when splitting at the median
            idxsbelow, idxsabove = unsafe_median_split!(idxs, predictions, argmax_var)

            # add a bin of all indices if one of the subsets is too small
            # can happen if there are many samples that are equal to the median
            if length(idxsbelow) < minsize || length(idxsabove) < minsize
                push!(bins, Bin(predictions[idxs], targets[idxs]))
                continue
            end

            for newidxs in (idxsbelow, idxsabove)
                if length(newidxs) < minsplit
                    # add a new bin if the subset can not be split further
                    push!(bins, Bin(predictions[newidxs], targets[newidxs]))
                else
                    # otherwise update the queue with the new subsets
                    max_var_newidxs, argmax_var_newidxs = max_argmax_var(
                        predictions, newidxs
                    )
                    enqueue!(queue, (newidxs, argmax_var_newidxs), max_var_newidxs)
                end
            end

            # in total one additional bin was created
            nbins += 1
        end

        # add remaining bins
        while !isempty(queue)
            # pop queue
            idxs, _ = dequeue!(queue)

            # create bin
            push!(bins, Bin(predictions[idxs], targets[idxs]))
        end
    end

    return bins
end

function max_argmax_var(x::AbstractVector{<:AbstractVector{<:Real}}, idxs)
    # compute variance along the first dimension
    maxvar = unsafe_variance_welford(x, idxs, 1)

    maxdim = 1
    for d in 2:length(x[1])
        # compute variance along the d-th dimension
        vard = unsafe_variance_welford(x, idxs, d)

        # update current optimum if required
        if vard > maxvar
            maxvar = vard
            maxdim = d
        end
    end

    return maxvar, maxdim
end

# use Welford algorithm to compute the unbiased sample variance
# taken from: https://github.com/JuliaLang/Statistics.jl/blob/da6057baf849cbc803b952ef7adf979ae3a9f9d2/src/Statistics.jl#L184-L199
# this function is unsafe since it does not perform any bounds checking
function unsafe_variance_welford(
    x::AbstractVector{<:AbstractVector{<:Real}}, idxs::Vector{Int}, dim::Int
)
    n = length(idxs)

    @inbounds begin
        M = x[idxs[1]][dim] / 1
        S = zero(M)
        for i in 2:n
            value = x[idxs[i]][dim]

            new_M = M + (value - M) / i
            S += (value - M) * (value - new_M)
            M = new_M
        end
    end

    return S / (n - 1)
end

# this function is unsafe since it leads to undefined behaviour if the
# outputs are accessed afer `idxs` has been garbage collected
function unsafe_median_split!(
    idxs::Vector{Int}, x::AbstractVector{<:AbstractVector{<:Real}}, dim::Int
)
    n = length(idxs)

    if length(idxs) < 2
        cutoff = 0
    else
        # partially sort the indices `idxs` according to the corresponding values in the
        # `d`th component of`x`
        m = div(n, 2) + 1
        f = let x = x, dim = dim
            idx -> x[idx][dim]
        end
        partialsort!(idxs, 1:m; by=f)

        # figure out all values < median
        # the median is `x[idxs[m]][dim]`` for vectors of odd length
        # and `(x[idxs[m - 1]][dim] + x[idxs[m]][dim]) / 2` for vectors of even length
        if x[idxs[m - 1]][dim] < x[idxs[m]][dim]
            cutoff = m - 1
        else
            # otherwise obtain the last value < median
            firstidxs = unsafe_wrap(Array, pointer(idxs, 1), m - 1)
            cutoff = searchsortedfirst(firstidxs, idxs[m]; by=f) - 1
        end
    end

    # create two new arrays that refer to the two subsets of indices
    idxsbelow = unsafe_wrap(Array, pointer(idxs, 1), cutoff)
    idxsabove = unsafe_wrap(Array, pointer(idxs, cutoff + 1), n - cutoff)

    return idxsbelow, idxsabove
end


text:
# ---
# title: 438. Find All Anagrams in a String
# id: problem438
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Hash Table
# link: <https://leetcode.com/problems/find-all-anagrams-in-a-string/description/>
# hidden: true
# ---
#
# Given a string **s** and a **non-empty** string **p** , find all the start
# indices of **p** 's anagrams in **s**.
#
# Strings consists of lowercase English letters only and the length of both
# strings **s** and **p** will not be larger than 20,100.
#
# The order of output does not matter.
#
# **Example 1:**
#
#
#
#     Input:
#     s: "cbaebabacd" p: "abc"
#
#     Output:
#     [0, 6]
#
#     Explanation:
#     The substring with start index = 0 is "cba", which is an anagram of "abc".
#     The substring with start index = 6 is "bac", which is an anagram of "abc".
#
#
# **Example 2:**
#
#
#
#     Input:
#     s: "abab" p: "ab"
#
#     Output:
#     [0, 1, 2]
#
#     Explanation:
#     The substring with start index = 0 is "ab", which is an anagram of "ab".
#     The substring with start index = 1 is "ba", which is an anagram of "ab".
#     The substring with start index = 2 is "ab", which is an anagram of "ab".
#
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
@testset "Array" begin
    @test size(Leaf(Tape(), ones(1, 2, 3, 4))) == (1, 2, 3, 4)
    @test length(Leaf(Tape(), 1:3)) == 3

    let rng = MersenneTwister(123456)
        x = randn(2, 10)
        f1 = x̂ -> reshape(x̂, 5, 4)
        f2 = x̂ -> reshape(x̂, (5, 4))
        @test check_errs(f1, f1(x), x, randn(size(x)...))
        @test check_errs(f2, f2(x), x, randn(size(x)...))
    end

    a = rand(3, 2); b = rand(3); c = rand(3, 3);
    f(a, b, c) = sum(hcat(2*a, 3*b, 4*c))
    @test ∇(f)(a,b,c) == (2*ones(3, 2), 3*ones(3), 4*ones(3, 3))
    a = rand(2, 4); b = rand(1, 4); c = rand(3, 4);
    f(a, b, c) = sum(vcat(2*a, 3*b, 4*c))
    @test ∇(f)(a,b,c) == (2*ones(2, 4), 3*ones(1, 4), 4*ones(3, 4))

    @test check_errs(x->fill(x, 4, 4), randn(4, 4), randn(), randn())
    @test check_errs(x->fill(x, (4, 4)), randn(4, 4), randn(), randn())
    x = Leaf(Tape(), 6)
    y = fill(x, 3)
    @test y isa Branch{Vector{Int}}
    @test getfield(y, :f) === Base.fill
end


text:
using LinearAlgebra
using StatsBase
using Plots

# m = minlp(false)
# optimize!(m)
# mcost = JuMP.getobjectivevalue(m)
# msol = solution(m)

gm = test_gqp()
set_optimizer(gm, SOLVER_SILENT)

# bbcs = [bbl for bbl in gm.bbls if bbl isa BlackBoxClassifier]
bbr = [bbl for bbl in gm.bbls if bbl isa BlackBoxRegressor][1]
uniform_sample_and_eval!(gm)
scatter(bbr.X[:,"x[1]"], bbr.X[:, "x[2]"], bbr.Y)

# # Learning BBCs just for fun.
# for bbc in bbcs
#     learn_constraint!(bbc)
#     update_tree_constraints!(gm, bbc)
# end

# # Different types of regression learners for testing purposes.
learn_constraint!(bbr, regression_sparsity=0, max_depth = 3, minbucket = 2*length(bbr.vars) + 2)
update_tree_constraints!(gm, bbr)
show_trees(bbr)
optimize!(gm)

# Plotting all of the different leaves
scatter()
idx = length(bbr.learners)
ul_data = bbr.ul_data[idx]
lnr = bbr.learners[idx]
leaf_binnings = IAI.apply(lnr, bbr.X)
for leaf in find_leaves(lnr)
    in_leaf = findall(x -> x == leaf, leaf_binnings)
    Plots.display(scatter!(bbr.X[in_leaf,"x[1]"], bbr.X[in_leaf, "x[2]"], bbr.Y[in_leaf], label = leaf))
end

# # Scoring different trees
# for tree in bbr.learners
#     println("Tree score: ", IAI.score(tree, bbr.X, bbr.Y))
# end
# DETECTING CONVEXITY

all_leaves = find_leaves(lnr)
points_in_leaves = Dict(leaf => [] for leaf in all_leaves)
for i=1:length(leaf_binnings)
    push!(points_in_leaves[leaf_binnings[i]], i)
end
centroids = Dict(leaf => [mean(col) for col in eachcol(
             bbr.X[points_in_leaves[leaf], :])] for leaf in all_leaves)

convex_bins = Any[[leaf] for leaf in all_leaves]
parents = unique([IAI.get_parent(lnr, leaf[1]) for leaf in convex_bins])

for par in parents
    lower_child = IAI.get_lower_child(lnr, par)
    upper_child = IAI.get_upper_child(lnr, par)
    if IAI.is_leaf(lnr, lower_child) && IAI.is_leaf(lnr, upper_child)
        α0, α = ul_data[lower_child]
        β0, β = ul_data[upper_child]
        lower_X = bbr.X[points_in_leaves[lower_child], :]
        lower_Y = bbr.Y[points_in_leaves[lower_child]]
        lower_predictions = Matrix(lower_X)*α .+ α0
        upper_X = bbr.X[points_in_leaves[upper_child], :]
        upper_Y = bbr.Y[points_in_leaves[upper_child]]
        upper_predictions = Matrix(upper_X)*β .+ β0

        l_in_l = lower_Y - Matrix(lower_X)*α .- α0
        u_in_u = upper_Y - Matrix(upper_X)*β .- β0
        l_in_u = lower_Y - Matrix(lower_X)*β .- β0
        u_in_l = upper_Y - Matrix(upper_X)*α .- α0
        if all(l_in_u .>= 0) && all(u_in_l .>= 0)
            push!(convex_bins, [par, [lower_child, upper_child]])

        else
            println("Violations ", (sum(l_in_u .>= 0) + sum(u_in_l .>= 0))/(length(l_in_u) + length(u_in_l)))
        end
    end
end


text:
module JLLWrappers

# We need to glue expressions together a lot
function excat(exs::Union{Expr,Nothing}...)
    ex = Expr(:block)
    for exn in exs
        if exn !== nothing
            append!(ex.args, exn.args)
        end
    end
    return esc(ex)
end

include("toplevel_generators.jl")
include("wrapper_generators.jl")
include("runtime.jl")

end # module


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

JLLWrappers.@generate_main_file_header("Jansson")
JLLWrappers.@generate_main_file("Jansson", UUID("83cbd138-b029-500a-bd82-26ec0fbaa0df"))
end  # module Jansson_jll


text:
###############################################################################
#
# GEinvestment.jl solves the KT general equlibrium model of firm investment with labor
#
# Shane McMiken
# November 2021
# Boston College
#
# CONSUMER PROBLEM:
#   max_{N_{t+s}, B_{t+s+1}} \sum_{s ≥ 0} β^s (log C_{t+s} - ϕ N_{t+s})
#   s.t. C_{t+s} + B_{t+s+1} = W * N_{t+s} + R * B_{t+s} + ∫ d(z,k) dF(z,k)
#
#   steady state: R = β^-1, W = ϕ * C
#
# RECURSIVE FIRM PROBLEM:
# V(Z,K) = max_K' {R(Z,K,K') + (1/(1+r))* E[V(Z',K')|Z] }
#
# WHERE R(z,k,k') = Z * K^α * N^ν - pK + p(1-δ)K - AC(Z,K,K')
#
################################################################################

#--------------------------------#
#         PREAMBLE
#--------------------------------#
using Plots # plots
using Distributions # cdf(normal())
using Parameters  # macros: @with_kw & @unpack
#using DelimitedFiles # writedlm

#--------------------------------#
#       TAUCHEN (1986)
#--------------------------------#

function MyTauchen(σ_a::Float64,ρ_a::Float64,na::Int64,m::Float64)

    # Discrete markov process approximation for AR(1)
    # y = ρ_a * y_  + σ_a * ϵ
    # m: number of standard deviations
    # σ_a: standard deviations
    # ρ_a: AR(1) parameter
    # na: number of grid points

    vecsize = na
    σ_y = σ_a / sqrt(1-(ρ_a^2)) # Stationary distribution variance
    step = 2 * m * σ_y / (vecsize-1)
    agrid  = [-m*σ_y + (i-1)*step for i∈1:na]

    # Calculate transition probabilities, AR(1) - Tauchen
    P = zeros(na, na)

    mm = agrid[2] - agrid[1]; # re-compute step
    for j = 1:na
      for k = 1:na
        if(k == 1) # first column
          P[j, k] = cdf(Normal(), (agrid[k] - ρ_a*agrid[j] + (mm/2))/σ_a)
        elseif(k == na) # second-last column
          P[j, k] = 1 - cdf(Normal(), (agrid[k] - ρ_a*agrid[j] - (mm/2))/σ_a)
        else
          P[j, k] = cdf(Normal(), (agrid[k] - ρ_a*agrid[j] + (mm/2))/σ_a) - cdf(Normal(), (agrid[k] - ρ_a*agrid[j] - (mm/2))/σ_a)
        end
      end #k
    end #j

    # agrid = exp.(agrid) # Exponentiate grid

    return agrid, P #return grid and Markov transition matrix
end

#--------------------------------#
#         PARAMETERIZATION
#--------------------------------#

@with_kw struct params

    # firm parameters
    α::Float64      = 0.25 # capital share
    δ::Float64      = 0.1 # depreciation rate
    r::Float64      = 0.04 # interest rate
    ρ::Float64      = 0.8 # technology persistence
    σ::Float64      = 0.2 # technology volatility
    γ_c::Float64    = 2.5 # variable investment adjustment cost
    γ_f::Float64    = 0.0 # fixed investment adjustment cost

    # household parameters
    β::Float64      = 1.0/0.04 # discount rate
    ϕ::Float64      = 2.0 #disutility of labor
    ν::Float64      = 0.5 #labor share

    # wage bisection bounds
    a::Float64      = 0.4 #wage lower bound
    c::Float64      = 1.8 #wage upper bound

end

#--------------------------------#
# INNER LOOP - VFI & ERGODIC DIST
#--------------------------------#

function vfivec(W::Float64, kkgr::Array, zzgr::Array, kgrid::Array, zgrid::Array;α = α, ν = ν, r = r, δ = δ, γ_f = γ_f, γ_c = γ_c, maxiter = 1000, tol = 1e-4, nfix = 10, iter = 0, verr = Inf)

    # objects
    nz          = length(zgrid)
    nk          = length(kgrid)
    labor       = zzgr.^(1.0/(1.0-ν)).*(ν/W).^(1.0/(1.0-ν)).*kkgr.^(α/(1.0-ν))
    output      = zzgr .* kkgr.^α .* labor.^ν
    investment  = kgrid' .- (1.0 - δ).*kkgr # kgrid - k' policy
    adjcost     = γ_c/2.0 .* (investment./kkgr).^2 .*kkgr + γ_f.*output.*(kgrid'.>(1-δ)*kkgr)
    dividend    = repeat(output,1,nk) - W * repeat(labor,1,nk) - investment - adjcost

    # steady state
    kbar        = (α/(r+δ))^((1.0-ν)/(1.0-α-ν))*(ν/W)^(ν/(1.0-α-ν))
    nbar        = (ν/W)^(1.0/(1.0-ν))*kbar^(α/(1.0-ν))
    Vbar        = ((1.0+r)/r)*(kbar^α*nbar^ν - δ*kbar -W*nbar)

    # initialise
    V_init  = fill(Vbar, (nz,nk))
    V_next  = Array{Float64,1}
    index   = Array{CartesianIndex{2},2}

    while (verr > tol && iter < maxiter) # inner-loop value function

        if mod(iter,nfix) == 0 # Howard acceleration

            # obj dims: V(z×k,k')
            obj     = dividend + (1.0/(1.0+r)).*repeat(P*V_init,nk) # compute values for V(z×k,k')
            V_next, index = findmax(obj, dims = 2) # get maximum

        else

            # Re-use conjectured index
            idx = [index[i][2] for i∈1:(nz*nk)]

            # E(V) condition on A(t) and each possible K(t+1)
            EVp = P*V_init
            evp_k = zeros(nk*nz,)

            for ind = 1:(nz*nk)
                zz = convert(Int, floor(mod(ind-0.05,nz))+1)
                evp_k[ind] = EVp[zz,idx[ind]]
            end

            # Compute return function using conjectured policy - idx
            investment_idx  = kgrid[idx] .- (1.0 - δ).*kkgr
            adjcost_idx     = γ_c/2.0 .* ((kgrid[idx] .- (1.0-δ)*kkgr)./kkgr).^2 .*kkgr + γ_f.*output.*(kgrid[idx].>(1-δ)*kkgr)
            dividend_idx    = output - W * labor - investment_idx - adjcost_idx

            # Re-compute value function under conjectured optimal policy
            V_next = dividend_idx + (1.0+r).^(-1) .* evp_k

        end

        # check tolerance
        V_next  = reshape(V_next,nz,nk)
        verr    = maximum(abs.(V_next .- V_init))
        V_init  = copy(V_next)
        iter   += 1
        #print("iteration: ",iter, ", critical val: ", verr*1.0e4," times 10^4", "\n")

    end # inner-loop: value function converges

    index = [index[i][2] for i∈1:(nz*nk)] # policy funciton
    return index, labor, output, investment, adjcost

end

function young(index::Vector{Int64}; maxiter = 1000, tol = 1e-4, iter = 0, crit = Inf, nz = nz, nk = nk)

    # initialise
    μ_init  = fill((1/(nz*nk)), (nz,nk))

    while (crit > tol && iter < maxiter)

        # intialise
        μ_next = zeros(nz,nk)

        # use policies from vfi to infer masses next period
        for ind = 1:(nz*nk)
        kk = convert(Int, floor((ind-0.05)/nz))+1
        zz = convert(Int, floor(mod(ind-0.05,nz))+1)
        μ_next[zz,index[ind]] =  μ_next[zz,index[ind]] + P[:,zz]'*μ_init[:,kk]
        end

        # check tolerance
        crit = maximum(abs.(μ_init .- μ_next))
        μ_init = μ_next
        iter += 1

        #print("iteration: ",iter, ", critical val: ", crit*1.0e5," times 10^5", "\n")

    end
    μ_ergodic = vec(μ_init)
    return μ_ergodic
end

#--------------------------------#
# OUTER LOOP - WAGE
#--------------------------------#

function geinvest(a_init::Float64, c_init::Float64, kgrid::Array,zgrid::Array; wmaxiter::Int64 = 50, tol::Float64 = 1e-4, nk::Int64 = nk, nz::Int64 = nz, α::Float64=α, δ::Float64 = δ, r::Float64 = r, γ_c::Float64 = γ_c, γ_f::Float64 = γ_f, ϕ::Float64 = ϕ, ν::Float64 = ν, iter = 0, err = Inf)

    # wage initialisation
    a = copy(a_init)
    c = copy(c_init)
    b = (a+c)/2.0
    W_next = (a_init + c_init)/2.0

    # useful grids
    kkgr = vec(repeat(kgrid',nz))
    zzgr = repeat(zgrid,nk,1)
    μ_ergodic = Vector{Float64}
    index = Vector{Int64}

    while (err > tol && iter < wmaxiter) #outer loop for wage

        # update prices
        W = copy(W_next)

        # value-function iteration
        index_stat, labor, output, investment, adjcost = vfivec(W,kkgr,zzgr,kgrid,zgrid)
        index = copy(index_stat)

        # obtain ergodic distribution
        μ = young(index)
        μ_ergodic = copy(μ)

        #------------------------------------------#
        #  COMPUTE AGGREGATES
        #------------------------------------------#

        Lagg        = sum(μ_ergodic.*labor)
        Yagg        = sum(μ_ergodic.*output)
        Kagg        = sum(μ_ergodic.*kkgr)
        Zagg        = sum(μ_ergodic.*zzgr)
        Iagg        = sum(μ_ergodic.*investment[index])
        ACagg       = sum(μ_ergodic.*adjcost[index])
        Cagg        = Yagg .- Iagg .- ACagg

        # compute error for prices outer-loop
        diff        = W - ϕ * Cagg
        err         = abs(diff)
        iter        += 1

        # update prices for outer loop
        if diff > 0
            c = copy(b)
        elseif diff < 0
            a = copy(b)
        end
        b = (a+c)/2.0
        W_next = b

        #print("iteration: ",iter, ", critical val: ", err*1.0e4," times 10^4 ", " wage ", W, " diff ", diff,   " Cagg ", Cagg,  "\n")

    end # outer-loop: wage converges

    return index,W_next,μ_ergodic
end

#--------------------------------#
#         RUN MODEL
#--------------------------------#

@unpack α, δ, r, ρ, σ, γ_c, γ_f, β, ϕ, ν, a, c= params()

# technology grid
cover = 5
nz    = 2*cover+1
m     = 2.5
zgrid, P = MyTauchen(σ,ρ,nz,m)
zgrid = exp.(zgrid) # exponentiate grid

# capital grid
nk = 100
kmax = (α/(r+δ))^((1.0-ν)/(1.0-α-ν))*(ν/a)^(ν/(1.0-α-ν))*maximum(zgrid)^(1.0/(1.0-α-ν))
kgrid = [kmax*(1-δ)^(nk-i-1) for i∈0:(nk-1)]

# Convex adjustment cost
@time index, W, μ_ergodic  = geinvest(a, c, kgrid, zgrid)

# policy function

#ergodicdist = heatmap(kgrid,zgrid,reshape(μ_ergodic,nz,nk), title = "Ergodic distribution")
#display(ergodicdist)
#savefig(ergodicdist,"./src/GEergodicdistribution.png")

policyfunc = heatmap(kgrid,zgrid,reshape(index,(nz,nk)), title = "k'(z,k) Policy")
display(policyfunc)
savefig(policyfunc,"./src/GEpolicyfunction_ConvexAC.png")

# Fixed Adjustment cost
γ_c, γ_f = 0.0, 0.05
@time index, W, μ_ergodic  = geinvest(a, c, kgrid, zgrid)

# policy function
policyfunc = heatmap(kgrid,zgrid,reshape(index,(nz,nk)), title = "k'(z,k) Policy")
display(policyfunc)
savefig(policyfunc,"./src/GEpolicyfunction_FixedAC.png")

max_stars_repo_name: shanemcmiken/PEinvestment

max_issues_repo_name: shanemcmiken/PEinvestment

max_forks_repo_name: shanemcmiken/PEinvestment


text:
using Revise
using JDF, CSV, DataFrames
@time a = gf()

# 2G file
# @time a = CSV.read(
#     "C:/data/Performance_All/Performance_2010Q3.txt",
#     delim = '|',
#     header = false
# );

b = type_compress(Array(a[!, :ORIGIN]));

io = iow()
@time metadata = compress_then_write(b, io)
close(io)

using Revise
using JDF, CSV, DataFrames
# using JLSO
# metadata = JLSO.load("C:/data/metatmp")["data"]

io = ior()
@time oo = column_loader(metadata.type, io, metadata);
close(io)

all(b.==oo)


text:
# methods common to IBD and BD

ninv(pa, ::S) = length(pa.sx.x)
ninv(pa, ::G) = length(pa.gx.x)

function add_gprecon!(pa, gprecon)

	copyto!(pa.gx.precon, gprecon)
	copyto!(pa.gx.preconI, pa.gx.precon)
	for i in eachindex(gprecon)
		if(!(iszero(gprecon[i])))
			pa.gx.preconI[i]=inv(pa.gx.precon[i])
		end
	end
end
function add_sprecon!(pa, sprecon)
	copyto!(pa.sx.precon, sprecon)
	copyto!(pa.sx.preconI, pa.sx.precon)
	for i in eachindex(sprecon)
		if(!(iszero(sprecon[i])))
			pa.sx.preconI[i]=inv(pa.sx.precon[i])
		end
	end
end
function add_gweights!(pa, gweights)
	(gweights===nothing) && (gweights=ones(pa.optm.cal.g))
	copyto!(pa.gx.weights, gweights)
end

# Remove preconditioners from pa
function remove_gprecon!(pa; including_zeros=false)
	for i in eachindex(pa.gx.precon)
		if((pa.gx.precon[i]≠0.0) || including_zeros)
			pa.gx.precon[i]=1.0
			pa.gx.preconI[i]=1.0
		end
	end
end

"""
Remove weights from pa
"""
function remove_gweights!(pa; including_zeros=false)
	for i in eachindex(pa.gx.weights)
		if((pa.gx.weights[i]≠0.0) || including_zeros)
			pa.gx.weights[i]=1.0
		end
	end
end

# these are not use at this moment

#=
function update_window!(paf::FourierConstraints, obs)
	pac=obs
	paf=paf
	Conv.pad!(pac.d, pac.dpad, pac.dlags[1], pac.dlags[2], pac.np2)
	mul!(pac.dfreq, pac.dfftp, pac.dpad)
	for i in eachindex(paf.window)
		paf.window[i]=complex(0.0,0.0)
	end

	nr=size(pac.dfreq,2)
	# stack the spectrum of data
	for j in 1:nr
		for i in eachindex(paf.window)
			paf.window[i] += (abs2(pac.dfreq[i,j]))
		end
	end

	normalize!(paf.window)
	paf.window = 10. * log10.(paf.window)

	# mute frequencies less than -40 dB
	for i in eachindex(paf.window)
		if(paf.window[i] ≤ -40)
			paf.window[i]=0.0
		else
			paf.window[i]=1.0
		end
	end
end

function apply_window_s!(s, pac, paf)
	Conv.pad!(s, pac.spad, pac.slags[1], pac.slags[2], pac.np2)
	mul!(pac.sfreq, pac.sfftp, pac.spad)
	for i in eachindex(pac.sfreq)
		pac.sfreq[i] *= paf.window[i]
	end
	mul!(pac.spad, pac.sifftp, pac.sfreq)
	Conv.truncate!(s, pac.spad, pac.slags[1], pac.slags[2], pac.np2)
end

function apply_window_g!(g, pac, paf)
	Conv.pad!(g, pac.gpad, pac.glags[1], pac.glags[2], pac.np2)
	mul!(pac.gfreq, pac.gfftp, pac.gpad)
	nr=size(pac.gfreq,2)
	for ir in 1:nr
		for i in size(pac.gfreq,1)
			pac.gfreq[i,ir] *= paf.window[i]
		end
	end
	mul!(pac.gpad, pac.gifftp, pac.gfreq)
	Conv.truncate!(g, pac.gpad, pac.glags[1], pac.glags[2], pac.np2)

end
=#


text:
# Create symbolic variables
@polyvar x y
p = 2*x^4 + 2*x^3*y - x^2*y^2 + 5*y^4

# We want constraint `p` to be a sum of squares
@constraint model p >= 0

optimize!(model)
# Solution status is `OPTIMAL` which means `p` is a sum of squares
@show termination_status(model)

max_stars_repo_name: sums-of-squares/tools

max_issues_repo_name: sums-of-squares/tools

max_forks_repo_name: sums-of-squares/tools


text:
using NeuralOcenanigans
using Test

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

max_stars_repo_name: drozzy/NeuralOceananigans

max_issues_repo_name: drozzy/NeuralOceananigans

max_forks_repo_name: drozzy/NeuralOceananigans


text:
export UndefToken, tisdefined, @tisdefined, ttypeof, tisa

"""
    UndefToken

A type that shows if the an evaluated `Token` is defined.
"""
struct UndefToken end

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

"""
    @tisdefined(t)
    tisdefined(modul, t)

Check if a `Token` that represent a `Symbol` is defined. Throws error if the `Token` is representing an `Expr`.

# Examples
```jldoctest
julia> t = collect(tokenize("Int64"))

julia> tisdefined(t)
true
```
"""
function tisdefined(modul, t::T) where {T <: Union{Token, Array{Token}, String}}
    pt = Meta.parse(t)
    return isdefined(modul,pt)
end

tisdefined(modul::Module, t::T) where {T <: Union{Symbol, Expr}} = isdefined(modul, t)

macro tisdefined(t)
    m = __module__
    return :(tisdefined($m, $t))
end

################################################################
"""
    ttypeof(t)
    ttypeof(t, check_isdefined::Bool = false)

Returns the type of an evaluated Token

If you set `check_isdefined` to `true`, and `t` is not defined in the scope it returns `UndefToken` instead of throwing an error.

# Examples
```jldoctest
julia> t = collect(tokenize("Int64"))

julia> ttypeof(t)
DataType
```
"""
function ttypeof(modul::Module, t::T, check_isdefined::Bool = false) where {T <: Union{Token, Array{Token}, String, Symbol, Expr}}
    tevaled = @teval(t, check_isdefined)
    typeof(tevaled)
    return
 end

"""
    tisa(t, Tspecified)
    tisa(t, Tspecified, check_isdefined::Bool = false)

Compares the specified type with the type of an evaluated Token

If you set `check_isdefined` to `true`, and `t` is not defined in the scope it returns `UndefToken` instead of throwing an error.

# Examples
```jldoctest
julia> t = collect(tokenize("Int64"))

julia> tisa(t, DataType)
true
```
"""
function tisa(t::T, Tspecified::Type, check_isdefined::Bool = false) where {T <: Union{Token, Array{Token}, String, Symbol, Expr}}
    tevaled = @teval(t, check_isdefined)
    isa(tevaled, Tspecified)
end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

version = v"1.8.5" # also change in raw script string

# Collection of sources required to build LibSpatialIndexBuilder
sources = [
    "http://download.osgeo.org/libspatialindex/spatialindex-src-$version.tar.bz2" =>
    "31ec0a9305c3bd6b4ad60a5261cba5402366dd7d1969a8846099717778e9a50a",
    "./patches"
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir

cd spatialindex-src-1.8.5/

patch < ${WORKSPACE}/srcdir/makefile.patch
rm Makefile.am.orig

if [ $target = "x86_64-w64-mingw32" ] || [ $target = "i686-w64-mingw32" ]; then
  patch < ${WORKSPACE}/srcdir/header-check.patch
fi

aclocal
autoconf
automake --add-missing --foreign

# Show options in the log
./configure --help

./configure --prefix=$prefix --host=$target
make
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()

# The products that we will ensure are always built
products(prefix) = [
    LibraryProduct(prefix, "libspatialindex_c", :libspatialindex_c),
    LibraryProduct(prefix, "libspatialindex", :libspatialindex)
]

# Dependencies that must be installed before this package can be built
dependencies = []

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, "SpatialIndex", version, sources, script, platforms, products, dependencies)


text:
using CSTParser
using Base.Test

import CSTParser: parse, remlineinfo!, check_base, span
# write your own tests here
include("parser.jl")
check_base()


text:
#####
##### Dependencies
#####

include("./kernel.jl")
include("./transform.jl")
include("./incremental.jl")
include("./system_data.jl")
include("./kernel_data.jl")

import LinearAlgebra: logabsdet

#####
##### Structs
#####

@doc raw"""
"""
struct KernelSystem{K<:AbstractCorrelation}
    sd::SystemData
    kd::KernelData{K}
end

#####
##### Construction
#####

#####
##### Functions
#####

function reduced_solve!(cx, ks, yu, yx)
    cx .= yx
    mul!(cx, ks.sd.W', yu, -1, 1)
    ldiv!(ks.kd.rKx, cx)
    return cx
end
function solve!(c, ks, y)
    p = size(ks.sd.W, 1)
    @views cu, cx = c[1:p], c[p+1:end]
    @views yu, yx = y[1:p], y[p+1:end]
    reduced_solve!(cx, ks, yu, yx)
    mul!(cu, ks.sd.W, cx, -1, 0)
    return cu, cx
end
function fit_tail!(d, ks, yu, cu, cx)
    d .= yu
    mul!(d, ks.kd.Ku, cu, -1, 1)
    mul!(d, ks.kd.Kux, cx, -1, 1)
    ldiv!(ks.sd.Fu', d)
    return d
end

#####
##### Prediction & Statistics
#####

function det(ks::KernelSystem)
    d, s = logabsdet(ks)
    return s * exp(d)
end
function logdet(ks::KernelSystem)
    d, s = logabsdet(ks)
    s >= 0 || throw(DomainError(s))
    return d
end
function logabsdet(ks::KernelSystem)
    s = (-1.0) ^ (size(ks.sd.U, 1) % 2)
    return 2 * logabsdet(ks.sd.Fu)[1] + logdet(ks.kd.rKx), s
end

function predict!(m, ks::KernelSystem, cu, cx, d, X′, Fx′) # this is m
    Kux′ = similar(Fx′)
    Kxx′ = Array{Float64}(undef, size(ks.sd.X, 2), size(X′, 2))

    cross_correlation!(Kux′, ks.kd.kern, ks.sd.U, X′)
    cross_correlation!(Kxx′, ks.kd.kern, ks.sd.X, X′)

    mul!(m, Kux′', cu, 1, 1)
    mul!(m, Kxx′', cx, 1, 1)
    mul!(m, Fx′', d, 1, 1)

    return m
end

 # This is q
energy(cu, cx, y) = @views dot(cu, y[1:length(cu)]) + dot(cx, y[length(cu)+1:end])

# Returns vectors m, (qC)^-1 for test points
function predict_dist!(m, σ, ks::KernelSystem, cu, cx, d, q, X′, Fx′)
    W′ = similar(Fx′)
    Kux′ = Array{Float64}(undef, size(Fx′)...)
    Kxx′ = Array{Float64}(undef, size(ks.sd.X, 2), size(X′, 2))
    Dx′ = similar(σ)

    ldiv!(W′, ks.sd.Fu, Fx′)
    cross_correlation!(Kux′, ks.kd.kern, ks.sd.U, X′)
    cross_correlation!(Kxx′, ks.kd.kern, ks.sd.X, X′)

    mul!(m, Kux′', cu, 1, 1)
    mul!(m, Kxx′', cx, 1, 1)
    mul!(m, Fx′', d, 1, 1)

    diag_reduced_schur!(Dx′, Kxx′, ks.kd, Kux′, ks.sd.W, W′)
    σ .= inv.(Dx′ .* q)

    return m, σ
end

max_stars_repo_name: dbindel/btg

max_issues_repo_name: dbindel/btg

max_forks_repo_name: dbindel/btg


text:
using Stheno: Parameter, value
using Bijectors
using Random

rng = MersenneTwister(123456)

@testset "parameter wrapper" begin
    @testset "bijector" begin
        # scalar
        l = rand(rng)
        # vector
        N = 5
        v = rand(rng, N)
        # matrix
        M, M′ = 4, 7
        W = rand(rng, M, M′)

        # identity bijector
        pl = Parameter(l)
        @test l == value(pl)
        pv = Parameter(v)
        @test v == value(pv)
        pW = Parameter(W)
        @test W == value(pW)

        # Exp/Log bijector
        pl_pos = Parameter(l, Val(:pos))
        @test first(pl_pos.x) == Bijectors.Log{0}()(l)
        @test value(pl_pos) == Bijectors.Exp{0}()(Bijectors.Log{0}()(l))
        pv_pos = Parameter(v, Val(:pos))
        @test pv_pos.x == Bijectors.Log{1}()(v)
        @test value(pv_pos) == Bijectors.Exp{1}()(Bijectors.Log{1}()(v))
        pW_pos = Parameter(W, Val(:pos))
        @test pW_pos.x == Bijectors.Log{2}()(W)
        @test value(pW_pos) == Bijectors.Exp{2}()(Bijectors.Log{2}()(W))

        # check input domain
        @test_throws DomainError Parameter(-3.0, Val(:pos))
        @test_throws DomainError Parameter(9.0, inv(Bijectors.Logit(1.3, 6.0)))
    end
end


text:
using PyPlot
using Optim

A(ϕ) = acosh( 2ϕ / ( (4ϕ-π) * sin(ϕ) ) )
B(ϕ) = ( π^2 - 2π*ϕ ) / A(ϕ)

res = optimize(ϕ -> -B(ϕ), π/4, π/2)
ϕ_opt = res.minimizer

δ = 1e-3
ϕ = range(π/4+δ, π/2-δ, length=201)
figure(1)
plot(ϕ, B.(ϕ), [ϕ_opt], [B(ϕ_opt)], "o")
xticks((π/4, 3π/8, π/2),(L"$\pi/4$", L"$3π/8$", L"$\pi/2$"),fontsize=12)
grid(true)
savefig("Bplot.pdf")

println("Optimal values:")
println("\tϕ = ", ϕ_opt)
println("\th = $(A(ϕ_opt)) / N")
println("\tμ = $((4π*ϕ_opt-π^2)/A(ϕ_opt)) N")
println("\tB(ϕ) = ", B(ϕ_opt))
println("\tDecay factor = $(exp(B(ϕ_opt)))^(-N)")


text: include("SAA/SAA.jl")


text:
using Documenter
using Glob

using Distributions
using HMMBase
using Random

Random.seed!(2019)

# https://discourse.julialang.org/t/deactivate-plot-display-to-avoid-need-for-x-server/19359
ENV["GKSwstype"] = "nul"

# Documenter wants relative paths
examples = map(glob("*.md", joinpath(@__DIR__, "src", "examples"))) do x
    joinpath("examples", basename(x))
end

makedocs(
    sitename = "HMMBase",
    modules = [HMMBase],
    format = Documenter.HTML(assets = ["assets/goatcounter.js"]),
    pages = [
        "index.md",
        "Manual" => ["basics.md", "models.md", "algorithms.md", "utilities.md"],
        "Examples" => examples,
        "internals.md",
        "migration.md",
        "_index.md",
    ],
)

deploydocs(repo = "github.com/maxmouchet/HMMBase.jl.git")


text:
@testset "Morphisms" begin
    @testset "Find Automorphisms" begin
        G, = generic_group([1, -1, im, -im], *)
        @test find_small_group(generic_group(automorphisms(G), *)[1])[1] == (2,1)
        @test find_small_group(generic_group(automorphisms(small_group(8,2)), *)[1])[1] == (8,3)
        @test find_small_group(generic_group(automorphisms(small_group(13,1)), *)[1])[1] == (12,2)
        @test find_small_group(generic_group(automorphisms(small_group(18,5)), *)[1])[1] == (48,29)
        @test find_small_group(generic_group(automorphisms(small_group(24,12)), *)[1])[1] == (24,12)
        @test find_small_group(generic_group(automorphisms(small_group(33,1)), *)[1])[1] == (20,5)
        @test find_small_group(generic_group(automorphisms(small_group(36,2)), *)[1])[1] == (12,5)
        @test find_small_group(generic_group(automorphisms(small_group(42,3)), *)[1])[1] == (36,12)

        @test length(automorphisms(small_group(27,3))) == 432

        for aut in automorphisms(small_group(25,2))
            @test is_bijective(aut)
        end

        for aut in automorphisms(small_group(34,2))
            @test is_bijective(aut)
        end

        G = small_group(16,4)
        auts = automorphisms(G)
        for aut in auts
            @test all([aut(g * h) == aut(g) * aut(h) for g in G for h in G]) == true
        end

        G = small_group(20,2)
        auts = automorphisms(G)
        for aut in auts
            @test all([aut(g * h) == aut(g) * aut(h) for g in G for h in G]) == true
        end

        G = small_group(39,1)
        auts = automorphisms(G)
        for aut in auts
            @test all([aut(g * h) == aut(g) * aut(h) for g in G for h in G]) == true
        end

        G = small_group(44,2)
        auts = automorphisms(G)
        for aut in auts
            @test all([aut(g * h) == aut(g) * aut(h) for g in G for h in G]) == true
        end
    end

    @testset "Find Morphisms" begin
        G, = generic_group([1, -1, im, -im], *)
        H = small_group(6,2)
        @test length(morphisms(G,H)) == 2

        @test length(morphisms(small_group(23,1), small_group(42,5))) == 1
        @test length(morphisms(small_group(24,4), small_group(24,4))) == 124
        @test length(morphisms(small_group(24,4), small_group(45,2))) == 1
        @test length(morphisms(small_group(12,2), small_group(6,1))) == 6
        @test length(morphisms(small_group(56,4), small_group(24,2))) == 8
        @test length(morphisms(small_group(24,2), small_group(56,4))) == 32
        @test length(morphisms(small_group(24,2), small_group(58,2))) == 2
        @test length(morphisms(small_group(32,2), small_group(60,2))) == 64
        @test length(morphisms(small_group(56,4), small_group(62,2))) == 4
        @test length(morphisms(small_group(36,2), small_group(63,2))) == 9

        G = small_group(16,3)
        morphs = morphisms(G, G)
        for mor in morphs
            @test all([mor(g * h) == mor(g) * mor(h) for g in G for h in G]) == true
        end

        G = small_group(48,3)
        H = small_group(58,2)
        morphs = morphisms(G, H)
        for mor in morphs
            @test all([mor(g * h) == mor(g) * mor(h) for g in G for h in G]) == true
        end

        G = small_group(50,2)
        H = small_group(61,1)
        morphs = morphisms(G, H)
        for mor in morphs
            @test all([mor(g * h) == mor(g) * mor(h) for g in G for h in G]) == true
        end
    end

    @testset "morphisms" begin
        G,AtoG,GtoA = generic_group([1, -1, im, -im], *)
        Hom = GrpGenToGrpGenMor(G,G,[G[1],G[1],G[1],G[1]])
        @test order(image(Hom)[1]) == 1
        @test order(kernel(Hom)[1]) == 4
        @test is_surjective(Hom) == false
        @test is_injective(Hom) == false
        @test is_bijective(Hom) == false

        Hom = GrpGenToGrpGenMor(G,G,[G[1],G[2],G[3],G[4]])
        @test order(image(Hom)[1]) == 4
        @test order(kernel(Hom)[1]) == 1
        @test is_surjective(Hom) == true
        @test is_injective(Hom) == true
        @test is_bijective(Hom) == true
    end
end


text:
using CustomVertexGraphs
using Graphs
using Test

@testset "CustomVertexGraphs.jl" begin
  g = set_vertices(grid((4,)), ["A", "B", "C", "D"])
  @test vertices(g[["A", "B"]]) == ["A", "B"]
  @test has_edge(g, "A" => "B")
  add_edge!(g, "A" => "C")
  @test has_edge(g, "A" => "C")
end


text:
export Household,
    Introverted,
    Social,
    Enforced,
    Pike,
    Nutrients,
    Constant,
    Dynamic,
    Noise,
    TransientUp,
    TransientDown,
    Decision,
    Municipality,
    Governance,
    Planting,
    WastewaterTreatment,
    Trawling,
    Angling,
    VegetationImbalance,
    lake_initial_state,
    Clear,
    Turbid,
    X1,
    X2,
    X3,
    S1,
    S2,
    S3,
    T1,
    T2,
    T3,
    Experiment,
    LakeParameters

abstract type Intervention end
abstract type Knowledge end
abstract type Status end
abstract type Threshold end
abstract type HouseOwner end

abstract type NutrientSeries end

mutable struct Household <: AbstractAgent
    id::Int # The identifier number of the agent
    pos::Tuple{Int,Int} # The x, y location of the agent on a 2D grid
    compliance::Float64 # From willingness to upgrade
    oss::Bool # Is the sewage system upgraded or not?
    information::Bool # Has this agent been told about the upgrade?
    implementation_lag::Int # How long does it take the agent to choose to upgrade?
    municipality::Int # ID of the municipality this household belongs to
end

mutable struct Municipality <: AbstractAgent
    id::Int # The identifier number of the agent
    pos::Tuple{Int,Int} # The x, y location of the agent on a 2D grid
    name::String
    information::Bool # information on state of the lake
    legislation::Bool # legislation power to enforce new rules
    regulate::Bool
    respond_direct::Bool
    threshold_variable::Threshold
    interventions::Dict{Integer,Vector{Intervention}} # Set of interventions municipality will act on
    policies::Dict{Type{<:Intervention},NamedTuple} # Decision parameters, subject to change
    agents_uniform::Bool # TODO: This is a poorly named bool. Point here is that if this is true, the agents willingness_to_upgrade will be pulled form a uniform distribution.
    houseowner_type::HouseOwner
    willingness_to_upgrade::Float64
    tolerance_level_affectors::Float64
    neighbor_distance::Int
    households::Int
    affectors::Int
end

@with_kw mutable struct Governance
    # Properties
    information::Bool = false # information on state of the lake
    legislation::Bool = false # legislation power to enforce new rules
    # Actions
    regulate::Bool = true
    respond_direct::Bool = false
    threshold_variable::Threshold = Nutrients()
    interventions::Dict{Integer,Vector{Intervention}} = Dict(-1 => [WastewaterTreatment()]) # Set of interventions municipality will act on
    # Additions for decision making. Subject to change
    policies::Dict{Type{<:Intervention},NamedTuple} =
        Dict{Type{<:Intervention},NamedTuple}()
    # Related to home owners
    agents_uniform::Bool = false # TODO: This is a poorly named bool. Point here is that if this is true, the agents willingness_to_upgrade will be pulled form a uniform distribution.
    houseowner_type::HouseOwner = Introverted()
    willingness_to_upgrade::Float64 = 0.2
    tolerance_level_affectors::Float64 = 50.0
    neighbor_distance::Int = 3
end

struct Idle <: Status end
struct Running <: Status end
struct Complete <: Status end

struct Pike <: Threshold end
struct Nutrients <: Threshold end

struct Introverted <: HouseOwner end
struct Social <: HouseOwner end
struct Enforced <: HouseOwner end

"""
    Constant()

Nutrient level remains constant at the level of `init_nutrients`.
"""
struct Constant <: NutrientSeries end

"""
    Dynamic()

Nutrient runoff is managed by the municipality by incentivising households to upgrade
sewage systems that seep P into the lake.
"""
struct Dynamic <: NutrientSeries end

"""
    Noise(process, min, max)

Noise process given by DiffEqNoiseProcess.jl. For the moment this does not connect to
the actual start time or `init_nutrients` value, so these must be manually duplicated
here. Will be fixed in the future. `min` and `max` nutrient values can also be applied.

!!! warning

    Care must be taken to verify that bounded processes are possible over certain time
    horizons. This is noticably the case with brownian bridges. If your σ is small such
    that the gap between your current position and the end position cannot be met within
    the timeframe specified, the solver will error. Unfortunately this cannot be handled
    in a softer manner.
"""
@with_kw struct Noise <: NutrientSeries
    process::NoiseProcess = GeometricBrownianMotionProcess(0.0,0.05,0.0,2.0)
    min::Float64 = 0.0
    max::Float64 = 20.0
end

"""
    TransientUp(;start_year = 11, post_target_series = Constant())

Synthetic nutrient profile that alters lake dynamics regardless of municipal management.

- `start_year`: year when nutrients begin to increase with a rate of `nutrient_change`.
- `post_target_series`: behaviour after `target_nutrients` value is reached.
    Default is `TransientDown(start_year = 0, post_target_series = Constant())`

!!! warning

    Post target series selection must include a final `Constant` phase, otherwise an
    infinite recursion cascade will occur.
"""
@with_kw struct TransientUp <: NutrientSeries
    start_year::Integer = 11
    post_target_series::NutrientSeries =
        TransientDown(start_year = 0, post_target_series = Constant())
end

"""
    TransientDown(;start_year = 11, post_target_series = Constant())

Synthetic nutrient profile that alters lake dynamics regardless of municipal management.

- `start_year`: year when nutrients begin to decrease with a rate of `nutrient_change`.
- `post_target_series`: behaviour after `target_nutrients` value is reached.
    Default is `TransientUp(start_year = 0, post_target_series = Constant())`

!!! warning

    Post target series selection must include a final `Constant` phase, otherwise an
    infinite recursion cascade will occur.
"""
@with_kw struct TransientDown <: NutrientSeries
    start_year::Integer = 11
    post_target_series::NutrientSeries =
        TransientUp(start_year = 0, post_target_series = Constant())
end

"""
## Keywords

A few keywords that can be sent to the `bboptimize` routine have been made available
here:

- `max_time = 0.0`, a hard time limit for the optimiser to run (off by default).
- `max_steps = 10_000_000`, a hard step limit for the optimiser to run.
- `trace_mode = :compact`, logging output control. Other options are `:silent` and
`:verbose`.

Note that the standard stopping condition we envoke is `MaxStepsWithoutProgress = 10000`.
"""
@with_kw mutable struct Decision{F<:Function}
    start::Int = 1 # year when first optimisation is completed
    every::Int = 5 # year when next optimisation is completed (if target not met)
    current_term_only::Bool = true # If true, only optimise the next X years
    objectives::NTuple{N,Tuple{Function,Float64}} where {N} = ((min_time, 1.0),)
    target::F = clear_state
    # Optimiser settings
    max_time::Float64 = 0.0
    max_steps::Int = 10_000_000
    trace_mode::Symbol = :compact
    opt_replicates::Int = 0
    # Paralellism in BBO is currently not functioning, we don't need the workstealer anyway, so set it to empty
    #opt_pool = Agents.Distributed.WorkerPool()
end

# Properties of the experiment. For now this is a drop in for GUI values
# TODO: Some of these may need to be put under Governance
@with_kw mutable struct Experiment
    identifier::String = "biggs-baseline" # String for the moment. Might want to reconsider this in the future.
    # Related to lake
    # TODO: Perhaps the Water Council can set these values? Then all municipalities adhere to them.
    pike_expectation::Float64 = 1.4
    target_nutrients::Float64 = 0.7
    nutrient_series::NutrientSeries = Constant()
    nutrient_change::Float64 = 0.1
    critical_nutrients::Float64 = 3.0
    recycling_rate::Float64 = 0.1
    max_sewage_water::Float64 = 0.1
    # Additions for decision making. Subject to change
    policy::Decision = Decision()
    knowledge::Vector{Knowledge} = Knowledge[] # Current options, :vegetation_imbalance
end

@with_kw_noshow mutable struct Outcomes
    @deftype Int
    year_when_informed = 0# remembers the first year when municipality was informed about critical lake status from monitoring
    year_when_pike_became_critical = 0 # remember the first year pike is below expectation level
    year_when_nutrients_became_critical = 0 # remember the first year nutrients are above expectation level
    year_when_desired_pike_is_back = 0
    year_when_desired_level_is_back = 0
    upgraded_households_sum = 0 # cumulatively aggregates the number of upgraded oss each year
    year_of_full_upgrade = 0 #remember the year when all households finished their update
end

struct WastewaterTreatment <: Intervention end
@with_kw_noshow mutable struct Planting <: Intervention
    rate::Float64 = 1e-4
    cost::Float64 = 1.0 # Cost is given here as a ratio. Default is 1:1.
end
@with_kw_noshow mutable struct Trawling <: Intervention
    rate::Float64 = 5e-4
    cost::Float64 = 1.0
end
@with_kw_noshow mutable struct Angling <: Intervention
    rate::Float64 = 2.25e-4 # 10% of default rate
end

@with_kw_noshow mutable struct VegetationImbalance <: Knowledge
    bream_density_flip::Float64 = 40.0
    year_target_reached::Integer = 0
end

# Lake dynamics
abstract type LakeModel end
abstract type LakeParameters end
abstract type LakeDefinition end

struct Clear <: LakeDefinition end
struct Turbid <: LakeDefinition end
struct X1 <: LakeDefinition end
struct X2 <: LakeDefinition end
struct X3 <: LakeDefinition end

struct S1 <: LakeDefinition end
struct S2 <: LakeDefinition end
struct S3 <: LakeDefinition end
struct T1 <: LakeDefinition end
struct T2 <: LakeDefinition end
struct T3 <: LakeDefinition end

"""
    lake_initial_state(<:LakeDefinition, <:LakeModel)

Provides initial conditions of the lake. Parameters are model dependent.
For example: `lake_initial_state(Clear, Martin)`.
"""
function lake_initial_state(
    ::Type{L},
    ::Type{M};
    kwargs...,
) where {L<:LakeDefinition,M<:LakeModel}
    nutrients, initial_state = preset_conditions(L, M)
    (initial_state, LakeParameters(M, nutrients; kwargs...))
end


text:
####################
## Phantom OBJECT ##
####################
"""
    Phantom(x,y,ρ,T2,Δw,Dλ1,Dλ2,Dθ,ux,uy)

Phantom object.

# Arguments
 - `name::String` := Name of the Phantom
 - `x`            := Spins x-coordinates.
 - `y`            := Spins y-coordinates.
 - `z`            := Spins z-coordinates.
 - `ρ::Matrix`    := Proton density.
 - `T1::Matrix`   := T1 map.
 - `T2::Matrix`   := T2 map.
 - `Δw::Matrix`   := Off-resonace map;
 - `D::Matrix`    := Diffusion model.
 - `ux::Function` := Displacement field x.
 - `uy::Function` := Displacement field y.
 - `uz::Function` := Displacement field z.
"""
@with_kw mutable struct Phantom # TODO: in the future use @with_kw
	name::String = "spin" #Name of the Phantom
	x::Vector #x-coordinates of spins
	y::Vector = zeros(size(x)) #y-coordinates of spins
	z::Vector = zeros(size(x)) #z-coordinates of spins
	ρ::Vector = ones(size(x)) #proton density
	T1::Vector = Inf*ones(size(x)) #T1 map
	T2::Vector = Inf*ones(size(x)) #T2 map
	T2s::Vector = T2
	#Off-resonance related
	Δw::Vector = zeros(size(x)) #Off-resonace map
	#χ::Vector{SusceptibilityModel}
	#Diffusion
	Dλ1::Vector = zeros(size(x))
	Dλ2::Vector = zeros(size(x))
	Dθ::Vector = zeros(size(x))
	#Diff::Vector{DiffusionModel}  #Diffusion map
	#Motion
	ux::Function = (x,y,z,t)->0 #Displacement field x
	uy::Function = (x,y,z,t)->0 #Displacement field y
	uz::Function = (x,y,z,t)->0 #Displacement field z
end

# Phantom() = Phantom(name="spin",x=zeros(1,1))
size(x::Phantom) = size(x.ρ)
"""Separate object spins in a sub-group."""
Base.getindex(obj::Phantom, p::UnitRange{Int}) = begin
	Phantom(name=obj.name,
			x=obj.x[p],
			y=obj.y[p],
			z=obj.z[p],
			ρ=obj.ρ[p],
			T1=obj.T1[p],
			T2=obj.T2[p],
			T2s=obj.T2s[p],
			Δw=obj.Δw[p],
			#Diff=obj.Diff[p], #TODO!
			Dλ1=obj.Dλ1[p],
			Dλ2=obj.Dλ2[p],
			Dθ=obj.Dθ[p],
			#Χ=obj.Χ[p], #TODO!
			ux=obj.ux,
			uy=obj.uy,
			uz=obj.uz
			)
end
# Compartment enabling routines:
# Addition of compartments
+(s1::Phantom,s2::Phantom) =begin
	Phantom(name=s1.name*"+"*s2.name,
		x=[s1.x;s2.x],
		y=[s1.y;s2.y],
		z=[s1.z;s2.z],
		ρ=[s1.ρ;s2.ρ],
		T1=[s1.T1;s2.T1],
		T2=[s1.T2;s2.T2],
		T2s=[s1.T2s;s2.T2s],
		Δw=[s1.Δw;s2.Δw],
		#Diff=obj.Diff[p], #TODO!
		Dλ1=[s1.Dλ1;s2.Dλ1],
		Dλ2=[s1.Dλ2;s2.Dλ2],
		Dθ=[s1.Dθ;s2.Dθ],
		#Χ=obj.Χ[p], #TODO!
		ux=s1.ux,
		uy=s1.uy,
		uz=s1.uz
	)
end
#Fraction of compartments
*(α::Real,obj::Phantom) = begin
	Phantom(name=obj.name,
		x=obj.x,
		y=obj.y,
		z=obj.z,
		ρ=α*obj.ρ, #Only affects the proton density
		T1=obj.T1,
		T2=obj.T2,
		T2s=obj.T2s,
		Δw=obj.Δw,
		#Diff=obj.Diff[p], #TODO!
		Dλ1=obj.Dλ1,
		Dλ2=obj.Dλ2,
		Dθ=obj.Dθ,
		#Χ=obj.Χ[p], #TODO!
		ux=obj.ux,
		uy=obj.uy,
		uz=obj.uz
	)
end

# Movement related commands
# StartAt(s::Phantom,t0::Float64) = Phantom(s.name,s.x,s.y,s.ρ,s.T2,s.Δw,s.Dλ1,s.Dλ2,s.Dθ,(x,y,t)->s.ux(x,y,t.+t0),(x,y,t)->s.uy(x,y,t.+t0))
# FreezeAt(s::Phantom,t0::Float64) = Phantom(s.name*"STILL",s.x.+s.ux(s.x,s.y,t0),s.y.+s.uy(s.x,s.y,t0),s.ρ,s.T2,s.Δw,s.Dλ1,s.Dλ2,s.Dθ,(x,y,t)->0,(x,y,t)->0)

#TODO: jaw-pitch-roll, expand, contract, functions

# Getting maps
# get_DxDy2D(obj::Phantom) = begin
# 	P(i) = rotz(obj.Dθ[i])[1:2,1:2];
# 	D(i) = [obj.Dλ1[i] 0;0 obj.Dλ2[i]]
# 	nx = [1;0]; ny = [0;1]
# 	Dx = [nx'*P(i)'*D(i)*P(i)*nx for i=1:prod(size(obj.Dλ1))]
# 	Dy = [ny'*P(i)'*D(i)*P(i)*ny for i=1:prod(size(obj.Dλ1))]
# 	Dx, Dy
# end

"""
Heart-like LV phantom. The variable `α` is for strech, `β` for contraction, and `γ` for rotation.
"""
heart_phantom(α=1,β=1,γ=1,fat_bool::Bool=false) = begin
	#PARAMETERS
	FOV = 10e-2 #m Diameter ventricule
	N = 21
	Δxr = FOV/(N-1) #Aprox rec resolution, use Δx_pix and Δy_pix
	Ns = 50 #number of spins per voxel
	Δx = Δxr/sqrt(Ns) #spin separation
	#POSITIONS
	x = y = -FOV/2:Δx:FOV/2-Δx #spin coordinates
	x, y = x .+ y'*0, x*0 .+ y' #grid points
	#PHANTOM
	⚪(R) =  (x.^2 .+ y.^2 .<= R^2)*1. #Circle of radius R
	v = FOV/4 #m/s 1/16 th of the FOV during acquisition
	ωHR = 2π/1 #One heart-beat in one second

	# θ(t) = -π/4*γ*(sin.(ωHR*t).*(sin.(ωHR*t).>0)+0.25.*sin.(ωHR*t).*(sin.(ωHR*t).<0) )
	ux(x,y,t) = begin
		strech = 0 #α * v * (x.^2 .+ y.^2) / (FOV/2)^2 .* sign.(x)
		contract = - β * v * x / (FOV/2)  #expand
		rotate = - γ * v * y / (FOV/2)
		def = (strech .+ contract .+ rotate) .* sin.(ωHR*t)
	end
	uy(x,y,t) = begin
		strech = 0 #α * v * (x.^2 .+ y.^2) / (FOV/2)^2 .* sign.(y)
		contract = - β * v * y / (FOV/2)
		rotate = γ * v * x / (FOV/2)
		def = (strech .+ contract .+ rotate) .* sin.(ωHR*t)
	end
	# Water spins
	R = 9/10*FOV/2
	r = 6/11*FOV/2
	ring = ⚪(R) .- ⚪(r)
	ρ = ⚪(r) .+ 0.9*ring #proton density
	# Diffusion tensor model
	D = 2e-9 #Diffusion of free water m2/s
	D1, D2 = D, D/20
	Dλ1 = D1*⚪(R) #Diffusion map
	Dλ2 = D1*⚪(r) .+ D2*ring #Diffusion map
	Dθ =  atan.(x,-y) .* ring #Diffusion map
	T1 = (1400*⚪(r) .+ 1026*ring)*1e-3 #Myocardial T1
	T2 = ( 308*⚪(r) .+ 42*ring  )*1e-3 #T2 map [s]
	# Generating Phantoms
	heart = Phantom("LeftVentricle",x,y,ρ,T2,Dλ1,Dλ2,Dθ,ux,uy)
	# Fat spins
	ring2 = ⚪(FOV/2) .- ⚪(R) #outside fat layer
	ρ_fat = .5*ρ.*ring2
	Δw_fat = 2π*220*ring2 #fat should be dependant on B0
	T1_fat = 800*ring2*1e-3
	T2_fat = 120*ring2*1e-3 #T2 map [s]
	fat = Phantom("fat",x,y,ρ_fat,T2_fat,Δw_fat)
	#Resulting phantom
	obj = fat_bool ? heart + fat : heart #concatenating spins
end

"""
B. Aubert-Broche, D.L. Collins, A.C. Evans: "A new improved version of the realistic digital brain phantom"
NeuroImage, in review - 2006.
B. Aubert-Broche, M. Griffin, G.B. Pike, A.C. Evans and D.L. Collins: "20 new digital brain phantoms for creation of validation image data bases"
IEEE TMI, in review - 2006

https://brainweb.bic.mni.mcgill.ca/brainweb
"""
function brain_phantom2D(;axis="axial",ss=4)
    path = @__DIR__
    data = MAT.matread(path*"/data/brain2D.mat")

    class = data[axis][1:ss:end,1:ss:end]
    Δx = 1e-3*ss
    M, N = size(class)
    FOVx = (M-1)*Δx #[m]
    FOVy = (N-1)*Δx #[m]
    x = -FOVx/2:Δx:FOVx/2 #spin coordinates
    y = -FOVy/2:Δx:FOVy/2 #spin coordinates
    x, y = x .+ y'*0, x*0 .+ y' #grid points

    T2 = (class.==23)*329 .+ #CSF
        (class.==46)*83 .+ #GM
        (class.==70)*70 .+ #WM
        (class.==93)*70 .+ #FAT1
        (class.==116)*47 .+ #MUSCLE
        (class.==139)*329 .+ #SKIN/MUSCLE
        (class.==162)*0 .+ #SKULL
        (class.==185)*0 .+ #VESSELS
        (class.==209)*70 .+ #FAT2
        (class.==232)*329 .+ #DURA
        (class.==255)*70 #MARROW
    T2s = (class.==23)*58 .+ #CSF
        (class.==46)*69 .+ #GM
        (class.==70)*61 .+ #WM
        (class.==93)*58 .+ #FAT1
        (class.==116)*30 .+ #MUSCLE
        (class.==139)*58 .+ #SKIN/MUSCLE
        (class.==162)*0 .+ #SKULL
        (class.==185)*0 .+ #VESSELS
        (class.==209)*61 .+ #FAT2
        (class.==232)*58 .+ #DURA
        (class.==255)*61 .+#MARROW
        (class.==255)*70 #MARROW
    T1 = (class.==23)*2569 .+ #CSF
        (class.==46)*833 .+ #GM
        (class.==70)*500 .+ #WM
        (class.==93)*350 .+ #FAT1
        (class.==116)*900 .+ #MUSCLE
        (class.==139)*569 .+ #SKIN/MUSCLE
        (class.==162)*0 .+ #SKULL
        (class.==185)*0 .+ #VESSELS
        (class.==209)*500 .+ #FAT2
        (class.==232)*2569 .+ #DURA
        (class.==255)*500 #MARROW
    ρ = (class.==23)*1 .+ #CSF
        (class.==46)*.86 .+ #GM
        (class.==70)*.77 .+ #WM
        (class.==93)*1 .+ #FAT1
        (class.==116)*1 .+ #MUSCLE
        (class.==139)*1 .+ #SKIN/MUSCLE
        (class.==162)*0 .+ #SKULL
        (class.==185)*0 .+ #VESSELS
        (class.==209)*.77 .+ #FAT2
        (class.==232)*1 .+ #DURA
        (class.==255)*.77 #MARROW
	Δw_fat = -220*2π/10
	Δw = (class.==93)*Δw_fat .+ #FAT1
		(class.==209)*Δw_fat    #FAT2
	T1 = T1*1e-3
	T2 = T2*1e-3
	T2s = T2s*1e-3
    phantom = Phantom(name="brain2D_"*axis,
					  x=x[ρ.!=0],
					  y=y[ρ.!=0],
					  z=0*x[ρ.!=0],
					  ρ=ρ[ρ.!=0],
					  T1=T1[ρ.!=0],
					  T2=T2[ρ.!=0],
					  T2s=T2s[ρ.!=0],
					  Δw=Δw[ρ.!=0])
	phantom
end


text:
function build_master(sub_problems::T where T <: Dict, tree::T where T <: AbstractTree, probabilities::Dict{AbstractTree,Float64}, solver, discount_factor::Float64,
  risk::Union{Risk,Array{Risk,1}}, sideconstraints)

	if typeof(solver) <: Tuple
		model = Model(solver[1])
	else
		model = Model(solver)
	end
	@objective(model,Min,0)

	model.ext[:columns] = Dict{AbstractTree,Array{Column,1}}()
	for node in collect(tree)
		model.ext[:columns][node]=Array{Column,1}()
	end
	leafs=Array{Leaf,1}()

	model.ext[:expansions] = Dict{AbstractTree,Dict{Symbol,Any}}()

	model.ext[:options]=Dict{Symbol,Tuple}()
	for (name,options) in sub_problems[tree].ext[:options]
		model.ext[:options][name]=options
	end

	form=nothing
	for (name,options) in model.ext[:options]
		if options[4]==:Bin
			if form==:continuous
				form=:mixed
				break
			else
				form=:binary
			end
		elseif options[4]==:Int
			form=:mixed
			break
		else
			if form==:binary
				form=:mixed
				break
			else
				form=:continuous
			end
		end
	end

	model.ext[:discrete_con]=Dict{AbstractTree,Dict{Int64,ConstraintRef}}()
	model.ext[:discrete_var]=Dict{AbstractTree,Dict{Int64,VariableRef}}()

	for (node,sp) in sub_problems
		if haskey(node.ext,:sum_max) || haskey(node.ext,:sum_min) || haskey(node.ext,:sum_set)
			sp.ext[:form]=:binarycolumns
		else
			sp.ext[:form]=form
		end
		if sp.ext[:form] ∉ [:binary,:binarycolumns]
			vars = all_variables(sp)
			if form==:continuous
				for var in vars
					if is_binary(var) || is_integer(var)
						sp.ext[:form]=:mixed
						break
					end
				end
			end

			if sp.ext[:form]==:mixed
				sp.ext[:discrete]=VariableRef[]
				sp.ext[:discrete_branch]=VariableRef[]
				ignore=VariableRef[]

				for (name,variable) in sp.ext[:expansions]
					if typeof(variable) <: AbstractArray
						for index in eachindex(variable)
							push!(ignore,variable[index])
						end
					else
						push!(ignore,variable)
					end
				end

				for var in vars
					if (is_binary(var) || is_integer(var)) && var ∉ ignore
						push!(sp.ext[:discrete],var)
						push!(sp.ext[:discrete_branch],var)
					end
				end
			end
		end
	end

	for (node,sp) in sub_problems
		if typeof(node)==Leaf
			push!(leafs,node)
		end
		model.ext[:expansions][node] = Dict{Symbol,Any}()

		for (name,variable) in sp.ext[:expansions]
			if sp.ext[:options][name][4]==:Bin
				model.ext[:expansions][node][name] = copy_variable!(model, variable, relaxbinary)
			elseif sp.ext[:options][name][4]==:Int
				model.ext[:expansions][node][name] = copy_variable!(model, variable, relaxinteger)
			else
				model.ext[:expansions][node][name] = copy_variable!(model, variable)
			end
			if model.ext[:options][name][5]!=nothing
				if typeof(variable) <: AbstractArray
					for i in eachindex(variable)
						set_lower_bound(model.ext[:expansions][node][name][i],model.ext[:options][name][5])
					end
				else
					set_lower_bound(model.ext[:expansions][node][name],model.ext[:options][name][5])
				end
			end
			if model.ext[:options][name][6]!=nothing
				if typeof(variable) <: AbstractArray
					for i in eachindex(variable)
						set_upper_bound(model.ext[:expansions][node][name][i],model.ext[:options][name][6])
					end
				else
					set_upper_bound(model.ext[:expansions][node][name],model.ext[:options][name][6])
				end
			end
		end
	end

	model.ext[:scenprofit_var] = Dict{Leaf,VariableRef}()
	model.ext[:scenprofit_con] = Dict{Leaf,ConstraintRef}()

	for leaf in leafs
		model.ext[:scenprofit_var][leaf] = @variable(model)
		model.ext[:scenprofit_con][leaf] = @constraint(model, 0 == model.ext[:scenprofit_var][leaf])

  		set_objective_coefficient(model, model.ext[:scenprofit_var][leaf], probabilities[leaf])
	end

	risk_objectives=AffExpr[]
	remain=1.0
	if typeof(risk)==Risk
		if (risk.α==1.0 || risk.λ==0.0) && risk.bound==nothing
			risk=[]
		else
			risk=[risk]
		end
	end
    for i in 1:length(risk)
		risk_objective=AffExpr(0.0)
		eta=@variable(model)
	    for leaf in leafs
			v = @variable(model)
			w = @variable(model)
			set_lower_bound(v,0.0)
			set_lower_bound(w,0.0)
			if risk[i].offset==nothing
				@constraint(model,v>=eta-model.ext[:scenprofit_var][leaf])
				@constraint(model,w>=model.ext[:scenprofit_var][leaf]-eta)
				add_to_expression!(risk_objective,(model.ext[:scenprofit_var][leaf])*probabilities[leaf])
			else
				@constraint(model,v>=eta-model.ext[:scenprofit_var][leaf]+risk[i].offset[leaf])
				@constraint(model,w>=model.ext[:scenprofit_var][leaf]-risk[i].offset[leaf]-eta)
				add_to_expression!(risk_objective,(model.ext[:scenprofit_var][leaf]-risk[i].offset[leaf])*probabilities[leaf])
			end
			add_to_expression!(risk_objective,probabilities[leaf]*(v+w/risk[i].α*(1-risk[i].α)))
		end
		remain-=risk[i].λ
		push!(risk_objectives,risk_objective)
	end

	for leaf in leafs
		nodes=history(leaf)
		for n in eachindex(nodes)
			node=nodes[n]
			sp=sub_problems[node]
			df=discount_factor^depth(node)
			for (name,variable) in sp.ext[:expansions]
				interval=max(1,n-sp.ext[:options][name][3]-sp.ext[:options][name][2]+1):n-sp.ext[:options][name][2]
				disc=Dict{Int64,Float64}()
				for i in interval
					disc[i]=df/discount_factor^(i-1)
				end
				if typeof(variable) <: AbstractArray
					for i in eachindex(variable)
						cost_coef=df*coef(sp.ext[:capitalcosts],variable[i])
						for j in interval
							cost_coef+=disc[j]*coef(sp.ext[:ongoingcosts],variable[i])
						end
						set_normalized_coefficient(model.ext[:scenprofit_con][leaf], model.ext[:expansions][node][name][i], cost_coef)
					end
				else
					cost_coef=df*coef(sp.ext[:capitalcosts],variable)
					for j in interval
						cost_coef+=disc[j]*coef(sp.ext[:ongoingcosts],variable)
					end
					set_normalized_coefficient(model.ext[:scenprofit_con][leaf], model.ext[:expansions][node][name], cost_coef)
				end
			end
		end
	end

	# create the cover constraints
	model.ext[:coverconstraint] = Dict{AbstractTree,Dict{Symbol,Any}}()
	for (node,sp) in sub_problems
		model.ext[:coverconstraint][node] = Dict{Symbol,Any}()
		past=history(node)
		for (name,variable) in sp.ext[:expansions]
			interval=sp.ext[:options][name][2]+1:min(sp.ext[:options][name][2]+sp.ext[:options][name][3],length(past))
			if typeof(variable) <: AbstractArray
				model.ext[:coverconstraint][node][name] = Dict()
				for i in eachindex(variable)
					if sp.ext[:options][name][1]==:shutdown
						model.ext[:coverconstraint][node][name][i] = @constraint(model, 0 >= sum(model.ext[:expansions][past[index]][name][i] for index in interval))
					elseif sp.ext[:options][name][1]==:expansion
						model.ext[:coverconstraint][node][name][i] = @constraint(model, 0 <= sum(model.ext[:expansions][past[index]][name][i] for index in interval))
					elseif sp.ext[:options][name][1]==:enforced
						model.ext[:coverconstraint][node][name][i] = @constraint(model, 0 == sum(model.ext[:expansions][past[index]][name][i] for index in interval))
					elseif sp.ext[:options][name][1]==:state
						if node.parent==nothing
							model.ext[:coverconstraint][node][name][i] = @constraint(model, 0 == -sp.ext[:options][name][7] + model.ext[:expansions][node][name][i])
						else
							model.ext[:coverconstraint][node][name][i] = @constraint(model, 0 == model.ext[:expansions][node][name][i] - model.ext[:expansions][node.parent][name][i])
						end
					end
				end
				# if typeof(node)==Leaf && !sp.ext[:options][name][4]
				# 	for i in eachindex(variable)
				# 		@constraint(model, sum(model.ext[:expansions][past][name][i] for past in history_function(node))<=1)
				# 	end
				# end
			else
				if sp.ext[:options][name][1]==:shutdown
					model.ext[:coverconstraint][node][name] = @constraint(model, 0 >= sum(model.ext[:expansions][past[index]][name] for index in interval))
				elseif sp.ext[:options][name][1]==:expansion
					model.ext[:coverconstraint][node][name] = @constraint(model, 0 <= sum(model.ext[:expansions][past[index]][name] for index in interval))
				elseif sp.ext[:options][name][1]==:enforced
					model.ext[:coverconstraint][node][name] = @constraint(model, 0 == sum(model.ext[:expansions][past[index]][name] for index in interval))
				elseif sp.ext[:options][name][1]==:state
					if node.parent==nothing
						model.ext[:coverconstraint][node][name] = @constraint(model, 0 == -sp.ext[:options][name][7] + model.ext[:expansions][node][name])
					else
						model.ext[:coverconstraint][node][name] = @constraint(model, 0 == model.ext[:expansions][node][name] - model.ext[:expansions][node.parent][name])
					end
				end
				# if typeof(node)==Leaf && !sp.ext[:options][name][4]
				# 	@constraint(model, sum(model.ext[:expansions][past][name] for past in history_function(node))<=1)
				# end
			end
		end
	end

	model.ext[:convexcombination] = Dict{AbstractTree,Array{ConstraintRef,1}}()
	model.ext[:mip]=false
	model.ext[:branch_cons]=ConstraintRef[]

	for node in keys(sub_problems)
		model.ext[:convexcombination][node]=ConstraintRef[]
		if haskey(node.ext,:sum_max)
			push!(model.ext[:convexcombination][node],@constraint(model, 0 <= node.ext[:sum_max]))
		end
		if haskey(node.ext,:sum_min)
			push!(model.ext[:convexcombination][node],@constraint(model, 0 >= node.ext[:sum_min]))
		end
		if haskey(node.ext,:sum_set)
			if length(model.ext[:convexcombination][node])>0
				error(":eq and :le/:ge specifed for node "*node.name)
			end
			push!(model.ext[:convexcombination][node],@constraint(model, 0 == node.ext[:sum_set]))
		end
		if length(model.ext[:convexcombination][node])==0
			push!(model.ext[:convexcombination][node],@constraint(model, 0 == 1))
		end
	end

	if typeof(sideconstraints) <: Function
		map(Main.eval,unpack_expansions(model.ext[:expansions])) #bring expansion variables into global scope
		sideconstraints(model,tree)
		map(Main.eval,clear_expansions(model.ext[:expansions]))
	end

	objective_fn=objective_function(model)*remain

	for i in 1:length(risk_objectives)
		objective_fn+=risk_objectives[i]*risk[i].λ
		if risk[i].bound!=nothing
			surplus=@variable(model)
			set_lower_bound(surplus,0)
			@constraint(model,risk_objectives[i]<=risk[i].bound+surplus)
			objective_fn+=risk[i].penalty*surplus
		end
	end

	set_objective_function(model,objective_fn)

	model
end

function solve_master_binary(judge::JuDGEModel,initial_time::Float64,termination::Termination,warm_starts::Bool,mp_callback,verbose)
	set=0
	current=nothing
	if (termination.allow_frac==:binary_solve || termination.allow_frac==:binary_solve_return_relaxation) #&& judge.log[end].int>0
		if verbose==2
			print("Solving master problem as MIP")
		end
		set=1
		if solve_binary(judge,termination.abstol,termination.reltol,warm_starts,mp_callback)
			current = ConvergenceState(judge.log[end].obj, judge.bounds.UB, judge.bounds.LB, time()-initial_time, judge.log[end].iter + 1, 0)
			if verbose>0
				if verbose==2
					overprint("")
				end
				display(current,relaxation=false)
			end
			push!(judge.log,current)
		else
			if verbose==2
				overprint("")
			end
		end
	end
	if termination.allow_frac==:binary_solve_return_relaxation && set==1
		remove_binary(judge)
		optimize!(judge.master_problem)
	end
	current
end

function solve_binary(judge::JuDGEModel,abstol::Float64,reltol::Float64,warm_starts::Bool,mp_callback)
	if typeof(judge.master_solver) <: Tuple
		JuMP.set_optimizer(judge.master_problem,judge.master_solver[2])
	end

	judge.master_problem.ext[:mip]=true

	for node in keys(judge.master_problem.ext[:expansions])
		if judge.sub_problems[node].ext[:form]==:binarycolumns
			for col in judge.master_problem.ext[:columns][node]
			 	if !haskey(col.node.ext,:col_max)
			 		set_binary(col.var)
				else
			 		set_integer(col.var)
			 	end
			 end
		 end

		for x in keys(judge.master_problem.ext[:expansions][node])
			var = judge.master_problem.ext[:expansions][node][x]
			if judge.sub_problems[judge.tree].ext[:options][x][4]==:Bin
				if typeof(var) <: JuMP.Containers.SparseAxisArray
					for key in keys(var.data)
						set_binary(var.data[key])
					end
				elseif typeof(var) <: AbstractArray
					for key in keys(var)
						set_binary(var[key])
					end
				else
					set_binary(var)
				end
			elseif judge.sub_problems[judge.tree].ext[:options][x][4]==:Int
				if typeof(var) <: JuMP.Containers.SparseAxisArray
					for key in keys(var.data)
						set_integer(var.data[key])
					end
				elseif typeof(var) <: AbstractArray
					for key in keys(var)
						set_integer(var[key])
					end
				else
					set_integer(var)
				end
			end
		end
		if judge.sub_problems[node].ext[:form]==:mixed
			for (i,var) in judge.master_problem.ext[:discrete_var][node]
				set_integer(var)
			end
		end
	end

	if typeof(mp_callback) <: Function
		mp_callback(judge,abstol,reltol)
		printright("Initialising")
	end
	optimize!(judge.master_problem)
	status=termination_status(judge.master_problem)
	if status!=MathOptInterface.INFEASIBLE_OR_UNBOUNDED && status!=MathOptInterface.INFEASIBLE && status!=MathOptInterface.DUAL_INFEASIBLE
		obj = objective_value(judge.master_problem)
		if obj<judge.bounds.UB
			judge.bounds.UB=obj
			if warm_starts
				vars=all_variables(judge.master_problem)
				set_start_value.(vars, JuMP.value.(vars))
			end
		end
		return true
	else
		return false
	end
end

function remove_binary(judge::JuDGEModel)
	# if judge.master_problem.ext[:binarycolumns]
	# 	for col in judge.master_problem.ext[:columns]
	# 		if !haskey(col.node.ext,:col_max)
	# 			unset_binary(col.var)
	# 		else
	# 			unset_integer(col.var)
	# 		end
	# 	end
	# end
	for node in keys(judge.master_problem.ext[:expansions])
		if judge.sub_problems[node].ext[:form]==:binarycolumns
			for col in judge.master_problem.ext[:columns][node]
			 	if !haskey(col.node.ext,:col_max)
			 		unset_binary(col.var)
				else
			 		unset_integer(col.var)
			 	end
			 end
		 end
		for x in keys(judge.master_problem.ext[:expansions][node])
			var = judge.master_problem.ext[:expansions][node][x]
			if judge.sub_problems[judge.tree].ext[:options][x][4]==:Bin
				if typeof(var) <: AbstractArray
					for key in eachindex(var)
						unset_binary(var[key])
					end
				else
					unset_binary(var)
				end
			elseif judge.sub_problems[judge.tree].ext[:options][x][4]==:Int
				if typeof(var) <: AbstractArray
					for key in eachindex(var)
						unset_integer(var[key])
					end
				else
					unset_integer(var)
				end
			end
		end
		if judge.sub_problems[node].ext[:form]==:mixed
			for (i,var) in judge.master_problem.ext[:discrete_var][node]
				unset_integer(var)
			end
		end
	end

	judge.master_problem.ext[:mip]=false

	if typeof(judge.master_solver) <: Tuple
		JuMP.set_optimizer(judge.master_problem,judge.master_solver[1])
	end
end


text:
module Common
  export hello
  hello() = error("Common stub method.")
end

module X
  importall Common
  hello(s1::AbstractString, s2::AbstractString) = "Hello $(s1) $(s2)"
end

module Y
  importall Common
  hello(s1::AbstractString) = "Hello $(s1)"
end

module All
  using Common, X, Y

  println(hello("Peter", "Griffin"))
end


text:
#=

  From BLOG examples/healtiness.blog:
  """
  healthiness model
  This BLOG program describes the model for healthiness in
  a person.
  See exercise 3.5 in:
  Probabilistic Graphical Models: Principles and Techniques
  Daphne Koller, Nir Friedman, MIT 2009
  Each person is described by whether they are health
  conscious, have free time, exercise, have a good diet,
  have normal weight, have high cholesterol, and whether
  they tested positive for high cholesterol.
  @author: jnieh
  @date: 2012-09-13
  """

  For evidence TestedHighCholesterol == false:
  Summary Statistics
             parameters      mean       std   naive_se      mcse          ess      rhat   ess_per_sec
                 Symbol   Float64   Float64    Float64   Float64      Float64   Float64       Float64

        HealthConscious    0.5614    0.4962     0.0050    0.0080    4159.9755    0.9999      405.4952
         LittleFreeTime    0.4870    0.4999     0.0050    0.0081    4242.7390    0.9999      413.5626
               Exercise    0.5430    0.4982     0.0050    0.0075    4070.4820    1.0002      396.7718
               GoodDiet    0.6641    0.4723     0.0047    0.0074    4392.4673    1.0001      428.1575
           NormalWeight    0.5620    0.4962     0.0050    0.0069    4156.5796    1.0010      405.1642
        HighCholesterol    0.0987    0.2983     0.0030    0.0043    4344.1532    0.9999      423.4480
  TestedHighCholesterol    0.0001    0.0100     0.0001    0.0001   10000.0008    1.0000      974.7540

  For evidence exercise = true
  Summary Statistics
             parameters      mean       std   naive_se      mcse         ess      rhat   ess_per_sec
                 Symbol   Float64   Float64    Float64   Float64     Float64   Float64       Float64

        HealthConscious    0.6805    0.4663     0.0047    0.0071   4317.9635    1.0006      425.9607
         LittleFreeTime    0.2970    0.4570     0.0046    0.0070   4361.5047    1.0000      430.2560
               Exercise    1.0000    0.0000     0.0000    0.0000         NaN       NaN           NaN
               GoodDiet    0.5774    0.4940     0.0049    0.0072   4022.1381    1.0000      396.7780
           NormalWeight    0.6787    0.4670     0.0047    0.0074   3908.1858    0.9999      385.5367
        HighCholesterol    0.4715    0.4992     0.0050    0.0078   4117.2828    1.0004      406.1638
  TestedHighCholesterol    0.4780    0.4995     0.0050    0.0078   4117.0265    0.9999      406.1386

  Cf ~/webppl/healthiness.wppl

=#

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

@model function healthiness()

    HealthConscious ~ flip(0.5)
    LittleFreeTime ~ flip(0.5)

    Exercise ~
            (HealthConscious==true && LittleFreeTime==true) ? flip(0.5) :
            (HealthConscious==true && LittleFreeTime== false) ? flip(0.9) :
            (HealthConscious==false && LittleFreeTime== true) ? flip(0.1) :
            (HealthConscious==false && LittleFreeTime== false) ? flip(0.5) : Dirac(false)

    GoodDiet ~ HealthConscious ? flip(0.7) : flip(0.3)

    NormalWeight ~ (GoodDiet==true && Exercise==true) ? flip(0.8) :
                   (GoodDiet==true && Exercise== false) ? flip(0.5) :
                   (GoodDiet==false && Exercise== true) ? flip(0.5) :
                   (GoodDiet==false && Exercise== false) ? flip(0.2) : Dirac(false)

    HighCholesterol ~ GoodDiet ? flip(0.3) : flip(0.7)

    TestedHighCholesterol ~ HighCholesterol ? flip(0.9) : flip(0.1)

    # Evidence
    true ~ Dirac(TestedHighCholesterol == false) #  Original evidence
    # true ~ Dirac(TestedHighCholesterol == true)
    # true ~ Dirac(Exercise == true)
    # true ~ Dirac(LittleFreeTime == true)

end

model = healthiness()

num_chains = 4

# chns = sample(model, Prior(), 10_000)
# chns = sample(model, MH(), 10_000)
chns = sample(model, PG(5), 10_000)
# chns = sample(model, SMC(), 10_000)
# chns = sample(model, IS(), 10_000)

display(chns)
# display(plot(chns))

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank