text:
using Setfield, Documenter, Literate

inputdir = joinpath(@__DIR__, "..", "examples")
outputdir = joinpath(@__DIR__, "src", "examples")
mkpath(outputdir)
for filename in readdir(inputdir)
    inpath = joinpath(inputdir, filename)
    Literate.markdown(inpath, outputdir; documenter=true)
end

makedocs(
         modules = [Setfield],
         sitename = "Setfield.jl",
         pages = [
            "Introduction" => "intro.md",
            "Docstrings" => "index.md",
            "Custom Macros" => "examples/custom_macros.md",
            ],
        strict = true,  # to exit with non-zero code on error
        )

deploydocs(
    repo = "github.com/jw3126/Setfield.jl.git",
)


text:
using JuMP, MathProgBase, CPLEX, Random, LinearAlgebra, Mosek, JLD, LaTeXStrings, DataFrames, Test, Suppressor, DelimitedFiles, CSV, StatsBase

include("types.jl")
include("miop_formulation.jl")
include("inner_problem.jl")
include("hillclimb.jl")
include("socp_relaxation.jl")
include("raw.jl")
include("kelley_primal.jl")

ext: jl


text:
@symbol_func function f_IN(cur_reactor::AbstractReactor)
  cur_f = 1.0

  cur_f -= f_BS(cur_reactor)

  cur_f -= f_CD(cur_reactor)

  cur_f
end


text:
# This file can be used to create the validation file for scc files. It will create
# a validation file containing all possibilities of parameter values defined in the
# specs dictionary below produce the same results. A subset of these values are
# tested in the "runtests.jl" file deployed by Travis, and the full set can be tested
# manually using the scc_validation_full.jl file

using MimiFUND
using DataFrames
using Query
using CSVFiles
using Test

specs = Dict([
    :gas => [:CO2, :CH4, :N2O, :SF6],
    :year => [2020, 2055],
    :eta => [0, 1.45],
    :prtp => [0.015, 0.03],
    :equity_weights => [true, false],
    :equity_weights_normalization_region => [0, 1, 10],
    :last_year => [2300, 3000],
    :pulse_size => [1, 1e7]
])

results = DataFrame(gas = [], year = [], eta = [], prtp = [], equity_weights = [], equity_weights_normalization_region = [], last_year = [], pulse_size = [], SC = [])

for gas in specs[:gas]
    for year in specs[:year]
        for eta in specs[:eta]
            for prtp in specs[:prtp]
                for equity_weights in specs[:equity_weights]
                    for equity_weights_normalization_region in specs[:equity_weights_normalization_region]
                        for last_year in specs[:last_year]
                            for pulse_size in specs[:pulse_size]
                                sc = MimiFUND.compute_sc(gas=gas, year=year, eta=eta, prtp=prtp, equity_weights=equity_weights, equity_weights_normalization_region=equity_weights_normalization_region, last_year=last_year, pulse_size=pulse_size)
                                push!(results, (gas, year, eta, prtp, equity_weights, equity_weights_normalization_region, last_year, pulse_size, sc))
                            end
                        end
                    end
                end
            end
        end
    end
end

path = joinpath(@__DIR__, "SC validation data", "deterministic_sc_values_v3-13-0.csv")
save(path, results)

## Compare v3.11.7 to v3.13.0 ... :CO2 and :CH4 should not change

old_sc = load(joinpath(@__DIR__, "SC validation data", "deterministic_sc_values_v3-11-7.csv")) |> DataFrame
new_sc = load(joinpath(@__DIR__, "SC validation data", "deterministic_sc_values_v3-13-0.csv")) |> DataFrame

filter!(:gas => x -> (x == "CO2" || x == "CH4"), old_sc)
filter!(:gas => x -> (x == "CO2" || x == "CH4"), new_sc)

# after running scc_validation_all.jl on Mimi v3.12.1 on Lisa Rennels' machine,
# we note that the differences between results and validation results pass only at
# a tolerance of 1.0e-1 for pulse_size of 1 and 1e-9 for pulse_size of 1.0e7, so
# we mirror that here

old_sc_bigpulse = filter(:pulse_size => x -> x == 1.0e7, old_sc)
new_sc_bigpulse = filter(:pulse_size => x -> x == 1.0e7, new_sc)
@test all(isapprox.(old_sc_bigpulse[!, :SC], new_sc_bigpulse[!, :SC], atol = 1e-9))

old_sc_smallpulse = filter(:pulse_size => x -> x == 1.0e7, old_sc)
new_sc_smallpulse = filter(:pulse_size => x -> x == 1.0e7, new_sc)
@test all(isapprox.(old_sc_smallpulse[!, :SC], new_sc_smallpulse[!, :SC], atol = 1e-1))


text:
export LShapedRoom

"""
`L-shaped cuboid room object`

Returns a `LShapedRoom` object containing the geometrical information of the room.

# Usage

* LShapedRoom(Lx::Float64,Ly::Float64,Lz::Float64,xedge::Float64,yedge::Float64,ξ::Array{Float64,1}, env::FDTDEnv)
* LShapedRoom(Nx::Int64,Ny::Int64,Nz::Int64,xedge::Int64,yedge::Int64,ξ::Array{Float64,1}, env::FDTDEnv)
  * `ξ::Array{Float64}`: must be 6 element Array containg acoustic impedance of the 6 walls
  * `xedge` and `yedge` indicates the position of the edges, can also be expressed in samples

"""
immutable LShapedRoom <: AbstractGeometry
	Lx::Float64   #x dimension
	Ly::Float64   #y dimension
	Lz::Float64   #z dimension

	Nx::Int64   #x-axis samples
	Ny::Int64   #y-axis samples
	Nz::Int64   #z-axis samples

	xedge::Float64  #edge x position
	yedge::Float64  #edge y position

	G::Array{Bool,3}      #Geometry
	ξ::Array{Float64,1}   #Impedance
	env::FDTDEnv          #FDTD acoustic env.
end

function LShapedRoom(Lx::Float64,Ly::Float64,Lz::Float64,
		     xedge::Float64,yedge::Float64,
		    ξ::Array{Float64,1}, env::FDTDEnv)

	if(any([Lx;Ly;Lz].< 0)) error("room dimensions L should be positive") end
	if(length(ξ)!= 6) error("length(ξ) must be 6") end

	Nx = round(Int64,Lx/env.X)
	Ny = round(Int64,Ly/env.X)
	Nz = round(Int64,Lz/env.X)

	Ex = round(Int64,xedge/env.X)
	Ey = round(Int64,yedge/env.X)

	G = CreateGeometryLshaped(Nx,Ny,Nz,Ex,Ey)

	LShapedRoom(Lx,Ly,Lz,Nx,Ny,Nz,xedge,yedge,G,ξ,env)
end

function LShapedRoom(Nx::Int64,Ny::Int64,Nz::Int64,
		     Ex::Int64,Ey::Int64,
		    ξ::Array{Float64,1}, env::FDTDEnv)

	if(any([Nx;Ny;Nz].< 0)) error("room dimensions L should be positive") end
	if(length(ξ)!= 6) error("length(ξ) must be 6") end

	Lx =Nx*env.X
	Ly =Ny*env.X
	Lz =Nz*env.X

	xedge =Ex*env.X
	yedge =Ey*env.X
	G = CreateGeometryLshaped(Nx,Ny,Nz,Ex,Ey)

	LShapedRoom(Lx,Ly,Lz,Nx,Ny,Nz,xedge,yedge,G,ξ,env)
end

function CreateGeometryLshaped(Nx::Int64,Ny::Int64,Nz::Int64,
			       Ex::Int64,Ey::Int64)
	#=
        Create cuboid geometry
        inputs: Nx, number of samples x direction
        Ny, number of samples y direction
        Nz, number of samples z direction
        output: G, 3D tensor containing 1 inside domain and 0 outside
	=#

	if(Ex>Nx || Ey>Ny) error("edge outside room!") end

	G = ones(Bool,Nx,Ny,Nz)                 # creates the tensor with 1 where domain is present
	G[end-Ex:end,end-Ey:end,:] = 0     # this change geometry from pure cuboid

	G2 = zeros(Bool,Nx+2,Ny+2,Nz+2)# puts a shell of zeros around the tensor
	G2[2:Nx+1, 2:Ny+1, 2:Nz+1] = G

	return G2
end

fun_name(f::LShapedRoom) = "L-Shaped room"
fun_dim( f::LShapedRoom) = @sprintf("Lx x Ly x Lz = %.2f x %.2f x %.2f m^3", f.Lx, f.Ly, f.Lz)
fun_edge(f::LShapedRoom) = @sprintf("(x,y) = (%.2f, %.2f) m", f.xedge, f.yedge)
fun_Nxyz(f::LShapedRoom) = @sprintf("Nx x Ny x Nz = %d x %d x %d", f.Nx, f.Ny, f.Nz)
fun_ξ(   f::LShapedRoom) =
@sprintf("frequency independent   \n\n |  ξx1 |  ξx2 |  ξy1 |  ξy2 |  ξz1 | ξz2  | \n | % 02.1f| % 02.1f| % 02.1f| % 02.1f| % 02.1f| % 02.1f|\n" , f.ξ[1],f.ξ[2],f.ξ[3],f.ξ[4],f.ξ[5],f.ξ[6])
fun_X(  f::LShapedRoom) = @sprintf(" %.2f cm", f.env.X*100)
fun_Fs( f::LShapedRoom) = @sprintf(" %.2f kHz", f.env.Fs/1000)

function Base.show(io::IO, f::LShapedRoom)

	println(io, "geometry     : ", fun_name(f))
	println(io, "dimensions   : ", fun_dim(f))
	println(io, "edge at      : ", fun_edge(f))
	println(io, "samples      : ", fun_Nxyz(f))
	println(io, "ξ            : ", fun_ξ(f))
	println(io, "spatial step : ", fun_X(f))
	println(io, "sampling frq.: ", fun_Fs(f))

end

max_stars_repo_name: nantonel/JuFdtd

max_issues_repo_name: nantonel/JuFdtd

max_forks_repo_name: nantonel/JuFdtd


text:
using OpenTelemetryAPI
using Test

@testset "OpenTelemetryAPI.jl" begin

include("context.jl")
include("propagator.jl")
include("trace.jl")

end


text:
using Printf
using HTTP
using JSON

#include("white_backend.jl")
include("psql_backend.jl")
using Main.Backend

url = "https://api.telegram.org/"
key=ENV["TELEGRAM_KEY"]

function send_message(params)
    req = HTTP.request("POST",string(url,key,"/sendMessage"),["Content-Type" => "application/json"],JSON.json(params))
end

function send_welcome(contact)
    id = contact["id"]
    name = contact["first_name"]
    if "language_code" in keys(contact)
        lc = contact["language_code"]
    else
        lc = "it"
    end
    if lc=="it"
        params = Dict("chat_id"=>id, "text"=>string("Ciao ", name, "! Sono Pizzabot\n
        Ti aiuterò a preparare pizze e focacce favolose!
        Digita /cuciniamo per iniziare e poi segui le istruzioni.

        Buon lavoro e buon appetito!!"))
    else
        params = Dict("chat_id"=>id, "text"=>string("Welcome ", name, " \nClick /cuciniamo to begin.."))
    end
    send_message(params)
end

function ask_program(contact)
    id = contact["id"]
    buttons = [Dict("text"=>"pizza"), Dict("text"=>"focaccia")]
    reply_markup = Dict("keyboard"=>[buttons], "one_time_keyboard"=>true)
    params = Dict("chat_id"=>id, "text"=>"Cosa vuoi cucinare? (pizza o focaccia)", "reply_markup"=>reply_markup)
    send_message(params)
end

function ask_people(message)
    id = message["from"]["id"]
    program = message["text"]

    if !(program  in ["pizza" "focaccia"])
        program = "pizza"
    end

    params = Dict("chat_id"=>id, "text"=>string("Bene allora ti darò istruzioni per fare la ", program))
    send_message(params)
    params = Dict("chat_id"=>id, "text"=>string("per quante persone vuoi fare la ", program, "?"))
    send_message(params)
    return program
end

function read_people(message, state)
    id = message["from"]["id"]
    people = message["text"]
    program = get_key(id, "program")
    try
        people = parse(Int, people)
        if people <= 0
            throw("too few people")
        end
        send_recipe(message, people)
        return state + 1
    catch
        params = Dict("chat_id"=>id, "text"=>string("Non ho capito. Per quante persone?"))
        send_message(params)
        return state
    end
end

function send_recipe(message, people)
    id = message["from"]["id"]
    program = get_key(id, "program")
    if program == "pizza"
        send_pizza(id, people)
    elseif program == "focaccia"
        send_focaccia(id, people)
    end
end

function send_pizza(id, people)
    pizza_recipe = [1000 100 500]
    recipe = pizza_recipe / 4. * people
    params = Dict("chat_id"=>id, "text"=>string("Ecco le dosi per $(people) persone:\n
    $(@sprintf("%.1f", recipe[1]/1000)) Kg di farina\n
    $(@sprintf("%.1f", recipe[2])) ml di olio\n
    $(@sprintf("%.1f", recipe[3])) ml di acqua\n
    lievito q.b.\n
    sale q.b.\n"))
    send_message(params)
end

function send_focaccia(id, people)
    focaccia_recipe = [660 50 440]
    recipe = focaccia_recipe / 4. * people

    number, rectangular_str, circular_str, rectangular, circular = get_teglia(sum(recipe))

    params = Dict("chat_id"=>id, "text"=>string("Ecco le dosi per $(people) persone:\n
    $(@sprintf("%.3f", recipe[1]/1000)) Kg di farina di Manitoba (W=280)\n
    $(@sprintf("%.1f", recipe[2])) g di olio\n
    $(@sprintf("%.1f", recipe[3])) ml di acqua\n
    lievito q.b.\n
    sale q.b.\n\n
    Puoi usare $(circular_str) di diametro: $(@sprintf("%.1f", circular))cm\n
    oppure $(rectangular_str) $(@sprintf("%.1f", rectangular[1])) cm X $(@sprintf("%.1f", rectangular[2])) cm"))
    send_message(params)
end

function get_teglia(impasto)
    area = impasto / 0.60
    rapp = 42 / 35
    area_max = 42 * 35
    number = Int(round(floor(area/area_max) + 1))

    area = area / number # area singola teglia

    h = sqrt(area / rapp)
    b = area / h
    rectangular = [b h]
    circular = sqrt(area/pi)*2

    if number > 1
        rectangular_str = "$(@sprintf("%d", number)) teglie rettangolari"
        circular_str = "$(@sprintf("%d", number)) teglie circolari"
    else
        rectangular_str = "una teglia rettangolare"
        circular_str = "una teglia circolare"
    end

    return (number, rectangular_str, circular_str, rectangular, circular)
end

function send_bye(contact)
    id = contact["id"]
    name = contact["first_name"]
    if "language_code" in keys(contact)
        lc = contact["language_code"]
    else
        lc = "it"
    end
    if lc=="it"
        params = Dict("chat_id"=>id, "text"=>string("Buon appetito ", name, "!"))
    else
        params = Dict("chat_id"=>id, "text"=>string("Enjoy your meal ", name, "!"))
    end
    send_message(params)
end

function parse_message(message)
    id = message["message"]["from"]["id"]
    update_contact_list(message)
    state = get_state(id)

    if (state == 0) | (message["message"]["text"]=="/start")
        send_welcome(message["message"]["from"])
        set_state(id, 1)
    elseif (state == 1) | (message["message"]["text"]=="/cuciniamo")
        ask_program(message["message"]["from"])
        set_state(id, 2)
    elseif state == 2
        program = ask_people(message["message"])
        set_state(id, 3)
        set_key(id,"program",program)
    elseif state == 3
        new_state = read_people(message["message"], 3)
        state = new_state
        set_state(id, state)
    end

    if state == 4
        send_bye(message["message"]["from"])
        set_state(id, 1)
    end
end

function run()
    offset = -1

    println("ready!")
    while true
    #for i in 1:100
        if offset >= 0
            params = Dict("offset"=>offset)
        else
            params = Dict()
        end
        req = HTTP.request("GET", string(url,key,"/getUpdates"), ["Content-Type" => "application/json"],JSON.json(params))
        body = String(req.body)
        result = JSON.Parser.parse(body)

        if length(result["result"]) < 1
            continue
        end
        offset = result["result"][end]["update_id"] + 1
        for message in result["result"]
           parse_message(message)
        end
        sleep(0.1)
    end
end

init_db()
conn = get_conn()

run()

@sync close(conn)


text:
using Flux
using Statistics: mean, std
using Random
#import MultiTheisDP.MultiTheis
#using PyPlot
#ioff()
#using PyCall
#gridspec = pyimport("matplotlib.gridspec")
#using LaTeXStrings
using JLD2
#include("plot_train.jl")
#include("theislike.jl")
push!(LOAD_PATH,"../../src")
push!(LOAD_PATH,"./")
import TheisLike

########################################################
# User input
# Location to monitor pressure at
mon_well_crds = [[-75.0,-75.0]]
# Location of injection and extraction wells
well_crds = [[-50.0,-50.0],
			[50.0,50.0]]
# Injection rate
Qinj = 0.031688 # [m^3/s] (1 MMT water/yr)
# Time to check pressure
t = 30.0*86400.0 # 30 days converted to seconds
# Create xy grid for plotting pressures
xs = ys = -100.0:1.0:100.0
# Training epochs
epochs = 1:120
########################################################

# Calculate distance between extraction and injection wells
rs = [sqrt((sum((cs-mon_well_crds[1]).^2))) for cs in well_crds]

# Input/output normalization function
lerp(x, lo, hi) = x*(hi-lo)+lo

global nevals = 0
# Wrapper around MultiTheis solve returning pressure head at a point
function head(T::Float64,S::Float64,Q1::Float64,Qinj::Float64,rs::Array{Float64,1},t::Float64;train=false)
	Qs = [Q1,Qinj]
	if train
		global nevals = nevals + 1
	end
	#MultiTheis.solve(Qs,T,S,t,rs)
	TheisLike.solve_numerical(Qs, T, S, t, rs)
end

# This helper function doesn't work with DPFEHM. Not really needed anyhow.
#function head(T::Number,S::Number,Q1::Number;train=true)
#	# Injection rate
#	Qinj = 0.031688 # [m^3/s] (1 MMT water/yr)
#	# Time to check pressure
#	t = 30.0*86400.0 # 30 days converted to seconds
#	# Calculate distance between extraction and injection wells
#	rs = [sqrt((sum((cs-mon_well_crds[1]).^2))) for cs in well_crds]
#	#Qs = [Q1,Qinj]
#	return head(T,S,Q1,Qinj,rs,t,train=train)
#end

# Set random seed for repeatability
Random.seed!(0)

# Neural network structure
model = Chain(Dense(3, 8, σ),
              Dense(8, 16, σ),
              Dense(16, 8, σ),
              Dense(8, 1)) |> f64

# Make neural network parameters trackable by Flux
θ = params(model)

# Neural network function (T, S, target)->(Q1,Q2)
# 'target' is the target pressure head at the monitoring location
function aim(T, S, target)
  Q1 = model([T, S, target])
  # Make is so that Q1 and Q2 are always negative for extraction
  Q1 = -softplus(Q1[1])
  #Q2 = -softplus(Q2)
  Q1
end

# Set ranges for P, T, S
# Fix pressure target for now. It should be removed, but when I removed it with MultiTheis, the PIML framework didn't work as well. Moral of the story, NN's are squirrelly things that don't always make sense.
PRES = (0.0,0.0) # Pressure target range [m]
TRANS = (-2.0,0.0) # log10 transmissivity [m^2/s]
# Assume that reservoir thickness is 100m (S = S_s*b)
STOR = (-4.0,-1.0) # log10 storativity

# Random T, S, and P generator
sample() = [10^lerp(rand(), TRANS...), 10^lerp(rand(), STOR...), lerp(rand(), PRES...)]

# Physics model vs NN
function mydiff(T, S, target; train=false)
    (head(T, S, aim(T, S, target), Qinj, rs, t, train=train) - target)
end
function loss(T, S, target; train=true)
	mydiff(T, S, target, train=train)^2
end
loss(x; train=true) = sum(map(x->loss(x[1],x[2],x[3];train=train),x))
mydiff(x; train=true) = sum(map(x->mydiff(x[1],x[2],x[3];train=train),x))

function cb2()
	# Terminal output
	loss_train = sum(map(x->loss(x,train=false),data_train_batch))
	loss_test = sum(map(x->loss(x,train=false),data_test))
	println(string("callback: train loss: ", loss_train, " test loss: ", loss_test))
end

global losses_train = Float64[]
global losses_test = Float64[]
global rmses_train = Float64[]
global rmses_test = Float64[]
global meandiffs_train = Float64[]
global stddiffs_train = Float64[]
global meandiffs_test = Float64[]
global stddiffs_test = Float64[]
global train_time = Float64[]

cb = Flux.throttle(cb2, 10)

opt = ADAM()


text:
@testset "ThermalGen data misspecification" begin
    # See https://discourse.julialang.org/t/how-to-use-test-warn/15557/5 about testing for warning throwing
    warn_message = "The data doesn't include devices of type ThermalStandard, consider changing the device models"
    model = DeviceModel(ThermalStandard, PSI.ThermalStandardUnitCommitment)
    op_model = OperationModel(TestOptModel, DCPPowerModel, c_sys5_re_only)
    @test_logs (:warn, warn_message) construct_device!(op_model, :Thermal, model)
end

################################### Unit Commitment tests ##################################
@testset "Thermal UC With DC - PF" begin
    bin_variable_names = [:ON_ThermalStandard,
                          :START_ThermalStandard,
                          :STOP_ThermalStandard]
    uc_constraint_names = [:ramp_up_ThermalStandard,
                           :ramp_dn_ThermalStandard,
                           :duration_up_ThermalStandard,
                           :duration_dn_ThermalStandard]
    model = DeviceModel(PSY.ThermalStandard, PSI.ThermalStandardUnitCommitment)

    @info "5-Bus testing"
    op_model = OperationModel(TestOptModel, DCPPowerModel, c_sys5_uc)
    construct_device!(op_model, :Thermal, model)
    moi_tests(op_model, false, 480, 0, 480, 120, 120, true)
    psi_constraint_test(op_model, uc_constraint_names)
    psi_checkbinvar_test(op_model, bin_variable_names)
    psi_checkobjfun_test(op_model, GAEVF)

    op_model = OperationModel(TestOptModel, DCPPowerModel, c_sys5_uc; parameters = true)
    construct_device!(op_model, :Thermal, model)
    moi_tests(op_model, true, 480, 0, 480, 120, 120, true)
    psi_constraint_test(op_model, uc_constraint_names)
    psi_checkbinvar_test(op_model, bin_variable_names)
    psi_checkobjfun_test(op_model, GAEVF)

    @info "14-Bus testing"
    for p in [true, false]
        op_model = OperationModel(TestOptModel, DCPPowerModel, c_sys14; parameters = p)
        construct_device!(op_model, :Thermal, model)
        moi_tests(op_model, p, 480, 0, 240, 120, 120, true)
        psi_checkbinvar_test(op_model, bin_variable_names)
        psi_checkobjfun_test(op_model, GQEVF)
    end
end

@testset "Thermal UC With AC - PF" begin
    bin_variable_names = [:ON_ThermalStandard,
                          :START_ThermalStandard,
                          :STOP_ThermalStandard]
    uc_constraint_names = [:ramp_up_ThermalStandard,
                           :ramp_dn_ThermalStandard,
                           :duration_up_ThermalStandard,
                           :duration_dn_ThermalStandard]
    model = DeviceModel(PSY.ThermalStandard, PSI.ThermalStandardUnitCommitment)

    @info "5-Bus testing"
    op_model = OperationModel(TestOptModel, ACPPowerModel, c_sys5_uc)
    construct_device!(op_model, :Thermal, model)
    moi_tests(op_model, false, 600, 0, 600, 240, 120, true)
    psi_constraint_test(op_model, uc_constraint_names)
    psi_checkbinvar_test(op_model, bin_variable_names)
    psi_checkobjfun_test(op_model, GAEVF)

    op_model = OperationModel(TestOptModel, ACPPowerModel, c_sys5_uc; parameters = true)
    construct_device!(op_model, :Thermal, model)
    moi_tests(op_model, true, 600, 0, 600, 240, 120, true)
    psi_constraint_test(op_model, uc_constraint_names)
    psi_checkbinvar_test(op_model, bin_variable_names)
    psi_checkobjfun_test(op_model, GAEVF)

    @info "14-Bus testing"
    for p in [true, false]
        op_model = OperationModel(TestOptModel, ACPPowerModel, c_sys14; parameters = p)
        construct_device!(op_model, :Thermal, model)
        moi_tests(op_model, p, 600, 0, 360, 240, 120, true)
        psi_checkbinvar_test(op_model, bin_variable_names)
        psi_checkobjfun_test(op_model, GQEVF)
    end
end

################################### Basic Unit Commitment tests ############################
@testset "Thermal Basic UC With DC - PF" begin
    bin_variable_names = [:ON_ThermalStandard,
                          :START_ThermalStandard,
                          :STOP_ThermalStandard]
    model = DeviceModel(PSY.ThermalStandard, PSI.ThermalBasicUnitCommitment)

    @info "5-Bus testing"
    op_model = OperationModel(TestOptModel, DCPPowerModel, c_sys5_uc)
    construct_device!(op_model, :Thermal, model)
    moi_tests(op_model, false, 480, 0, 240, 120, 120, true)
    psi_checkbinvar_test(op_model, bin_variable_names)
    psi_checkobjfun_test(op_model, GAEVF)

    op_model = OperationModel(TestOptModel, DCPPowerModel, c_sys5_uc; parameters = true)
    construct_device!(op_model, :Thermal, model)
    moi_tests(op_model, true, 480, 0, 240, 120, 120, true)
    psi_checkbinvar_test(op_model, bin_variable_names)
    psi_checkobjfun_test(op_model, GAEVF)

    @info "14-Bus testing"
    for p in [true, false]
        op_model = OperationModel(TestOptModel, DCPPowerModel, c_sys14; parameters = p)
        construct_device!(op_model, :Thermal, model)
        moi_tests(op_model, p, 480, 0, 240, 120, 120, true)
        psi_checkbinvar_test(op_model, bin_variable_names)
        psi_checkobjfun_test(op_model, GQEVF)
    end
end

@testset "Thermal Basic UC With AC - PF" begin
    bin_variable_names = [:ON_ThermalStandard,
                          :START_ThermalStandard,
                          :STOP_ThermalStandard]
     model = DeviceModel(PSY.ThermalStandard, PSI.ThermalBasicUnitCommitment)

    @info "5-Bus testing"
    op_model = OperationModel(TestOptModel, ACPPowerModel, c_sys5_uc)
    construct_device!(op_model, :Thermal, model)
    moi_tests(op_model, false, 600, 0, 360, 240, 120, true)
    psi_checkbinvar_test(op_model, bin_variable_names)
    psi_checkobjfun_test(op_model, GAEVF)

    op_model = OperationModel(TestOptModel, ACPPowerModel, c_sys5_uc; parameters = true)
    construct_device!(op_model, :Thermal, model)
    moi_tests(op_model, true, 600, 0, 360, 240, 120, true)
    psi_checkbinvar_test(op_model, bin_variable_names)
    psi_checkobjfun_test(op_model, GAEVF)

    @info "14-Bus testing"
    for p in [true, false]
        op_model = OperationModel(TestOptModel, ACPPowerModel, c_sys14; parameters = p)
        construct_device!(op_model, :Thermal, model)
        moi_tests(op_model, p, 600, 0, 360, 240, 120, true)
        psi_checkbinvar_test(op_model, bin_variable_names)
        psi_checkobjfun_test(op_model, GQEVF)
    end
end

################################### Basic Dispatch tests ###################################
@testset "Thermal Dispatch With DC - PF" begin
    model = DeviceModel(PSY.ThermalStandard, PSI.ThermalDispatch)
    @info "5-Bus testing"
    for p in [true, false]
        op_model = OperationModel(TestOptModel, DCPPowerModel, c_sys5; parameters = p)
        construct_device!(op_model, :Thermal, model)
        moi_tests(op_model, p, 120, 120, 0, 0, 0, false)
        psi_checkobjfun_test(op_model, GAEVF)
    end

    @info "14-Bus testing"
    for p in [true, false]
        op_model = OperationModel(TestOptModel, DCPPowerModel, c_sys14; parameters = p)
        construct_device!(op_model, :Thermal, model)
        moi_tests(op_model, p, 120, 120, 0, 0, 0, false)
        psi_checkobjfun_test(op_model, GQEVF)
    end
end

@testset "Thermal Dispatch With AC - PF" begin
    model = DeviceModel(PSY.ThermalStandard, PSI.ThermalDispatch)
    @info "5-Bus testing"
    for p in [true, false]
        op_model = OperationModel(TestOptModel, ACPPowerModel, c_sys5; parameters = p)
        construct_device!(op_model, :Thermal, model)
        moi_tests(op_model, p, 240, 240, 0, 0, 0, false)
        psi_checkobjfun_test(op_model, GAEVF)
    end

    @info "14-Bus testing"
    for p in [true, false]
        op_model = OperationModel(TestOptModel, ACPPowerModel, c_sys14; parameters = p)
        construct_device!(op_model, :Thermal, model)
        moi_tests(op_model, p, 240, 240, 0, 0, 0, false)
        psi_checkobjfun_test(op_model, GQEVF)
    end
end

################################### No Minimum Dispatch tests ##############################

@testset "Thermal Dispatch NoMin With DC - PF" begin
    model = DeviceModel(PSY.ThermalStandard, PSI.ThermalDispatchNoMin)
    @info "5-Bus testing"
    for p in [true, false]
        op_model = OperationModel(TestOptModel, DCPPowerModel, c_sys5; parameters = p)
        construct_device!(op_model, :Thermal, model)
        moi_tests(op_model, p, 120, 120, 0, 0, 0, false)
        moi_lbvalue_test(op_model, :activerange_ThermalStandard, 0.0)
        psi_checkobjfun_test(op_model, GAEVF)
    end

    @info "14-Bus testing"
    for p in [true, false]
        op_model = OperationModel(TestOptModel, DCPPowerModel, c_sys14; parameters = p)
        construct_device!(op_model, :Thermal, model)
        moi_tests(op_model, p, 120, 120, 0, 0, 0, false)
        moi_lbvalue_test(op_model, :activerange_ThermalStandard, 0.0)
        psi_checkobjfun_test(op_model, GQEVF)
    end
end

@testset "Thermal Dispatch NoMin With AC - PF" begin
    model = DeviceModel(PSY.ThermalStandard, PSI.ThermalDispatchNoMin)
    @info "5-Bus testing"
    for p in [true, false]
        op_model = OperationModel(TestOptModel, ACPPowerModel, c_sys5; parameters = p)
        construct_device!(op_model, :Thermal, model)
        moi_tests(op_model, p, 240, 240, 0, 0, 0, false)
        moi_lbvalue_test(op_model, :activerange_ThermalStandard, 0.0)
        psi_checkobjfun_test(op_model, GAEVF)
    end

    @info "14-Bus testing"
    for p in [true, false]
        op_model = OperationModel(TestOptModel, ACPPowerModel, c_sys14; parameters = p)
        construct_device!(op_model, :Thermal, model)
        moi_tests(op_model, p, 240, 240, 0, 0, 0, false)
        moi_lbvalue_test(op_model, :activerange_ThermalStandard, 0.0)
        psi_checkobjfun_test(op_model, GQEVF)
    end
end

################################### Ramp Limited Testing ##################################
@testset "Thermal Ramp Limited Dispatch With DC - PF" begin
    constraint_names = [:ramp_up_ThermalStandard, :ramp_dn_ThermalStandard]
    model = DeviceModel(PSY.ThermalStandard, PSI.ThermalRampLimited)
    @info "5-Bus testing"
    for p in [true, false]
        op_model = OperationModel(TestOptModel, DCPPowerModel, c_sys5_uc; parameters = p)
        construct_device!(op_model, :Thermal, model)
        moi_tests(op_model, p, 120, 120, 96, 0, 0, false)
        psi_constraint_test(op_model, constraint_names)
        psi_checkobjfun_test(op_model, GAEVF)
    end

    @info "14-Bus testing"
    for p in [true, false]
        op_model = OperationModel(TestOptModel, DCPPowerModel, c_sys14; parameters = p)
        construct_device!(op_model, :Thermal, model)
        moi_tests(op_model, p, 120, 120, 0, 0, 0, false)
        psi_checkobjfun_test(op_model, GQEVF)
    end
end

@testset "Thermal Ramp Limited Dispatch With AC - PF" begin
    constraint_names = [:ramp_up_ThermalStandard, :ramp_dn_ThermalStandard]
    model = DeviceModel(PSY.ThermalStandard, PSI.ThermalRampLimited)
    @info "5-Bus testing"
    for p in [true, false]
        op_model = OperationModel(TestOptModel, ACPPowerModel, c_sys5_uc; parameters = p)
        construct_device!(op_model, :Thermal, model)
        moi_tests(op_model, p, 240, 240, 96, 0, 0, false)
        psi_constraint_test(op_model, constraint_names)
        psi_checkobjfun_test(op_model, GAEVF)
    end

    @info "14-Bus testing"
    for p in [true, false]
        op_model = OperationModel(TestOptModel, ACPPowerModel, c_sys14; parameters = p)
        construct_device!(op_model, :Thermal, model)
        moi_tests(op_model, p, 240, 240, 0, 0, 0, false)
        psi_checkobjfun_test(op_model, GQEVF)
    end
end

############################# UC validation tests ##########################################
# Common System Information
node = Bus(1,"nodeA", "PV", 0, 1.0, (min = 0.9, max=1.05), 230);
load = PowerLoad("Bus1", true, node,nothing, 0.4, 0.9861, 1.0, 2.0);

branches = Dict{Symbol, PSI.DeviceModel}()
services = Dict{Symbol, PSI.ServiceModel}()
ED_devices = Dict{Symbol, DeviceModel}(:Generators => PSI.DeviceModel(PSY.ThermalStandard, PSI.ThermalRampLimited),
                                        :Loads =>  PSI.DeviceModel(PSY.PowerLoad, PSI.StaticPowerLoad))
UC_devices = Dict{Symbol, DeviceModel}(:Generators => DeviceModel(PSY.ThermalStandard, PSI.ThermalStandardUnitCommitment),
                                        :Loads =>  DeviceModel(PSY.PowerLoad, PSI.StaticPowerLoad))
# Testing Ramping Constraint
@testset "Solving UC with CopperPlate for testing Ramping Constraints" begin
    DA_ramp = collect(DateTime("1/1/2024  0:00:00", "d/m/y  H:M:S"):
                        Hour(1):
                        DateTime("1/1/2024  4:00:00", "d/m/y  H:M:S"))
    gen_ramp = [ThermalStandard("Alta", true, node,0.20, 0.010,
            TechThermal(0.5,PSY.PrimeMovers(6),PSY.ThermalFuels(6),
                            (min=0.0, max=0.40), nothing,
                            nothing, nothing),
            ThreePartCost((0.0, 1400.0), 0.0, 4.0, 2.0)
            ),
            ThermalStandard("Park City", true, node,0.70,0.20,
                TechThermal(2.0,PSY.PrimeMovers(6),PSY.ThermalFuels(6),
                            (min=0.7, max=2.20), nothing,
                            (up = 0.010625, down= 0.010625),nothing),
                ThreePartCost((0.0, 1500.0), 0.0, 1.5, 0.75)
            )];
    ramp_load = [ 0.9, 1.1, 2.485, 2.175, 0.9];
    load_forecast_ramp = Deterministic(load, "scalingfactor", TimeArray(DA_ramp, ramp_load))
    ramp_test_sys = PSY.System(100.0);
    add_component!(ramp_test_sys, node)
    add_component!(ramp_test_sys, load)
    add_component!(ramp_test_sys, gen_ramp[1])
    add_component!(ramp_test_sys, gen_ramp[2])
    add_forecasts!(ramp_test_sys, [load_forecast_ramp])

    model_ref = ModelReference(CopperPlatePowerModel, ED_devices, branches, services)
    ED = OperationModel(TestOptModel, model_ref,
                        ramp_test_sys; optimizer = Cbc_optimizer,
                        parameters = true)
    psi_checksolve_test(ED, [MOI.OPTIMAL], 11191.00)
    moi_tests(ED, true, 10, 10, 10, 0, 5, false)
end

# Testing Duration Constraints
@testset "Solving UC with CopperPlate for testing Duration Constraints" begin
    DA_dur  = collect(DateTime("1/1/2024  0:00:00", "d/m/y  H:M:S"):
                        Hour(1):
                        DateTime("1/1/2024  6:00:00", "d/m/y  H:M:S"))
    gens_dur = [ThermalStandard("Alta", true, node,0.40, 0.010,
            TechThermal(0.5,PSY.PrimeMovers(6),PSY.ThermalFuels(6),
                            (min=0.3, max=0.9), nothing,
                            nothing, (up=4, down=2)),
            ThreePartCost((0.0, 1400.0), 0.0, 4.0, 2.0)
            ),
            ThermalStandard("Park City", true, node,1.70,0.20,
                TechThermal(2.2125,PSY.PrimeMovers(6),PSY.ThermalFuels(6),
                                (min=0.7, max=2.2), nothing,
                                nothing, (up=6, down=4)),
                ThreePartCost((0.0, 1500.0), 0.0, 1.5, 0.75)
            )];

    duration_load = [0.3, 0.6, 0.8, 0.7, 1.7, 0.9, 0.7]
    load_forecast_dur = Deterministic(load, "scalingfactor", TimeArray(DA_dur, duration_load))
    duration_test_sys = PSY.System(100.0);
    add_component!(duration_test_sys, node)
    add_component!(duration_test_sys, load)
    add_component!(duration_test_sys, gens_dur[1])
    add_component!(duration_test_sys, gens_dur[2])
    add_forecasts!(duration_test_sys, [load_forecast_dur]);

    status = [1.0,0.0]
    up_time = [2.0,0.0]
    down_time = [0.0,3.0]

    alta = gens_dur[1]
    init_cond = PSI.DICKDA()
    init_cond[PSI.ICKey(PSI.DeviceStatus,typeof(alta))] = build_init(gens_dur, status)
    init_cond[PSI.ICKey(PSI.TimeDurationON,typeof(alta))] = build_init(gens_dur, up_time)
    init_cond[PSI.ICKey(PSI.TimeDurationOFF,typeof(alta))] = build_init(gens_dur, down_time)

    model_ref = ModelReference(CopperPlatePowerModel, UC_devices, branches, services);
    UC = OperationModel(TestOptModel, model_ref,
                        duration_test_sys; optimizer = Cbc_optimizer,
                        parameters = true, initial_conditions = init_cond);
    psi_checksolve_test(UC, [MOI.OPTIMAL], 8223.50)
    moi_tests(UC, true, 56, 0, 56, 14, 21, true)
end

## PWL linear Cost implementation test
@testset "Solving UC with CopperPlate testing Linear PWL" begin
    gens_cost = [ThermalStandard("Alta", true, node,0.52, 0.010,
            TechThermal(0.5,PSY.PrimeMovers(6),PSY.ThermalFuels(6),
                            (min = 0.22, max = 0.55), nothing,
                            nothing, nothing),
            ThreePartCost([ (589.99, 0.220),(884.99, 0.33)
                    ,(1210.04, 0.44),(1543.44, 0.55)],532.44, 5665.23, 0.0)
            ),
            ThermalStandard("Park City", true, node,0.62,0.20,
                TechThermal(2.2125,PSY.PrimeMovers(6),PSY.ThermalFuels(6),
                            (min = 0.62, max = 1.55), nothing,
                            nothing, nothing),
                ThreePartCost([   (1264.80, 0.62),(1897.20, 0.93),
                    (2594.4787, 1.24),(3433.04, 1.55)   ], 235.397, 5665.23, 0.0)
            )];
    DA_cost  = collect(DateTime("1/1/2024  0:00:00", "d/m/y  H:M:S"):
                        Hour(1):
                        DateTime("1/1/2024  1:00:00", "d/m/y  H:M:S"))
    cost_load = [1.3,2.1];
    load_forecast_cost = Deterministic(load, "scalingfactor", TimeArray(DA_cost, cost_load))
    cost_test_sys = PSY.System(100.0);
    add_component!(cost_test_sys, node)
    add_component!(cost_test_sys, load)
    add_component!(cost_test_sys, gens_cost[1])
    add_component!(cost_test_sys, gens_cost[2])
    add_forecasts!(cost_test_sys, [load_forecast_cost])

    model_ref = ModelReference(CopperPlatePowerModel, UC_devices, branches, services);
    UC = OperationModel(TestOptModel, model_ref,
                        cost_test_sys; optimizer = Cbc_optimizer,
                        parameters = true);
    psi_checksolve_test(UC, [MOI.OPTIMAL], 9336.736919354838)
    moi_tests(UC, true, 32, 0, 8, 4, 10, true)
end

## PWL SOS-2 Cost implementation test
@testset "Solving UC with CopperPlate testing SOS2 implementation" begin
    gens_cost_sos = [ThermalStandard("Alta", true, node,0.52, 0.010,
            TechThermal(0.5,PSY.PrimeMovers(6),PSY.ThermalFuels(6),
                            (min = 0.22, max = 0.55), nothing,
                            nothing, nothing),
            ThreePartCost([ (1122.43, 0.22),(1417.43, 0.33),
                    (1742.48, 0.44),(2075.88, 0.55) ],0.0, 5665.23, 0.0)
            ),
            ThermalStandard("Park City", true, node,0.62,0.20,
                TechThermal(2.2125,PSY.PrimeMovers(6),PSY.ThermalFuels(6),
                            (min = 0.62, max = 1.55), nothing,
                            nothing, nothing),
                ThreePartCost([ (1500.19, 0.62),(2132.59, 0.929),
                    (2829.875, 1.24),(3668.444, 1.55)], 0.0, 5665.23, 0.0)
            )];
    DA_cost_sos   = collect(DateTime("1/1/2024  0:00:00", "d/m/y  H:M:S"):
                        Hour(1):
                        DateTime("1/1/2024  1:00:00", "d/m/y  H:M:S"))
    cost_sos_load = [1.3,2.1];
    load_forecast_cost_sos  = Deterministic(load, "scalingfactor", TimeArray(DA_cost_sos, cost_sos_load))
    cost_test_sos_sys = PSY.System(100.0);
    add_component!(cost_test_sos_sys, node)
    add_component!(cost_test_sos_sys, load)
    add_component!(cost_test_sos_sys, gens_cost_sos[1])
    add_component!(cost_test_sos_sys, gens_cost_sos[2])
    add_forecasts!(cost_test_sos_sys, [load_forecast_cost_sos])

    model_ref = ModelReference(CopperPlatePowerModel, UC_devices, branches, services);
    UC = OperationModel(TestOptModel, model_ref,
                        cost_test_sos_sys; optimizer = Cbc_optimizer,
                        parameters = true);
    psi_checksolve_test(UC, [MOI.OPTIMAL], 9336.736919,10.0)
    moi_tests(UC, true, 32, 0, 8, 4, 14, true)
end


text:
function recursive_factorial(n)
    if n == 0
        return 1
    else
        return n*recursive_factorial(n-1)
    end
end

recursive_factorial(7)

function recursive_factorial_info(n)
    if n == 0
        println("n == 0, returning 1")
        return 1
    else
        println("n == ", n, ", calling itself with parameter ", n-1)
        output = n*recursive_factorial_info(n-1)
        println("n == ", n, ", finished calling itself, multiplying and returning ", n, "! = ", output)
        return output
    end
end

recursive_factorial_info(7)

function ack(m,n)
    if m == 0
        return n + 1
    elseif m > 0 && n == 0
        return ack(m-1,1)
    else
        return ack(m-1, ack(m,n-1))
    end
end

ack(3,4)

function ack_info(m,n)
    function printme()
        print("ack(", m, ",", n, "): ")
    end
    if m == 0
        output = n + 1
        printme()
        println("Case 1: returning n + 1 = ", output)
        return output
    elseif m > 0 && n == 0
        printme()
        println("Case 2: calling itself with parameters m-1,1 == ", m-1, ",", 1)
        output = ack_info(m-1,1)
        printme()
        println("Case 2: finished calling itself, returning with output ", output)
        return output
    else
        printme()
        println("Case 3: calling itself for new n-value with parameters m,n-1 == ", m, ",", n-1)
        newn = ack_info(m,n-1)
        printme()
        println("Case 3: finished calling itself for new n-value == ", newn)
        printme()
        println("Case 3: calling itself with parameters m-1, A(m,n-1) == ", m-1, ",", newn)
        output = ack_info(m-1,newn)
        printme()
        println("Case 3: finished calling itself, returning ", output)
        return output
    end
end

ack_info(2,1)

function my_gcd(a,b)
    if a == 0
        return b
    elseif b == 0
        return a
    else
        return my_gcd(b, a % b)
    end
end

factor = 123_456_789
prime1 = 67_867_979
prime2 = 86_028_121
my_gcd(prime1*factor, prime2*factor)

using PyPlot

function drawTriangle(x, y, level)
    # Draw recursively colored triangles.
    # x,y are 3-vectors that define the vertices of a triangle.

    if level == 0
        # Recursion limit (depth) reached
        fill(x, y, "y") # Color whole triangle yellow
    else
        # Draw the triangle...
        plot(x[[1,2,3,1]], y[[1,2,3,1]], "k", linewidth=0.5)
        # Determine the midpoints...
        a = (x + x[[2,3,1]]) / 2
        b = (y + y[[2,3,1]]) / 2
        # Draw and color the interior triangle mauve
        fill(a, b, "m")
        # Apply the process to the three "corner" triangles...
        newx = [x a a[[3,1,2]]]
        newy = [y b b[[3,1,2]]]
        for i = 1:3
            drawTriangle(newx[i,:], newy[i,:], level - 1)
        end
    end
end

# Equilateral triangle
x = [0, 1, 0.5]
y = [0, 0, 1/sqrt(2)]

drawTriangle(x, y, 5)

function mergeLR!(L, R, x)
    # Merge the *already sorted arrays* L and R into a sorted array x
    i = j = k = 1

    # Merge L and R into x
    while i <= length(L) && j <= length(R)
        if L[i] < R[j]
            x[k] = L[i]
            i += 1
        else
            x[k] = R[j]
            j += 1
        end
        k += 1
    end

    # Copy remaining elements
    while i <= length(L)
        x[k] = L[i]
        i += 1
        k += 1
    end
    while j <= length(R)
        x[k] = R[j]
        j += 1
        k += 1
    end
end

function mergesort!(x)
    # Sort the elements of the array x using the Mergesort algorithm
    if length(x) <= 1
        return x
    else
        mid = length(x) ÷ 2   # Find the midpoint of the array
        L = x[1:mid]          # Divide array into 2 halves
        R = x[mid+1:end]

        mergesort!(L)         # Sort first half
        mergesort!(R)         # Sort second half

        mergeLR!(L, R, x)
    end
end

# Example: Sort random integers
x = rand(1:1000, 10)
println(x)
mergesort!(x)
println(x)

for array_size = Int64[1e3, 1e4, 1e5, 1e6, 1e7]
    x = rand(array_size)   # Random floating point numbers
    println("n = ", array_size)
    @time mergesort!(x)
end

function Mvalues(n)
    returned_values = Int64[]

    function M(n)
        if n > 100
            newval = n - 10
        else
            newval = M(M(n + 11))
        end
        push!(returned_values, newval)
        return newval
    end

    M(n)
    return returned_values
end

Mvalues(105)     # Easy - terminates immediately, M(105) = 95

Mvalues(97)      # More complex - finally returns M(97) = 91


text:
using QuadGK

function main()

    f(x) = x^2
    g(x) = x^4
    h(x) = (x+1)^(-1)

    a = 0
    b = 2
    n = 10

    println(simpsons(f, a, b, n))
    println(quadgk(f, a, b, rtol=1e-5))

    println(simpsons(g, a, b, n))
    println(quadgk(g, a, b, rtol=1e-5))

    println(simpsons(h, a, b, n))
    println(quadgk(h, a, b, rtol=1e-5))
end

function simpsons(df, a, b, n)

    # n should be even
    h = (b - a) / n

    xi0 = df(a) + df(b)         # end points
    xi1 = 0                     # odd points
    xi2 = 0                     # even points

    # performing integration
    for i in 1:n-1

        x = a + i*h

        # if even add to xi2
        if i%2 == 0
            xi2 += df(x)
        else
            # else odd add to xi1
            xi1 += df(x)
        end
    end

    # add everything togeteher
    xi = h*(xi0 + 2*xi2 + 4*xi1)/3

    return xi

end

main()

max_stars_repo_name: nunesmelo/djs-office-hours

max_issues_repo_name: nunesmelo/djs-office-hours

max_forks_repo_name: nunesmelo/djs-office-hours


text:
using Test
using Money: Monetary, M,
             AssetMismatch,
             convert, decimals, name

@test convert(Monetary{:EUR,18}, 10) == M(:EUR, "10")
@test convert(Monetary{:USD,18}, 10) == M(:USD, "10")
@test convert(Monetary{:EUR,18}, M(:USD,10).amount) == M(:EUR, 10)
@test_throws InexactError convert(Monetary{:EUR,18}, 10.01)

@test promote_type(Monetary{:EUR,1},Monetary{:EUR,18}) == Monetary{:EUR,18}
@test promote_type(Monetary{:EUR,18},Monetary{:EUR,1}) == Monetary{:EUR,18}
@test promote_type(Monetary{:EUR,18},String) == Monetary{:EUR,18}
@test promote_type(Monetary{:EUR,18},Integer) == Monetary{:EUR,18}

@test promote_type(Monetary{:EUR,10},Monetary{:USD,18}) == Monetary{name,18} where name

@test_throws InexactError promote_type(Monetary{:EUR,18},Real)
@test_throws InexactError promote_type(Monetary{:EUR,18},AbstractFloat)
@test_throws InexactError promote_type(Monetary{:EUR,18},Float32)
@test_throws InexactError promote_type(Monetary{:EUR,18},Float64)
@test_throws InexactError promote_type(Monetary{:EUR,18},BigFloat)

@test promote(M(:EUR, 10), 20) == (M(:EUR, 10), M(:EUR, 20))
@test promote(M(:EUR, 10), "20.02") == (M(:EUR, 10), M(:EUR, "20.02"))
@test_throws ErrorException promote(M(:EUR, 10), M(:USD, 10))
@test_throws InexactError promote(M(:EUR, 10), 20.33)
@test_throws InexactError promote(M(:EUR, 10), BigFloat(20.33))

EUR_10_57 = M(:EUR, "10.57")
EUR_0_533 = M(:EUR, "0.533")
USD_N_0_533 = M(:USD, "-0.533")

@test decimals(EUR_10_57) == 18
@test name(EUR_10_57) == :EUR

@test EUR_10_57 + "10.01" == M(:EUR, "20.58")
@test EUR_10_57 + 10 == M(:EUR, "20.57")

@test_throws InexactError EUR_10_57 + 10.01
@test_throws InexactError 10.01 + EUR_10_57
@test_throws InexactError EUR_10_57 - 10.01
@test_throws InexactError 10.01 - EUR_10_57
@test_throws InexactError EUR_10_57 * 10.01
@test_throws InexactError 10.01 * EUR_10_57
@test_throws InexactError EUR_10_57 / 10.01
@test_throws InexactError 10.01 / EUR_10_57

@test_throws ErrorException EUR_10_57 + USD_N_0_533
@test_throws ErrorException USD_N_0_533 <= EUR_10_57

@test_throws ArgumentError EUR_10_57 + "abc"

# unary operators
@test abs(USD_N_0_533).amount == EUR_0_533.amount
@test sign(EUR_0_533) == 1
@test sign(USD_N_0_533) == -1

@test round(M(:EUR, "0.40")) == M(:EUR, "0.0")
@test round(M(:EUR, "0.50")) == M(:EUR, "0.0")
@test round(M(:EUR, "0.500000000000000001")) == M(:EUR, "1.0")  # 18 decimals
@test round(M(:EUR, "0.5000000000000000001")) == M(:EUR, "0.0") # 19 decimals
@test trunc(M(:EUR, "0.6")) == M(:EUR, "0.0")
@test floor(M(:EUR, "0.6")) == M(:EUR, "0.0")
@test ceil(M(:EUR, "0.6")) == M(:EUR, "1.0")

@test isapprox(M(:EUR, 5), M(:EUR, "5.0")) == true
@test isapprox(M(:EUR, 5), M(:EUR, "5.1")) == false
@test isapprox(M(:EUR, 5), M(:EUR, "5.000000000000000001")) == false
@test isapprox(M(:EUR, 5), M(:EUR, "5.0000000000000000001")) == true

# binary operators
@test EUR_10_57 + EUR_0_533 == M(:EUR, "11.103")
@test EUR_10_57 - EUR_0_533 == M(:EUR, "10.037")
@test EUR_10_57 * M(:EUR, "2.0") == EUR_10_57 + EUR_10_57
@test -1 * EUR_10_57 == M(:EUR, "-10.57")
@test EUR_10_57 * -1 == M(:EUR, "-10.57")
@test "-1" * EUR_10_57 == M(:EUR, "-10.57")
@test EUR_10_57 * "-1" == M(:EUR, "-10.57")
@test (EUR_10_57 / M(:EUR, "0.5")) / M(:EUR, "2.0") == EUR_10_57
@test (USD_N_0_533 / M(:USD, "0.5")) / M(:USD, "2.0") == USD_N_0_533
@test (M(:EUR, "11.103") / M(:EUR, "0.5")) / M(:EUR, "2.0") == M(:EUR, "11.103")

@test EUR_10_57 == EUR_10_57
@test !(EUR_10_57 == EUR_0_533)
@test EUR_10_57 != EUR_0_533
@test EUR_10_57 > EUR_0_533
@test EUR_0_533 < EUR_10_57
@test EUR_0_533 <= EUR_10_57
@test EUR_0_533 <= EUR_0_533

@test min(EUR_10_57, EUR_0_533) == EUR_0_533
@test max(EUR_10_57, EUR_0_533) == EUR_10_57
@test minmax(EUR_10_57, EUR_0_533) == (EUR_0_533, EUR_10_57)

@test div(M(:EUR, "9"), M(:EUR, "5")) == M(:EUR, "1.0")
@test div(M(:EUR, "10"), M(:EUR, "5")) == M(:EUR, "2.0")
@test div(M(:EUR, "11"), M(:EUR, "5")) == M(:EUR, "2.0")
@test rem(M(:EUR, "7"), M(:EUR, "5")) == M(:EUR, "2.0")
@test rem(M(:EUR, "10"), M(:EUR, "5")) == M(:EUR, "0.0")
@test rem(M(:EUR, "12"), M(:EUR, "5")) == M(:EUR, "2.0")
@test divrem(M(:EUR, "11"), M(:EUR, "5")) == (M(:EUR, "2.0"),M(:EUR, "1.0"))

@test flipsign(M(:EUR, 5),   1)  == M(:EUR,  "5.0")
@test flipsign(M(:EUR, 5),  "1") == M(:EUR,  "5.0")
@test flipsign(M(:EUR, 5),  -1)  == M(:EUR, "-5.0")
@test flipsign(M(:EUR, 5), "-1") == M(:EUR, "-5.0")


text:
using TcpInstruments
using Test

scope = initialize(AgilentDSOX4034A)
@info "Successfully connected $(scope.model) at $(scope.address)"

"""
Spec:

> scope_h = initialize("192.168.1.15")

Grab data from channel 1
> data_struct = get_data(scope_h, 1)

Grab data from channel 2 and 4
> data_struct = get_data(scope_h, [2,4])

Low Pass Filter
Turn on Low Pass Filter 25
> lpf_on(scope)

Check if low pass filter is on
> get_lpf_state(scope) == "1"

Turn on Low Pass Filter 25MHz
> lpf_off(scope)
> get_lpf_state(scope) == "0"

Impedance

> set_impedance_1Mohm(scope_h)

> get_impedance(scope_h) == ONEM

> set_impedance_50ohm(scope_h)

> get_impedance(scope_h) == FIFT

Terminate TCP connection
> terminate(scope_h)

"""

data = get_data(scope, 1)

data2 = get_data(scope, [1,2])
@info typeof(data2)

@test !isempty(data.volt)
@test !isempty(data.time)

@info get_lpf_state(scope)
lpf_on(scope)
@info get_lpf_state(scope)
lpf_off(scope)
@info get_lpf_state(scope)

@info get_impedance(scope)
set_impedance_1Mohm(scope)
@info get_impedance(scope)
set_impedance_50ohm(scope)
@info get_impedance(scope)

# plot(data)

terminate(scope)
@info "Successfully disconnected"
@info "Goodbye"


text:
using Cosmology
using Test

# values from http://icosmos.co.uk/

dist_rtol = 1e-6
age_rtol = 2e-4

@testset "FlatLCDM" begin
    c = cosmology(h=0.7, OmegaM=0.3, OmegaR=0)
    @test angular_diameter_dist_mpc(c,1,rtol=dist_rtol) ≈ 1651.9145 rtol = dist_rtol
    @test comoving_radial_dist_mpc(c,1,rtol=dist_rtol) ≈ 3303.829 rtol = dist_rtol
    @test comoving_volume_gpc3(c,1,rtol=dist_rtol) ≈ 151.0571 rtol = dist_rtol
    @test luminosity_dist_mpc(c,1,rtol=dist_rtol) ≈ 6607.6579 rtol = dist_rtol
    @test distmod(c,1,rtol=dist_rtol) ≈ 44.1002 rtol = dist_rtol
    @test age_gyr(c,0,rtol=age_rtol) ≈ 13.4694 rtol = age_rtol
    @test age_gyr(c,1,rtol=age_rtol) ≈ 5.7527 rtol = age_rtol
    @test lookback_time_gyr(c,1,rtol=age_rtol) ≈ 13.4694-5.7527 rtol = age_rtol
end

@testset "OpenLCDM" begin
    c = cosmology(h=0.7, OmegaK=0.1, OmegaM=0.3, OmegaR=0)
    @test angular_diameter_dist_mpc(c,1,rtol=dist_rtol) ≈ 1619.9588 rtol = dist_rtol
    @test comoving_radial_dist_mpc(c,1,rtol=dist_rtol) ≈ 3209.784 rtol = dist_rtol
    @test comoving_volume_gpc3(c,1,rtol=dist_rtol) ≈ 140.0856 rtol = dist_rtol
    @test luminosity_dist_mpc(c,1,rtol=dist_rtol) ≈ 6479.8352 rtol = dist_rtol
    @test distmod(c,1,rtol=dist_rtol) ≈ 44.0578 rtol = dist_rtol
    @test age_gyr(c,0,rtol=age_rtol) ≈ 13.064 rtol = age_rtol
    @test age_gyr(c,1,rtol=age_rtol) ≈ 5.5466 rtol = age_rtol
    @test lookback_time_gyr(c,1,rtol=age_rtol) ≈ 13.064-5.5466 rtol = age_rtol
end

@testset "ClosedLCDM" begin
    c = cosmology(h=0.7, OmegaK=-0.1, OmegaM=0.3, OmegaR=0)
    @test angular_diameter_dist_mpc(c,1,rtol=dist_rtol) ≈ 1686.5272 rtol = dist_rtol
    @test comoving_radial_dist_mpc(c,1,rtol=dist_rtol) ≈ 3408.937 rtol = dist_rtol
    @test comoving_volume_gpc3(c,1,rtol=dist_rtol) ≈ 163.8479 rtol = dist_rtol
    @test luminosity_dist_mpc(c,1,rtol=dist_rtol) ≈ 6746.1088 rtol = dist_rtol
    @test distmod(c,1,rtol=dist_rtol) ≈ 44.1453 rtol = dist_rtol
    @test age_gyr(c,0,rtol=age_rtol) ≈ 13.925 rtol = age_rtol
    @test age_gyr(c,1,rtol=age_rtol) ≈ 5.9868 rtol = age_rtol
    @test lookback_time_gyr(c,1,rtol=age_rtol) ≈ 13.925-5.9868 rtol = age_rtol
end

@testset "FlatWCDM" begin
    c = cosmology(h=0.7, OmegaM=0.3, OmegaR=0, w0=-0.9, wa=0.1)
    @test angular_diameter_dist_mpc(c,1,rtol=dist_rtol) ≈ 1612.0585 rtol = dist_rtol
    @test comoving_radial_dist_mpc(c,1,rtol=dist_rtol) ≈ 3224.1169 rtol = dist_rtol
    @test comoving_volume_gpc3(c,1,rtol=dist_rtol) ≈ 140.3851 rtol = dist_rtol
    @test luminosity_dist_mpc(c,1,rtol=dist_rtol) ≈ 6448.2338 rtol = dist_rtol
    @test distmod(c,1,rtol=dist_rtol) ≈ 44.0472 rtol = dist_rtol
    @test age_gyr(c,0,rtol=age_rtol) ≈ 13.1915 rtol = age_rtol
    @test age_gyr(c,1,rtol=age_rtol) ≈ 5.6464 rtol = age_rtol
    @test lookback_time_gyr(c,1,rtol=age_rtol) ≈ 13.1915-5.6464 rtol = age_rtol
end

@testset "OpenWCDM" begin
    c = cosmology(h=0.7, OmegaK=0.1, OmegaM=0.3, OmegaR=0, w0=-0.9, wa=0.1)
    @test angular_diameter_dist_mpc(c,1,rtol=dist_rtol) ≈ 1588.0181 rtol = dist_rtol
    @test comoving_radial_dist_mpc(c,1,rtol=dist_rtol) ≈ 3147.6227 rtol = dist_rtol
    @test comoving_volume_gpc3(c,1,rtol=dist_rtol) ≈ 132.0466 rtol = dist_rtol
    @test luminosity_dist_mpc(c,1,rtol=dist_rtol) ≈ 6352.0723 rtol = dist_rtol
    @test distmod(c,1,rtol=dist_rtol) ≈ 44.0146 rtol = dist_rtol
    @test age_gyr(c,0,rtol=age_rtol) ≈ 12.8488 rtol = age_rtol
    @test age_gyr(c,1,rtol=age_rtol) ≈ 5.4659 rtol = age_rtol
    @test lookback_time_gyr(c,1,rtol=age_rtol) ≈ 12.8488-5.4659 rtol = age_rtol
end

@testset "ClosedWCDM" begin
    c = cosmology(h=0.7, OmegaK=-0.1, OmegaM=0.3, OmegaR=0, w0=-0.9, wa=0.1)
    @test angular_diameter_dist_mpc(c,1,rtol=dist_rtol) ≈ 1637.5993 rtol = dist_rtol
    @test comoving_radial_dist_mpc(c,1,rtol=dist_rtol) ≈ 3307.9932 rtol = dist_rtol
    @test comoving_volume_gpc3(c,1,rtol=dist_rtol) ≈ 149.8301 rtol = dist_rtol
    @test luminosity_dist_mpc(c,1,rtol=dist_rtol) ≈ 6550.3973 rtol = dist_rtol
    @test distmod(c,1,rtol=dist_rtol) ≈ 44.0813 rtol = dist_rtol
    @test age_gyr(c,0,rtol=age_rtol) ≈ 13.5702 rtol = age_rtol
    @test age_gyr(c,1,rtol=age_rtol) ≈ 5.8482 rtol = age_rtol
    @test lookback_time_gyr(c,1,rtol=age_rtol) ≈ 13.5702-5.8482 rtol = age_rtol
end

@testset "Non-Float64" begin
    # Test that FlatLCDM works with non-Float64 (BigFloat in this example)
    c = cosmology(h=0.7, OmegaM=big(0.3), OmegaR=0)
    @test angular_diameter_dist_mpc(c,1,rtol=dist_rtol) ≈ 1651.9145 rtol = dist_rtol
    # Test that FlatWCDM works with non-Float64 (BigFloat in this example)
    c = cosmology(h=big(0.7), OmegaM=0.3, OmegaR=0, w0=-0.9, wa=0.1)
    @test angular_diameter_dist_mpc(c,1,rtol=dist_rtol) ≈ 1612.0585 rtol = dist_rtol
end


text:
# Update_mixcomps.jl

export rfullcond_fstar!;

function update_μ_σ2!(mixcomp::MixComponentNormal, suffstat::Dict{Symbol, T},
    prior::PriorMixcomponent_Normal) where T <: Real

    ## mean
    v1 = 1.0 / ( 1.0 / prior.μ.σ2 + suffstat[:sumWinv] / mixcomp.σ2 )
    sd1 = sqrt(v1)
    μ1 = v1 * ( prior.μ.μ / prior.μ.σ2 + suffstat[:sumWinv] * suffstat[:yhat] / mixcomp.σ2 )

    mixcomp.μ = rand(mixcomp.rng, Distributions.Normal(μ1, sd1) )

    ## variance
    a1 = 0.5 * ( prior.σ2.ν + length(mixcomp.ζon_indx) )
    ss = suffstat[:sumWinv]*(suffstat[:yhat] - mixcomp.μ)^2 + suffstat[:s2]
    b1 = 0.5 * ( prior.σ2.ν * prior.σ2.s0 + ss )

    mixcomp.σ2 = rand(mixcomp.rng, Distributions.InverseGamma(a1, b1) )

    return nothing
end

# function rpost_f!(y::Vector{T}, μ::T, σ2::T, κ::T, Cor::Symmetric{Float64, Matrix{Float64}},
#     rng::Union{MersenneTwister, Threefry4x}) where T <: Real
#     # conjugate update for MvNormal mean with prior mean μ, known general
#     #   correlation matrix Cor, signal-to-noise ratio κ, and indep. noisy obs with variance σ2.
#
#     Corinv = inv(Cor)
#     J = PDMat_adj((Corinv ./ κ + I) ./ σ2, 1.0e-3)
#
#     h = ( y .- μ ) ./ σ2
#
#     return rand(rng, Distributions.MvNormalCanon(h, J) )
# end
# function rpost_f!(y::Vector{T}, μ::T, σ2::T, κ::T, Cor::Symmetric{Float64, Matrix{Float64}},
#     rng::Union{MersenneTwister, Threefry4x}) where T <: Real
#     # conjugate update for MvNormal mean with prior mean μ, known general
#     #   correlation matrix Cor, signal-to-noise ratio κ, and indep. noisy obs with variance σ2.
#
#     J = inv(PDMat_adj(κ*σ2*Cor)) + I/σ2
#     h = ( y .- μ ) ./ σ2
#
#     return rand(rng, Distributions.MvNormalCanon(h, J) )
# end
function rpost_f!(y::Vector{T}, μ::T, σ2::T, κ::T, Cor::Symmetric{Float64, Matrix{Float64}},
    rng::Union{MersenneTwister, Threefry4x}) where T <: Real
    # conjugate update for MvNormal mean with prior mean μ, known general
    #   correlation matrix Cor, signal-to-noise ratio κ, and indep. noisy obs with variance σ2.

    # Rasmussen & Williams '06, p.46
    K = κ * Cor
    K1 = (K + I) \ K
    Σ0 = K * (I - K1)
    Σ1 = PDMat_adj(Symmetric( σ2*Σ0 ))
    μ1 = Σ1 * ( y .- μ ) ./ σ2

    return rand(rng, Distributions.MvNormal(μ1, Σ1) )
end

function rfullcond_fstar!(f_ell::T, Celel::T, Cstst::T, Cstel::T,
    rng::Union{MersenneTwister, Threefry4x}, tol=1.0e-10) where T <: Real

    ## here C is the **Covariance** matrix

    μ = Cstel * ( Celel \ f_ell )
    Σ = Cstst - Cstel * ( Celel \ Cstel )

    if Σ < 0.0 && abs(Σ) < tol
        Σ = abs(Σ)
    end

    return sqrt(Σ) * randn(rng) + μ
end
function rfullcond_fstar!(f_ell::T, Celel::T,
    Cstst::Matrix{T}, Cstel::Vector{T},
    rng::Union{MersenneTwister, Threefry4x}, tol=1.0e-10) where T <: Real

    ## here C is the **Covariance** matrix

    μ = Cstel .* f_ell ./ Celel
    Σ = PDMat_adj(Symmetric(Cstst - Cstel * Cstel' ./ Celel), tol)

    return rand(rng, Distributions.MvNormal( μ, Σ ) )
end
function rfullcond_fstar!(f_ell::Vector{T}, Celel::Matrix{T},
    Cstst::T, Cstel::Vector{T},
    rng::Union{MersenneTwister, Threefry4x}, tol=1.0e-10) where T <: Real

    ## here C is the **Covariance** matrix

    Celel = PDMat_adj(Celel, tol)

    μ = Cstel' * ( Celel \ f_ell )
    Σ = Cstst - Cstel' * ( Celel \ Cstel )

    if Σ < 0.0 && abs(Σ) < tol
        Σ = abs(Σ)
    end

    return sqrt(Σ) * randn(rng) + μ
end
# function rfullcond_fstar!(f_ell::Vector{T}, Celel::Matrix{T},
#     Cstst::Matrix{T}, Cstel::Matrix{T},
#     rng::Union{MersenneTwister, Threefry4x}) where T <: Real
#
#     ## here C is the **Covariance** matrix
#
#     μ = Cstel * ( Celel \ f_ell )
#     Σ = PDMat_adj(Symmetric(Cstst - Cstel * ( Celel \ Cstel' )))
#
#     return rand(rng, Distributions.MvNormal( μ, Σ ) )
# end
function rfullcond_fstar!(f_ell::Vector{T}, Celel::Matrix{T},
    Cstst::Matrix{T}, Cstel::Matrix{T},
    rng::Union{MersenneTwister, Threefry4x}, tol=1.0e-10) where T <: Real

    ## here C is the **Covariance** matrix

    Celel = PDMat_adj(Celel, tol)

    μ = Cstel * ( Celel \ f_ell )
    Σ = PDMat_adj(Symmetric(Cstst - X_invA_Xt(Celel, Cstel) ), tol)

    return rand(rng, Distributions.MvNormal( μ, Σ ) )
end

# function rfullcond_fstar!_canon(fx::Vector{T}, ## already indexed
#     Cov::PDMat, on_indx::Vector{Int}, off_indx::Vector{Int},
#     rng::Union{MersenneTwister, Threefry4x}) where T <: Real
#
#     ## assumes zero mean GP
#
#     Λ = inv(Cov)
#     J = PDMat_adj( Λ.mat[off_indx, off_indx] )
#
#     η = - vec( Λ.mat[off_indx, on_indx] * fx )
#
#     return rand(rng, Distributions.MvNormalCanon( η, J ) )
# end

### worker function that operates on its own
function update_mixcomp!(mixcomp::MixComponentNormal, prior::PriorMixcomponent_Normal,
    y::Vector{T}, update::Vector{Symbol}=[:μ, :σ2, :κ, :corParams, :fx], tol_posdef=1.0e-9) where T <: Real

    nζ = length(mixcomp.ζon_indx)

    if (:κ in update) || (:corParams in update)
        suffstat = update_Cov!(mixcomp, y, prior, SNR_Hyper_ScInvChiSq(), MaternHyper_ScInvChiSq())
    elseif nζ > 0
        W_now, sumWinv_now = getW(mixcomp.κ, mixcomp.D[mixcomp.ζon_indx, mixcomp.ζon_indx], mixcomp.corParams)
        yhat_now, s2_now = getSufficients_W(y[mixcomp.ζon_indx], W_now, sumWinv_now)
        suffstat = Dict(:sumWinv => sumWinv_now, :yhat => yhat_now, :s2 => s2_now)
    end

    if nζ > 0

        if (:μ in update) || (:σ2 in update)
            update_μ_σ2!(mixcomp, suffstat, prior)
        end

        if (:fx in update)
            mixcomp.fx[mixcomp.ζon_indx] = rpost_f!(y[mixcomp.ζon_indx], mixcomp.μ, mixcomp.σ2,
                mixcomp.κ, Symmetric(mixcomp.Cor[mixcomp.ζon_indx, mixcomp.ζon_indx]), mixcomp.rng)

            ζoff_indx = setdiff( 1:length(y) , mixcomp.ζon_indx )
            Cov = PDMat_adj(mixcomp.Cor .* (mixcomp.κ * mixcomp.σ2))

            # mixcomp.fx[ζoff_indx] = rfullcond_fstar_canon(mixcomp.fx[mixcomp.ζon_indx], Cov,
            #     mixcomp.ζon_indx, ζoff_indx, mixcomp.rng)

            mixcomp.fx[ζoff_indx] = rfullcond_fstar!(mixcomp.fx[mixcomp.ζon_indx],
                Cov.mat[mixcomp.ζon_indx, mixcomp.ζon_indx], Cov.mat[ζoff_indx, ζoff_indx],
                Cov.mat[ζoff_indx, mixcomp.ζon_indx], mixcomp.rng, tol_posdef)
        end

    else

        if (:μ in update) || (:σ2 in update)
            mixcomp.μ = rand(mixcomp.rng, Normal(prior.μ.μ, sqrt(prior.μ.σ2)))
            mixcomp.σ2 = rand(mixcomp.rng, InverseGamma(0.5*prior.σ2.ν, 0.5*prior.σ2.ν*prior.σ2.s0))
        end

        if (:fx in update)
            Cov = PDMat_adj(mixcomp.Cor .* (mixcomp.κ * mixcomp.σ2), tol_posdef)
            mixcomp.fx = rand(mixcomp.rng, Distributions.MvNormal( Cov.mat ) )
        end
    end

    return mixcomp
end

### master function that calls parallel updates
function update_mixcomps!(state::State_GPMTD, prior::Prior_GPMTD, y::Vector{T},
    update::Vector{Symbol}=[:μ, :σ2, :κ, :corParams, :fx];
    n_procs::Int=1, tol_posdef=1.0e-9) where T <: Real

    L = length(state.mixcomps)

    for ℓ = 1:L
        state.mixcomps[ℓ].adapt = state.adapt
    end

    if n_procs == 1
        for ℓ = 1:L
            state.mixcomps[ℓ] = update_mixcomp!(state.mixcomps[ℓ], prior.mixcomps[ℓ], y, update, tol_posdef)
        end
    elseif n_procs > 1
        state.mixcomps = pmap(update_mixcomp!, state.mixcomps, prior.mixcomps,
            fill(y, L), fill(update, L), fill(tol_posdef, L)) # pmap will update in place on each worker, but we also need to return the mixcomps
    end

    state.accpt = [ deepcopy(state.mixcomps[ℓ].accpt) for ℓ = 1:L ]
    state.adapt_iter = deepcopy(state.mixcomps[1].adapt_iter)

    return nothing
end


text:
d = map(x -> parse.(Int, x), split.(readlines("input.txt"), "/"))

bridge = [[], 0]
result = []
function calc(b, data)
    c = [x for x in data if b[2] in x]
    if length(c) == 0
        push!(result, b)
    else
        b[1] != [] ? push!(result, b) : nothing
        for x in c
            data_left = copy(data)
            filter!(v -> v != x, data_left)
            target = b[2] == x[1] ? x[2] : x[1]
            calc([[b[1]; [x]], target], data_left)
        end
    end
end
calc(bridge, d)

# Part 1
strength(x) = sum(vcat(x[1]...))
println(maximum((map(strength, result))))
# Part 2
max_length = maximum(map(x -> length(x[1]), result))
println(maximum(map(strength, filter(x -> length(x[1]) == max_length, result))))


text:
#Demo 8: Glyph test
#-------------------------------------------------------------------------------
using InspectDR
using Colors
import Graphics: width, height

#==Input
===============================================================================#

#Constants
#-------------------------------------------------------------------------------
black = RGB24(0, 0, 0)
white = RGB24(1, 1, 1)
red = RGB24(1, 0, 0)
green = RGB24(0, 1, 0)
blue = RGB24(0, 0, 1)

#Input data
#-------------------------------------------------------------------------------
x = collect(1:10)
y0 = x.*0
_glyphs = InspectDR.SUPPORTED_GLYPHS
#_colors = Colors.colormap("Oranges", length(_glyphs))
_colors = Colors.distinguishable_colors(length(_glyphs))

#==Generate plot
===============================================================================#
mplot = InspectDR.Multiplot(title="Glyph Tests")
mplot.layout[:ncolumns] = 1

plot = add(mplot, InspectDR.Plot2D())
	strip = plot.strips[1]
	plot.xext_full = InspectDR.PExtents1D(min=0, max=11)
	strip.yext_full = InspectDR.PExtents1D(min=0, max=length(_glyphs)+1)
	plot.layout[:enable_legend] = true

	a = plot.annotation
#	a.title = "Glyph Test"

let wfrm #HIDEWARN_0.7
	for (i, g) in enumerate(_glyphs)
		wfrm = add(plot, x, y0 .+ (length(_glyphs)+1-i), id="$g")
		wfrm.line = line(style=:solid, color=_colors[i], width=2)
		wfrm.glyph = glyph(shape=g, size=10, color=black, fillcolor=_colors[i])
	end
end

gplot = display(InspectDR.GtkDisplay(), mplot)

:DONE


text:
#File treating all the prediction functions

const pred_nodes, pred_weights = (x -> (x[1] .* sqrt2, x[2] ./ sqrtπ))(gausshermite(100))

"""
    predict_f(m::AbstractGP, X_test, cov::Bool=true, diag::Bool=true)

Compute the mean of the predicted latent distribution of `f` on `X_test` for the variational GP `model`

Return also the diagonal variance if `cov=true` and the full covariance if `diag=false`
"""
predict_f

function _predict_f(
    m::GP{T}, X_test::AbstractVector; cov::Bool=true, diag::Bool=true
) where {T}
    k_star = kernelmatrix(kernel(m.f), X_test, input(m))
    μf = k_star * mean(m.f) # k * α
    if !cov
        return (μf,)
    end
    if !diag
        k_starstar = kernelmatrix(kernel(m.f), X_test) + T(jitt) * I
        covf = Symmetric(k_starstar - k_star' * (AGP.cov(m.f) \ k_star)) # k** - k* Σ⁻¹ k*
        return μf, covf
    else
        k_starstar = kernelmatrix_diag(kernel(m.f), X_test) .+ T(jitt)
        varf = k_starstar - diag_ABt(k_star / AGP.cov(m.f), k_star)
        return μf, varf
    end
end

@traitfn function _predict_f(
    m::TGP, X_test::AbstractVector; cov::Bool=true, diag::Bool=true
) where {T,TGP<:AbstractGP{T};!IsMultiOutput{TGP}}
    k_star = kernelmatrix.(kernels(m), [X_test], Zviews(m))
    μf = k_star .* (pr_covs(m) .\ means(m))
    if !cov
        return (μf,)
    end
    A = pr_covs(m) .\ (Ref(I) .- covs(m) ./ pr_covs(m))
    if !diag
        k_starstar = kernelmatrix.(kernels(m), Ref(X_test)) .+ T(jitt) * [I]
        Σf = Symmetric.(k_starstar .- k_star .* A .* transpose.(k_star))
        return μf, Σf
    else
        k_starstar =
            kernelmatrix_diag.(kernels(m), Ref(X_test)) .+
            Ref(T(jitt) * ones(T, size(X_test, 1)))
        σ²f = k_starstar .- diag_ABt.(k_star .* A, k_star)
        return μf, σ²f
    end
end

@traitfn function _predict_f(
    m::TGP, X_test::AbstractVector; cov::Bool=true, diag::Bool=true
) where {T,TGP<:AbstractGP{T};IsMultiOutput{TGP}}
    k_star = kernelmatrix(kernels(m), [X_test], Zviews(m))
    μf = k_star .* (pr_covs(m) .\ means(m))

    μf = [[sum(m.A[i][j] .* μf) for j in 1:m.nf_per_task[i]] for i in 1:nOutput(m)]
    if !cov
        return (μf,)
    end
    A = pr_covs(m) .\ ([I] .- covs(m) ./ pr_covs(m))
    if !diag
        k_starstar = (kernelmatrix.(kernels(m), [X_test]) .+ T(jitt) * [I])
        Σf = k_starstar .- k_star .* A .* transpose.(k_star)
        Σf = [[sum(m.A[i][j] .^ 2 .* Σf) for j in 1:m.nf_per_task[i]] for i in 1:nOutput(m)]
        return μf, Σf
    else
        k_starstar =
            kernelmatrix_diag.(kernels(m), [X_test]) .+ [T(jitt) * ones(T, length(X_test))]
        σ²f = k_starstar .- diag_ABt.(k_star .* A, k_star)
        σ²f = [
            [sum(m.A[i][j] .^ 2 .* σ²f) for j in 1:m.nf_per_task[i]] for i in 1:nOutput(m)
        ]
        return μf, σ²f
    end
end

function _predict_f(
    m::MCGP{T}, X_test::AbstractVector; cov::Bool=true, diag::Bool=true
) where {T}
    k_star = kernelmatrix.(kernels(m), [X_test], Zviews(m))
    f = _sample_f(m, X_test, k_star)
    μf = mean.(f)
    if !cov
        return (μf,)
    end
    if !diag
        k_starstar = kernelmatrix.(kernels(m), [X_test]) + T(jitt)
        Σf = Symmetric.(k_starstar .- invquad.(pr_covs(m), k_star) .+ StatsBase.cov.(f))
        return μf, Σf
    else
        k_starstar =
            kernelmatrix_diag.(kernels(m), [X_test]) .+ [T(jitt) * ones(T, length(X_test))]
        σ²f = k_starstar .- diag_ABt.(k_star ./ pr_covs(m), k_star) .+ StatsBase.var.(f)
        return μf, σ²f
    end
end

function _sample_f(
    m::MCGP{T,<:AbstractLikelihood,<:GibbsSampling},
    X_test::AbstractVector,
    k_star=kernelmatrix.(kernels(m), [X_test], Zviews(m)),
) where {T}
    return f = [
        Ref(k_star[k] / pr_cov(m.f[k])) .* getindex.(inference(m).sample_store, k) for
        k in 1:nLatent(m)
    ]
end

## Wrapper for vector input ##
function predict_f(
    model::AbstractGP, X_test::AbstractVector{<:Real}; cov::Bool=false, diag::Bool=true
)
    return predict_f(model, reshape(X_test, :, 1); cov=cov, diag=diag)
end

function predict_f(
    model::AbstractGP,
    X_test::AbstractMatrix;
    cov::Bool=false,
    diag::Bool=true,
    obsdim::Int=1,
)
    return predict_f(
        model, KernelFunctions.vec_of_vecs(X_test; obsdim=obsdim); cov=cov, diag=diag
    )
end

##
@traitfn function predict_f(
    model::TGP, X_test::AbstractVector; cov::Bool=false, diag::Bool=true
) where {TGP; !IsMultiOutput{TGP}}
    return first.(_predict_f(model, X_test; cov=cov, diag=diag))
end

@traitfn function predict_f(
    model::TGP, X_test::AbstractVector; cov::Bool=false, diag::Bool=true
) where {TGP; IsMultiOutput{TGP}}
    return _predict_f(model, X_test; cov=cov, diag=diag)
end

## Wrapper to predict vectors ##
function predict_y(model::AbstractGP, X_test::AbstractVector{<:Real})
    return predict_y(model, reshape(X_test, :, 1))
end

function predict_y(model::AbstractGP, X_test::AbstractMatrix; obsdim::Int=1)
    return predict_y(model, KernelFunctions.vec_of_vecs(X_test; obsdim=obsdim))
end

"""
    predict_y(model::AbstractGP, X_test::AbstractVector)
    predict_y(model::AbstractGP, X_test::AbstractMatrix; obsdim = 1)

Return
    - the predictive mean of `X_test` for regression
    - 0 or 1 of `X_test` for classification
    - the most likely class for multi-class classification
    - the expected number of events for an event likelihood
"""
@traitfn function predict_y(
    model::TGP, X_test::AbstractVector
) where {TGP <: AbstractGP; !IsMultiOutput{TGP}}
    return predict_y(likelihood(model), first(_predict_f(model, X_test; cov=false)))
end

@traitfn function predict_y(
    model::TGP, X_test::AbstractVector
) where {TGP <: AbstractGP; IsMultiOutput{TGP}}
    return predict_y.(likelihood(model), _predict_f(model, X_test; cov=false))
end

function predict_y(l::MultiClassLikelihood, μs::AbstractVector{<:AbstractVector{<:Real}})
    return [l.class_mapping[argmax([μ[i] for μ in μs])] for i in 1:length(μs[1])]
end

function predict_y(
    l::MultiClassLikelihood, μs::AbstractVector{<:AbstractVector{<:AbstractVector{<:Real}}}
)
    return predict_y(l, first(μs))
end

predict_y(l::EventLikelihood, μ::AbstractVector{<:Real}) = expec_count(l, μ)

function predict_y(l::EventLikelihood, μ::AbstractVector{<:AbstractVector})
    return expec_count(l, first(μ))
end

## Wrapper to return proba on vectors ##
function proba_y(model::AbstractGP, X_test::AbstractVector{<:Real})
    return proba_y(model, reshape(X_test, :, 1))
end

function proba_y(model::AbstractGP, X_test::AbstractMatrix; obsdim::Int=1)
    return proba_y(model, KernelFunctions.vec_of_vecs(X_test; obsdim=obsdim))
end

"""
    proba_y(model::AbstractGP, X_test::AbstractVector)
    proba_y(model::AbstractGP, X_test::AbstractMatrix; obsdim = 1)

Return the probability distribution p(y_test|model,X_test) :

    - Tuple of vectors of mean and variance for regression
    - Vector of probabilities of y_test = 1 for binary classification
    - Dataframe with columns and probability per class for multi-class classification
"""

@traitfn function proba_y(
    model::TGP, X_test::AbstractVector
) where {TGP <: AbstractGP; !IsMultiOutput{TGP}}
    μ_f, Σ_f = _predict_f(model, X_test; cov=true)
    return pred = compute_proba(model.likelihood, μ_f, Σ_f)
end

@traitfn function proba_y(
    model::TGP, X_test::AbstractVector
) where {TGP <: AbstractGP; IsMultiOutput{TGP}}
    return proba_multi_y(model, X_test)
end

function proba_multi_y(model::AbstractGP, X_test::AbstractVector)
    μ_f, Σ_f = _predict_f(model, X_test; cov=true)
    return preds = compute_proba.(likelihood(model), μ_f, Σ_f)
end

function compute_proba(
    l::AbstractLikelihood,
    μ::AbstractVector{<:AbstractVector},
    σ²::AbstractVector{<:AbstractVector},
)
    return compute_proba(l, first(μ), first(σ²))
end

function StatsBase.mean_and_var(lik::AbstractLikelihood, fs::AbstractMatrix)
    vals = lik.(eachcol(fs))
    return StatsBase.mean(vals), StatsBase.var(vals)
end

function proba_y(model::MCGP, X_test::AbstractVector{<:Real}; nSamples::Int=100)
    return proba_y(model, reshape(X_test, :, 1); nSamples)
end

function proba_y(model::MCGP, X_test::AbstractMatrix; nSamples::Int=100, obsdim=1)
    return proba_y(model, KernelFunctions.vec_of_vecs(X_test; obsdim); nSamples)
end

function proba_y(model::MCGP, X_test::AbstractVector; nSamples::Int=200)
    nTest = length(X_test)
    f = first(_sample_f(model, X_test))
    return proba, sig_proba = mean_and_var(compute_proba_f.(likelihood(model), f))
end

function compute_proba_f(l::AbstractLikelihood, f::AbstractVector{<:Real})
    return compute_proba.(l, f)
end

function compute_proba(
    l::AbstractLikelihood, ::AbstractVector, ::AbstractVector
) where {T<:Real}
    return error("Non implemented for the likelihood $l")
end

for f in [:predict_f, :predict_y, :proba_y]
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

# To propose an update to the Cbc binaries please follow these steps:
# 1) Fork this CbcBuilder repository and create a new branch
# 2) Modify the build script (e.g., changing the portion marked by START/END-EASY-CHANGE-BLOCK) and test that it compiles
#    for all architectures by running `julia --color=yes build_tarballs.jl --verbose`. (Note that this build make take a
#    long time. Also, adding the --debug argument to the command will drop you into the build environment for debugging)
# 3) Create a github release in your fork of CbcBuilder with the generated tarballs and build_CbcBuilder.vX.X.X.jl
#    (https://github.com/tcnksm/ghr is usefull to automate this release creation)
# 4) Fork the Cbc.jl repository (https://github.com/JuliaOpt/Cbc.jl) and create a new branch
# 5) Update the portion marked by ## START-VERSION-UPDATE-BLOCK / ## END-VERSION-UPDATE-BLOCK in Cbc.jl/deps/build.jl
#    on your branch of Cbc.jl with the corresponding lines in the build_CbcBuilder.vX.X.X.jl file you generated in steps 2)/3)
# 6) Check that tests pass on your branch of Cbc (e.g. `]test Cbc`)
# 7) Create a PR for CbcBuilder.jl (No need to create a PR for Cbc.jl)

# Ideally, any update would need changes only on the EASY-CHANGE-BLOCK below
# Note that if two sources have the same version then they need to have
# different extension so that the corresponding filenames are different.
# For instance, consider the following case where COINMumps and COINLapack
# both use versio 1.6.0:
#    COINMumps_version = v"1.6.0"
#    COINMumps_extension = "tar.gz"
#    COINLapack_version = v"1.6.0"
#    COINLapack_extension = "zip"
# In this case, we need one of the libraries to use the "tar.gz" extension and
# the other to use the "zip" extension. The reason for this is that the file names
# for github releases for COINMumps and COINLapack for the configuration above will
# be respectively `1.6.0.tar.gz` and `1.6.0.zip`. If we instead had used the same
# extensions, one of the files would have been overwritten.

##START-EASY-CHANGE-BLOCK
Cbc_version = v"2.10.3"
Cbc_extension = "tar.gz"
Cbc_hash = "944c1f717e653aa04f4d6c775aa91d749ab5990d6ee9c3e4f76b9df5bc6e2a8d"
Cgl_version = v"0.60.2"
Cgl_extension = "tar.gz"
Cgl_hash = "c965b7173ab0e1a6c6d4be785813311e7869ce5333b952f70837309abd82e659"
Clp_version = v"1.17.3"
Clp_extension = "tar.gz"
Clp_hash = "25f0692fe1daa492e7801770af6991506ae9a8c34a4cae358d017400a02dfcf8"
Osi_version = v"0.108.5"
Osi_extension = "tar.gz"
Osi_hash = "c9a6f098e2824883bb3ec1f12df5987b7a8da0f1241988a5dd4663ac362e6381"
CoinUtils_version = v"2.11.3"
CoinUtils_extension = "tar.gz"
CoinUtils_hash = "7c4753816e765974941db75ec89f8855e56b86959f3a5f068fdf95b0003be61c"
COINMumps_version = v"1.6.0"
COINMumps_extension = "tar.gz"
COINMumps_hash = "3f2bb7d13333e85a29cd2dadc78a38bbf469bc3920c4c0933a90b7d8b8dc798a"
COINMetis_version = v"1.3.5"
COINMetis_extension = "tar.gz"
COINMetis_hash = "98a6110d5d004a16ad42ee26cfac508477f44aa6fe296b90a6413fe0273ebe24"
COINLapack_version = v"1.6.0"
COINLapack_extension = "zip"
COINLapack_hash = "227969f240176c8e1f391548f8f854bf81ac13c9c3f9803b345eaa052a399b3a"
COINBLAS_version = v"1.4.6"
COINBLAS_extension = "tar.gz"
COINBLAS_hash = "f9601efb98f04fdba220d49d5bda98d2a5a5e2ed7564df339bc7149b0c303f0c"
ASL_version = v"3.1.0"
ASL_extension = "tar.gz"
ASL_hash = "587c1a88f4c8f57bef95b58a8586956145417c8039f59b1758365ccc5a309ae9"
##END-EASY-CHANGE-BLOCK

name = "CbcBuilder"
# Collection of sources required to build CbcBuilder
sources = [
    "./bundled",
    "https://github.com/coin-or/Cbc/archive/releases/$(Cbc_version).$(Cbc_extension)" =>
    Cbc_hash,
    "https://github.com/coin-or/Cgl/archive/releases/$(Cgl_version).$(Cgl_extension)" =>
    Cgl_hash,
    "https://github.com/coin-or/Clp/archive/releases/$(Clp_version).$(Clp_extension)" =>
    Clp_hash,
    "https://github.com/coin-or/Osi/archive/releases/$(Osi_version).$(Osi_extension)" =>
    Osi_hash,
    "https://github.com/coin-or/CoinUtils/archive/releases/$(CoinUtils_version).$(CoinUtils_extension)" =>
    CoinUtils_hash,
    "https://github.com/coin-or-tools/ThirdParty-Mumps/archive/releases/$(COINMumps_version).$(COINMumps_extension)" =>
    COINMumps_hash,
    "https://github.com/coin-or-tools/ThirdParty-Metis/archive/releases/$(COINMetis_version).$(COINMetis_extension)" =>
    COINMetis_hash,
    "https://github.com/coin-or-tools/ThirdParty-Lapack/archive/releases/$(COINLapack_version).$(COINLapack_extension)" =>
    COINLapack_hash,
    "https://github.com/coin-or-tools/ThirdParty-Blas/archive/releases/$(COINBLAS_version).$(COINBLAS_extension)"  =>
    COINBLAS_hash,
    "https://github.com/ampl/mp/archive/$(ASL_version).tar.gz"  =>
    ASL_hash,
]

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

### Preliminaries

# Standard BB updates
update_configure_scripts

# Fix some paths
for path in ${LD_LIBRARY_PATH//:/ }; do
    for file in $(ls $path/*.la); do
        echo "$file"
        baddir=$(sed -n "s|libdir=||p" $file)
        sed -i~ -e "s|$baddir|'$path'|g" $file
    done
done

# Osi, Cgl, Clp, CoinUtils, mumps, metis, lapack, blas, asl
if [ $target = "x86_64-apple-darwin14" ]; then
  export AR=/opt/x86_64-apple-darwin14/bin/x86_64-apple-darwin14-ar
fi

### Compile ASL

# Use staticfloat's cross-compile trick for ASL https://github.com/ampl/mp/issues/115
cd $WORKSPACE/srcdir/mp-*
rm -rf thirdparty/benchmark
patch -p1 < $WORKSPACE/srcdir/asl-extra/no_benchmark.patch
# Build ASL
mkdir build
cd build
if [ $target = "x86_64-w64-mingw32" ] || [ $target = "i686-w64-mingw32" ]; then
   cmake -DCMAKE_C_FLAGS='-fPIC -DPIC' -DBUILD_SHARED_LIBS=OFF -DCMAKE_INSTALL_PREFIX=$prefix -DCMAKE_TOOLCHAIN_FILE=/opt/$target/$target.toolchain -DRUN_HAVE_STD_REGEX=0 -DRUN_HAVE_STEADY_CLOCK=0 -DHAVE_ACCESS_DRIVER_EXITCODE=0 -DHAVE_EXCEL_DRIVER_EXITCODE=0 -DHAVE_ODBC_TEXT_DRIVER_EXITCODE=0    ../
else
   cmake -DCMAKE_C_FLAGS='-fPIC -DPIC' -DBUILD_SHARED_LIBS=OFF -DCMAKE_INSTALL_PREFIX=$prefix  -DCMAKE_TOOLCHAIN_FILE=/opt/$target/$target.toolchain       -DRUN_HAVE_STD_REGEX=0       -DRUN_HAVE_STEADY_CLOCK=0       ../
fi
# Copy over pregenerated files after building arithchk, so as to fake out cmake,
# because cmake will delete our arith.h
## If this fails it is ok, we just want to prevend cmake deleting arith.h
set +e
make arith-h VERBOSE=1
set -e
mkdir -p src/asl
cp -v $WORKSPACE/srcdir/asl-extra/expr-info.cc ../src/expr-info.cc
cp -v $WORKSPACE/srcdir/asl-extra/arith.h.${target} src/asl/arith.h
# Build and install ASL
make -j${nproc} VERBOSE=1
make install VERBOSE=1

### Compile Lapack
cd $WORKSPACE/srcdir
cd ThirdParty-Blas-releases-*/
./get.Blas
mkdir build
cd build/

# configure, make and install
../configure --prefix=$prefix --with-pic --disable-pkg-config --host=${target} --disable-shared --enable-static --enable-dependency-linking lt_cv_deplibs_check_method=pass_all
make -j${nproc}
make install

### Compile Lapack
cd $WORKSPACE/srcdir
cd ThirdParty-Lapack-releases-*/
./get.Lapack
mkdir build
cd build/

# configure, make and install
../configure --prefix=$prefix --with-pic --disable-pkg-config  --host=${target} --disable-shared --enable-static --enable-dependency-linking lt_cv_deplibs_check_method=pass_all \
--with-blas-lib="-L$prefix/lib -lcoinblas"
make -j${nproc}
make install

### Compile Metis
cd $WORKSPACE/srcdir
cd ThirdParty-Metis-releases-*/
./get.Metis
mkdir build
cd build/

# configure, make and install
../configure --prefix=$prefix --with-pic --disable-pkg-config --host=${target} --disable-shared --enable-static --enable-dependency-linking lt_cv_deplibs_check_method=pass_all
make -j${nproc}
make install

### Compile Mumps
cd $WORKSPACE/srcdir
cd ThirdParty-Mumps-releases-*/
./get.Mumps
patch -p1 < $WORKSPACE/srcdir/mumps-extra/quiet.diff
mkdir build
cd build/

# configure, make and install
../configure --prefix=$prefix --with-pic --disable-pkg-config --host=${target} --disable-shared --enable-static --enable-dependency-linking lt_cv_deplibs_check_method=pass_all \
--with-metis-lib="-L${prefix}/lib -lcoinmetis" --with-metis-incdir="$prefix/include/coin/ThirdParty" \
--with-blas-lib="-L$prefix/lib -lcoinblas"
make -j${nproc}
make install

### Compile CoinUtils
cd $WORKSPACE/srcdir
cd CoinUtils-releases-*/
mkdir build
cd build/

# Set env vars: CPPFLAGS
if [ $target = "aarch64-linux-gnu" ] || [ $target = "arm-linux-gnueabihf" ]; then
   export CPPFLAGS="-std=c++11 -D_GLIBCXX_USE_CXX11_ABI=1"
fi

# configure, make and install
../configure --prefix=$prefix --with-pic --disable-pkg-config  --host=${target} --disable-shared --enable-static --enable-dependency-linking lt_cv_deplibs_check_method=pass_all \
--with-blas-lib="-L${prefix}/lib -lcoinblas" \
--with-lapack-lib="-L${prefix}/lib -lcoinlapack"
make -j${nproc}
make install

# Clear env vars: CPPFLAGS
export CPPFLAGS=""

### Compile Osi

cd $WORKSPACE/srcdir
cd Osi-releases-*/
mkdir build
cd build/

# Set env vars: CXXFLAGS and CPPFLAGS
export CXXFLAGS="-std=c++11"
if [ $target = "aarch64-linux-gnu" ] || [ $target = "arm-linux-gnueabihf" ]; then
   export CPPFLAGS="-std=c++11"
fi

# configure, make and install
../configure --prefix=$prefix --with-pic --disable-pkg-config --host=${target} --disable-shared --enable-static --enable-dependency-linking lt_cv_deplibs_check_method=pass_all \
--with-coinutils-lib="-L${prefix}/lib -lCoinUtils" --with-coinutils-incdir="$prefix/include/coin" \
--with-lapack-lib="-L${prefix}/lib -lcoinlapack" \
--with-blas-lib="-L${prefix}/lib -lcoinblas"
make -j${nproc}
make install

# Clear env vars: CXXFLAGS and CPPFLAGS
export CXXFLAGS=""
export CPPFLAGS=""

### Compile Cgl and Clp

## Compile Clp

cd $WORKSPACE/srcdir
cd Clp-releases-*/
mkdir build
cd build/

# Set env vars: CPPFLAGS
if [ $target = "aarch64-linux-gnu" ]; then
   export CPPFLAGS="-DNDEBUG -w -DCOIN_USE_MUMPS_MPI_H -D__arm__ -std=c++11 -D_GLIBCXX_USE_CXX11_ABI=1"
elif [ $target = "arm-linux-gnueabihf" ]; then
   export CPPFLAGS="-DNDEBUG -w -DCOIN_USE_MUMPS_MPI_H -D__arm__ -mfpu=neon -std=c++11 -D_GLIBCXX_USE_CXX11_ABI=1"
else
   export CPPFLAGS="-DNDEBUG -w -DCOIN_USE_MUMPS_MPI_H"
fi

# configure, make and install
if [ $target = "x86_64-w64-mingw32" ] || [ $target = "i686-w64-mingw32" ]; then
  ../configure --prefix=$prefix --with-pic --disable-pkg-config --host=${target} --disable-shared --enable-static \
  --enable-dependency-linking lt_cv_deplibs_check_method=pass_all \
  --with-asl-lib="-L${prefix}/lib -lasl" --with-asl-incdir="$prefix/include/asl" \
  --with-blas="-L${prefix}/lib -lcoinblas" \
  --with-lapack="-L${prefix}/lib -lcoinlapack" \
  --with-metis-lib="-L${prefix}/lib -lcoinmetis" --with-metis-incdir="$prefix/include/coin/ThirdParty" \
  --with-mumps-lib="-L${prefix}/lib -lcoinmumps" --with-mumps-incdir="$prefix/include/coin/ThirdParty" \
  --with-coinutils-lib="-L${prefix}/lib -lCoinUtils" --with-coinutils-incdir="$prefix/include/coin" \
  --with-osi-lib="-L${prefix}/lib -lOsi" --with-osi-incdir="$prefix/include/coin"
else
  ../configure --prefix=$prefix --with-pic --disable-pkg-config --host=${target} --disable-shared --enable-static \
  --enable-dependency-linking lt_cv_deplibs_check_method=pass_all \
  --with-asl-lib="-L${prefix}/lib -lasl" --with-asl-incdir="$prefix/include/asl" \
  --with-blas="-L${prefix}/lib -lcoinblas" \
  --with-lapack="-L${prefix}/lib -lcoinlapack" \
  --with-metis-lib="-L${prefix}/lib -lcoinmetis" --with-metis-incdir="$prefix/include/coin/ThirdParty" \
  --with-mumps-lib="-L${prefix}/lib -lcoinmumps" --with-mumps-incdir="$prefix/include/coin/ThirdParty" \
  --with-coinutils-lib="-L${prefix}/lib -lCoinUtils" --with-coinutils-incdir="$prefix/include/coin" \
  --with-osi-lib="-L${prefix}/lib -lOsi" --with-osi-incdir="$prefix/include/coin" \
  LDFLAGS=-ldl;
fi
make -j${nproc}
make install

# Clear env vars: CPPFLAGS
export CPPFLAGS=""

## Compile Cgl

cd $WORKSPACE/srcdir
cd Cgl-releases-*/
mkdir build
cd build/

# Set env vars: CPPFLAGS
if [ $target = "aarch64-linux-gnu" ]; then
   export CPPFLAGS="-DNDEBUG -w -DCOIN_USE_MUMPS_MPI_H -D__arm__ -std=c++11 -D_GLIBCXX_USE_CXX11_ABI=1"
elif [ $target = "arm-linux-gnueabihf" ]; then
   export CPPFLAGS="-DNDEBUG -w -DCOIN_USE_MUMPS_MPI_H -D__arm__ -mfpu=neon -std=c++11 -D_GLIBCXX_USE_CXX11_ABI=1"
else
   export CPPFLAGS="-DNDEBUG -w -DCOIN_USE_MUMPS_MPI_H"
fi

# configure, make and install
if [ $target = "x86_64-w64-mingw32" ] || [ $target = "i686-w64-mingw32" ]; then
  ../configure --prefix=$prefix --disable-pkg-config --with-pic --host=${target} --disable-shared --enable-static \
  --enable-dependency-linking lt_cv_deplibs_check_method=pass_all \
  --with-coinutils-lib="-L${prefix}/lib -lCoinUtils" --with-coinutils-incdir="$prefix/include/coin" \
  --with-osi-lib="-L${prefix}/lib -lOsi" --with-osi-incdir="$prefix/include/coin" \
  --with-osiclp-lib="-L${prefix}/lib -lOsiClp" --with-osiclp-incdir="$prefix/include/coin"
else
  ../configure --prefix=$prefix --disable-pkg-config --with-pic --host=${target} --disable-shared --enable-static \
  --enable-dependency-linking lt_cv_deplibs_check_method=pass_all \
  --with-coinutils-lib="-L${prefix}/lib -lCoinUtils" --with-coinutils-incdir="$prefix/include/coin" \
  --with-osi-lib="-L${prefix}/lib -lOsi" --with-osi-incdir="$prefix/include/coin" \
  --with-osiclp-lib="-L${prefix}/lib -lOsiClp" --with-osiclp-incdir="$prefix/include/coin" \
  LDFLAGS=-ldl;
fi
make -j${nproc}
make install

# Clear env vars: CPPFLAGS
export CPPFLAGS=""

### Compile Cbc
cd $WORKSPACE/srcdir
cd Cbc-releases-*/
mkdir build
cd build/

# Apply patch related to https://github.com/JuliaOpt/Cbc.jl/issues/117 and https://github.com/coin-or/Cbc/issues/267
cd ../Cbc/src
mv Cbc_C_Interface.cpp Cbc_C_Interface.cpp.orig
patch -l Cbc_C_Interface.cpp.orig $WORKSPACE/srcdir/cbc-extra/no_lp.patch  -o Cbc_C_Interface.cpp
cd ../../build

# force only exporting symbols related to Clp
sed -i~ -e 's|LT_LDFLAGS="-no-undefined"|LT_LDFLAGS="-no-undefined -export-symbols-regex \\"Cbc\\""|g' ../configure
sed -i~ -e 's|LT_LDFLAGS="-no-undefined"|LT_LDFLAGS="-no-undefined -export-symbols-regex \\"Cbc\\""|g' ../Cbc/configure

# configure, make and install
if [ $target = "x86_64-apple-darwin14" ]; then
    # For mac we need more to hide symbols
    sed -i~ -e "s|~nmedit -s \$output_objdir/\${libname}-symbols.expsym \${lib}| -exported_symbols_list \$output_objdir/\${libname}-symbols.expsym|g" ../configure

    ../configure --prefix=$prefix --with-pic --disable-pkg-config --host=${target} --enable-shared --disable-static \
    --enable-dependency-linking lt_cv_deplibs_check_method=pass_all \
    --with-cgl-lib="-L${prefix}/lib -lCgl" --with-cgl-incdir="$prefix/include/coin" \
    --with-asl-lib="-L${prefix}/lib -lasl" --with-asl-incdir="$prefix/include/asl" \
    --with-blas-lib="-L${prefix}/lib -lcoinblas" \
    --with-lapack-lib="-L${prefix}/lib -lcoinlapack" \
    --with-metis-lib="-L${prefix}/lib -lcoinmetis" --with-metis-incdir="$prefix/include/coin/ThirdParty" \
    --with-mumps-lib="-L${prefix}/lib -lcoinmumps" --with-mumps-incdir="$prefix/include/coin/ThirdParty" \
    --with-coinutils-lib="-L${prefix}/lib -lCoinUtils" --with-coinutils-incdir="$prefix/include/coin" \
    --with-osi-lib="-L${prefix}/lib -lOsi" --with-osi-incdir="$prefix/include/coin" \
    --with-clp-lib="${prefix}/lib/libClp.a ${prefix}/lib/libOsiClp.a" --with-clp-incdir="$prefix/include/coin" \
    --with-coindepend-lib="-L${prefix}/lib -lCgl -lOsi ${prefix}/lib/libClp.a -lCoinUtils -lcoinmumps -lcoinmetis -lbz2 -lz -lcoinlapack -lcoinblas" \
    --enable-cbc-parallel \
    LDFLAGS=-ldl;
elif [ $target = "x86_64-w64-mingw32" ] || [ $target = "i686-w64-mingw32" ]; then
    ../configure --prefix=$prefix --with-pic --disable-pkg-config --host=${target} --enable-shared --disable-static \
    --enable-dependency-linking lt_cv_deplibs_check_method=pass_all \
    --with-cgl-lib="-L${prefix}/lib -lCgl" --with-cgl-incdir="$prefix/include/coin" \
    --with-asl-lib="-L${prefix}/lib -lasl" --with-asl-incdir="$prefix/include/asl" \
    --with-blas-lib="-L${prefix}/lib -lcoinblas" \
    --with-lapack-lib="-L${prefix}/lib -lcoinlapack" \
    --with-metis-lib="-L${prefix}/lib -lcoinmetis" --with-metis-incdir="$prefix/include/coin/ThirdParty" \
    --with-mumps-lib="-L${prefix}/lib -lcoinmumps" --with-mumps-incdir="$prefix/include/coin/ThirdParty" \
    --with-coinutils-lib="-L${prefix}/lib -lCoinUtils" --with-coinutils-incdir="$prefix/include/coin" \
    --with-osi-lib="-L${prefix}/lib -lOsi" --with-osi-incdir="$prefix/include/coin" \
    --with-clp-lib="${prefix}/lib/libClp.a ${prefix}/lib/libOsiClp.a" --with-clp-incdir="$prefix/include/coin" \
    --with-coindepend-lib="-L${prefix}/lib -lCgl -lOsi ${prefix}/lib/libClp.a -lcoinmumps -lcoinmetis -lCoinUtils -lcoinlapack -lcoinblas -lgfortran" \
    --enable-cbc-parallel
else
    ../configure --prefix=$prefix --with-pic --disable-pkg-config --host=${target} --enable-shared --disable-static \
    --with-asl-lib="${prefix}/lib/libasl.a" --with-asl-incdir="$prefix/include/asl" \
    --with-blas-lib="${prefix}/lib/libcoinblas.a -lgfortran" \
    --with-lapack-lib="${prefix}/lib/libcoinlapack.a" \
    --with-metis-lib="${prefix}/lib/libcoinmetis.a" --with-metis-incdir="$prefix/include/coin/ThirdParty" \
    --with-mumps-lib="${prefix}/lib/libcoinmumps.a -lgfortran" --with-mumps-incdir="$prefix/include/coin/ThirdParty" \
    --with-coinutils-lib="${prefix}/lib/libcoinlapack.a ${prefix}/lib/libcoinblas.a ${prefix}/lib/libCoinUtils.a -lbz2 -lz" --with-coinutils-incdir="$prefix/include/coin" \
    --with-osi-lib="${prefix}/lib/libOsi.a" --with-osi-incdir="$prefix/include/coin" \
    --with-cgl-lib="${prefix}/lib/libCgl.a" --with-cgl-incdir="$prefix/include/coin" \
    --with-clp-lib="${prefix}/lib/libClp.a ${prefix}/lib/libOsiClp.a" --with-clp-incdir="$prefix/include/coin" \
    --with-coindepend-lib="${prefix}/lib/libCgl.a ${prefix}/lib/libOsi.a ${prefix}/lib/libClp.a ${prefix}/lib/libcoinmumps.a ${prefix}/lib/libcoinmetis.a ${prefix}/lib/libCoinUtils.a ${prefix}/lib/libcoinlapack.a ${prefix}/lib/libcoinblas.a -lgfortran" \
    lt_cv_deplibs_check_method=pass_all \
    --enable-cbc-parallel \
    LDFLAGS=-ldl;
fi
make -j${nproc}

# Clean-up bin directory before installing
rm ${prefix}/bin/*

# Install
make install

# Clean-up lib directory
rm ${prefix}/lib/*.a
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Linux(:i686, libc=:glibc),
    Linux(:x86_64, libc=:glibc),
    Linux(:aarch64, libc=:glibc),
    Linux(:armv7l, libc=:glibc, call_abi=:eabihf),
    MacOS(:x86_64),
    Windows(:i686),
    Windows(:x86_64)
]
platforms = expand_gcc_versions(platforms)
# To fix gcc4 bug in Windows
platforms = setdiff(platforms, [Windows(:x86_64, compiler_abi=CompilerABI(:gcc4)), Windows(:i686, compiler_abi=CompilerABI(:gcc4))])
push!(platforms, Windows(:i686,compiler_abi=CompilerABI(:gcc6)))
push!(platforms, Windows(:x86_64,compiler_abi=CompilerABI(:gcc6)))

# It seems Clp in ARM now needs C++11, so gcc4 is not working.
# A possible fix is to switch to gcc6 as for windows, but this will require updating all dependencies below too.
platforms = setdiff(platforms, [Linux(:armv7l, libc=:glibc, call_abi=:eabihf, compiler_abi=CompilerABI(:gcc4)), Linux(:aarch64, libc=:glibc, compiler_abi=CompilerABI(:gcc4))])

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

# Dependencies that must be installed before this package can be built
dependencies = [
#    "https://github.com/JuliaOpt/CglBuilder/releases/download/v0.60.2-static/build_CglBuilder.v0.60.2.jl",
  #  "https://github.com/JuliaOpt/ClpBuilder/releases/download/v1.16.11-1-static/build_ClpBuilder.v1.16.11.jl",
#    "https://github.com/JuliaOpt/OsiBuilder/releases/download/v0.108.5-static/build_OsiBuilder.v0.108.5.jl",
   # "https://github.com/JuliaOpt/CoinUtilsBuilder/releases/download/v2.11.3-static/build_CoinUtilsBuilder.v2.11.3.jl",
#    "https://github.com/JuliaOpt/COINMumpsBuilder/releases/download/v1.6.0-1-static/build_COINMumpsBuilder.v1.6.0.jl",
#    "https://github.com/JuliaOpt/COINMetisBuilder/releases/download/v1.3.5-1-static/build_COINMetisBuilder.v1.3.5.jl",
#    "https://github.com/JuliaOpt/COINLapackBuilder/releases/download/v1.6.0-static/build_COINLapackBuilder.v1.6.0.jl",
#    "https://github.com/JuliaOpt/COINBLASBuilder/releases/download/v1.4.6-1-static/build_COINBLASBuilder.v1.4.6.jl",
#    "https://github.com/JuliaOpt/ASLBuilder/releases/download/v3.1.0-1-static/build_ASLBuilder.v3.1.0.jl"
]

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


text:
using KernelFunctions
using Zygote, ForwardDiff, Tracker
using Test, LinearAlgebra

dims = [10,5]

A = rand(dims...)
B = rand(dims...)
K = [zeros(dims[1],dims[1]),zeros(dims[2],dims[2])]
kernels = [SqExponentialKernel,MaternKernel,Matern32Kernel,Matern52Kernel]
l = 2.0
vl = l*ones(dims[1])
testfunction(k,A,B) = det(kernelmatrix(k,A,B))
testfunction(k,A) = det(kernelmatrix(k,A))

##TODO Eventually store real results in file
@testset "Zygote Automatic Differentiation test" begin
    @testset "ARD" begin
        for k in kernels
            @testset "$k" begin
                @test all(isapprox.(Zygote.gradient(x->testfunction(k(x),A,B),vl)[1], ForwardDiff.gradient(x->testfunction(k(x),A,B),vl)))
                @test  all(isapprox.(Zygote.gradient(x->testfunction(k(vl),x,B),A)[1],ForwardDiff.gradient(x->testfunction(k(vl),x,B),A)))
                @test all(isapprox.(Zygote.gradient(x->testfunction(k(x),A),vl)[1],ForwardDiff.gradient(x->testfunction(k(x),A),vl)))
                @test all(isapprox.(Zygote.gradient(x->testfunction(k(vl),x),A)[1],ForwardDiff.gradient(x->testfunction(k(vl),x),A)))
            end
        end
    end
    @testset "ISO" begin
        for k in kernels
            @testset "$k" begin
                @test all(isapprox.(Zygote.gradient(x->testfunction(k(x),A,B),l)[1],ForwardDiff.gradient(x->testfunction(k(x[1]),A,B),[l])[1]))
                @test all(isapprox.(Zygote.gradient(x->testfunction(k(l),x,B),A)[1],ForwardDiff.gradient(x->testfunction(k(l),x,B),A)))
                @test all(isapprox.(Zygote.gradient(x->testfunction(k(x),A),l)[1],ForwardDiff.gradient(x->testfunction(k(x[1]),A),[l])))
                @test all(isapprox.(Zygote.gradient(x->testfunction(k(l),x),A)[1],ForwardDiff.gradient(x->testfunction(k(l[1]),x),A)))
            end
        end
    end
end

@testset "ForwardDiff AutomaticDifferentation test" begin
    @testset "ARD" begin
        for k in kernels
            @test_nowarn ForwardDiff.gradient(x->testfunction(k(x),A,B),vl)
            @test_nowarn ForwardDiff.gradient(x->testfunction(k(vl),x,B),A)
            @test_nowarn ForwardDiff.gradient(x->testfunction(k(x),A),vl)
            @test_nowarn ForwardDiff.gradient(x->testfunction(k(vl),x),A)
        end
    end
    @testset "ISO" begin
        for k in kernels
            @test_nowarn ForwardDiff.gradient(x->testfunction(k(x[1]),A,B),[l])
            @test_nowarn ForwardDiff.gradient(x->testfunction(k(l),x,B),A)
            @test_nowarn ForwardDiff.gradient(x->testfunction(k(x[1]),A),[l])
            @test_nowarn ForwardDiff.gradient(x->testfunction(k(l[1]),x),A)
        end
    end
end

@testset "Tracker AutomaticDifferentation test" begin
    @testset "ARD" begin
        for k in kernels
            @test_broken all(Tracker.gradient(x->testfunction(k(x),A,B),vl)[1] .≈ ForwardDiff.gradient(x->testfunction(k(x),A,B),vl))
            @test_broken all(Tracker.gradient(x->testfunction(k(vl),x,B),A)[1] .≈ ForwardDiff.gradient(x->testfunction(k(vl),x,B),A))
            @test_broken all(Tracker.gradient(x->testfunction(k(x),A),vl)[1] .≈  ForwardDiff.gradient(x->testfunction(k(x),A),vl))
            @test_broken all.(Tracker.gradient(x->testfunction(k(vl),x),A) .≈ ForwardDiff.gradient(x->testfunction(k(vl),x),A))
        end
    end
    @testset "ISO" begin
        for k in kernels
            @test_broken Tracker.gradient(x->testfunction(k(x[1]),A,B),[l])
            @test_broken Tracker.gradient(x->testfunction(k(l),x,B),A)
            @test_broken Tracker.gradient(x->testfunction(k(x[1]),A),[l])
            @test_broken Tracker.gradient(x->testfunction(k(l),x),A)

        end
    end
end


text:
function main(r::Robot,n::Int)
    x,y=go_to_left_bottom_corner_and_get_x_y(r)
    go_and_get_num_steps(r,Nord)
    all=go_and_get_num_steps(r,Sud)+1
    side=Nord
    t=0
    k=n
    kk=2n
    while isborder(r,Ost)==false
        if t>=k && t<kk
            for i in 1:n
                if isborder(r,side)==false
                    move!(r,side)
                end
            end
        end
        put_markers_in_direction(r,side,n,all)
        move!(r,Ost)
        side=inverse(side)
        go_and_get_num_steps(r,side)
        side=inverse(side)
        t=t+1
        if t>kk
            k=k+2n
            kk=kk+2n
        end
    end
    if t>=k && t<kk
        for i in 1:n
            if isborder(r,side)==false
                move!(r,side)
            end
        end
    end
    put_markers_in_direction(r,side,n,all)
    go_to_left_bottom_corner_and_get_x_y(r)
    go_back(r,x,y)
end

function put_markers_in_direction(r::Robot,side::HorizonSide,n::Int,all::Int)
    k=0
    ii=0
    while isborder(r,side)==false
        ii=0
        for i in 1:n
            if isborder(r,side)==false
                putmarker!(r)
                move!(r,side)
                k=k+1
                ii=ii+1
            end
        end
        if isborder(r,side)==true && ii!=n
            putmarker!(r)
        end
        for j in 1:n
            if isborder(r,side)==false
                move!(r,side)
                k=k+1
            end
        end
    end
    mm=2n+1
    while all>=mm
        if all==mm
            move!(r,HorizonSide(mod(Int(side)+2,4)))
            if ismarker(r)==false
                move!(r,side)
                putmarker!(r)
            end
        end
        mm=mm+3n
    end

    return k
end

function go_to_left_bottom_corner_and_get_x_y(r::Robot)
    x=0
    y=0
    while isborder(r,West)==false
        move!(r,West)
        x=x+1
    end
    while isborder(r,Sud)==false
        move!(r,Sud)
        y=y+1
    end
    return (x,y)
end

function go_and_get_num_steps(r::Robot,side::HorizonSide) :: Int
    n=0
    while isborder(r,side)==false
        move!(r,side)
        n=n+1
    end
    return n
end

function go_back(r,x,y)
    for i in 1:x
        move!(r,Ost)
    end
    for j in 1:y
        move!(r,Nord)
    end
end

inverse(side::HorizonSide) = HorizonSide(mod(Int(side)+2, 4))

max_stars_repo_name: VladislavBukin/mirea-progs

max_issues_repo_name: VladislavBukin/mirea-progs

max_forks_repo_name: VladislavBukin/mirea-progs