text:
module Form

using Gtk, Gtk.ShortNames, Gtk.GConstants
using ..Ahorn

abstract type Option

end

struct Section
    name::String
    dataName::String

    options::Array{Option, 1}
    fieldOrder::Array{String, 1}

    function Section(name::String, options::Array{Option, 1}=Option[]; dataName::String=name, fieldOrder::Array{String, 1}=String[])
        return new(name, dataName, options, fieldOrder)
    end
end

# ListOption Helpers

function columnEditCallback(store::GtkListStore, col::Integer, row::Integer, value::Any)
    if typeof(store[row, col]) <: Number
        if Ahorn.isNumber(value)
            store[row, col] = Ahorn.parseNumber(value)
        end

    else
        store[row, col] = string(value)
    end
end

function setListButtonSensitivity(option::Option)
    rows, cols = Base.size(option.container.store)

    GAccessor.sensitive(option.addButton, rows < option.maxRows || option.maxRows == -1)
    GAccessor.sensitive(option.removeButton, rows > option.minRows)
end

function addListRow(option::Option)
    rows, cols = Base.size(option.container.store)

    if rows < option.maxRows || option.maxRows == -1
        push!(option.container.store, tuple([typ == String ? "0" : zero(typ) for typ in option.container.dataType.parameters]...))
        Ahorn.selectRow!(option.container, rows + 1)

        setListButtonSensitivity(option)
    end
end

function deleteListRow(option::Option)
    rows, cols = Base.size(option.container.store)

    if hasselection(option.container.selection)
        if rows > option.minRows || option.minRows == -1
            row = selected(option.container.selection)
            deleteat!(option.container.store, row)

            setListButtonSensitivity(option)
        end
    end
end

# Entry Constructors and functions

struct TextEntryOption <: Option
    name::String
    dataName::String

    label::Gtk.GtkLabel
    entry::Gtk.GtkEntry

    as::Type

    function TextEntryOption(name::String; dataName::String=name, tooltip::String="", value::String="", as::Type=typeof(value))
        label = Label(name, xalign=0.0, margin_start=8, tooltip_text=tooltip)
        entry = Entry(text=string(value))

        return new(name, dataName, label, entry, as)
    end
end

Base.size(option::TextEntryOption) = (2, 1)
getValue(option::TextEntryOption) = Ahorn.getEntryText(option.entry, option.as)
setValue!(option::TextEntryOption, value::Any) = Ahorn.setEntryText!(option.entry, string(convert(option.as, value)))
getGroup(option::TextEntryOption) = 0
setGtkProperty!(option::TextEntryOption, field::Symbol, value::Any) = set_gtk_property!(option.entry, field, value)
getGtkProperty(option::TextEntryOption, field::Symbol, typ::DataType) = get_gtk_property(option.entry, field, typ)

function addToGrid!(grid::Gtk.GtkGrid, option::TextEntryOption, col::Integer=0, row::Integer=0)
    grid[col, row] = option.label
    grid[col + 1, row] = option.entry
end

struct TextViewOption <: Option
    name::String
    dataName::String

    textView::Gtk.GtkTextView
    scrollable::Gtk.GtkScrolledWindow

    function TextViewOption(name::String; dataName::String=name, tooltip::String="", value::String="", monospace::Bool=true, vexpand::Bool=true, hexpand::Bool=true)
        textView = TextView(
            vexpand=vexpand,
            hexpand=hexpand,
            monospace=monospace,
            tooltip_text=tooltip
        )

        Ahorn.setTextViewText!(textView, value)

        scrollable = ScrolledWindow(vexpand=vexpand, hexpand=hexpand)
        push!(scrollable, textView)

        return new(name, dataName, textView, scrollable)
    end
end

Base.size(option::TextViewOption) = (4, 2)
getValue(option::TextViewOption) = Ahorn.getTextViewText(option.textView)
setValue!(option::TextViewOption, value::String) = Ahorn.setTextViewText!(option.textView, value)
getGroup(option::TextViewOption) = 4
setGtkProperty!(option::TextViewOption, field::Symbol, value::Any) = set_gtk_property!(option.textView, field, value)
getGtkProperty(option::TextViewOption, field::Symbol, typ::DataType) = set_gtk_property!(option.textView, field, typ)

function addToGrid!(grid::Gtk.GtkGrid, option::TextViewOption, col::Integer=0, row::Integer=0)
    grid[col:col + 3, row:row + 1] = option.scrollable
end

mutable struct TextChoiceOption <: Option
    name::String
    dataName::String

    label::Gtk.GtkLabel
    combobox::Gtk.GtkComboBoxText

    value::String
    choices::Array{String, 1}

    as::Type

    function TextChoiceOption(name::String, options::Array{T, 1}; dataName::String=name, tooltip::String="", value::H="", editable::Bool=false, sortChoices::Bool=false) where {T, H}
        label = Label(name, xalign=0.0, margin_start=8, tooltip_text=tooltip)
        combobox = ComboBoxText(editable)

        choices = string.(options)

        if sortChoices
            sort!(choices)
        end

        append!(combobox, choices)
        Ahorn.setComboIndex!(combobox, options, value)

        return new(name, dataName, label, combobox, string(value), choices, typeof(value))
    end
end

Base.size(option::TextChoiceOption) = (2, 1)
getValue(option::TextChoiceOption) = typeof(option.value) === String ? Ahorn.convertString(option.as, option.value) : option.value
setValue!(option::TextChoiceOption, value::Any) = Ahorn.setComboIndex!(option.combobox, option.choices, string(value))
getGroup(option::TextChoiceOption) = 1
setGtkProperty!(option::TextChoiceOption, field::Symbol, value::Any) = set_gtk_property!(option.combobox, field, value)
getGtkProperty(option::TextChoiceOption, field::Symbol, typ::DataType) = set_gtk_property!(option.combobox, field, typ)

function addToGrid!(grid::Gtk.GtkGrid, option::TextChoiceOption, col::Integer=0, row::Integer=0)
    grid[col, row] = option.label
    grid[col + 1, row] = option.combobox

    @guarded signal_connect(option.combobox, "changed") do args...
        option.value = Gtk.GLib.@sigatom Gtk.bytestring(GAccessor.active_text(option.combobox))
    end
end

mutable struct DictionaryChoiceOption <: Option
    name::String
    dataName::String

    options::Dict{String, Any}

    label::Gtk.GtkLabel
    combobox::Gtk.GtkComboBoxText

    value::Any
    as::Type

    function DictionaryChoiceOption(name::String, options::Dict{String, T}; dataName::String=name, tooltip::String="", value::H=nothing, editable::Bool=false, sortChoices::Bool=false) where {T, H}
        label = Label(name, xalign=0.0, margin_start=8, tooltip_text=tooltip)
        combobox = ComboBoxText(editable)

        choices = collect(options)

        if sortChoices
            sort!(choices, by=p -> p[1])
        end

        keys = getindex.(choices, 1)
        values = getindex.(choices, 2)

        append!(combobox, keys)

        index = Ahorn.setComboIndex!(combobox, values, value)

        return new(name, dataName, options, label, combobox, value, typeof(value))
    end
end

Base.size(option::DictionaryChoiceOption) = (2, 1)
getValue(option::DictionaryChoiceOption) = option.value
setValue!(option::DictionaryChoiceOption, value::Any) = Ahorn.setComboIndex!(option.combobox, collect(values(option.options)), value, allowCustom=false)
getGroup(option::DictionaryChoiceOption) = 1
setGtkProperty!(option::DictionaryChoiceOption, field::Symbol, value::Any) = set_gtk_property!(option.combobox, field, value)
getGtkProperty(option::DictionaryChoiceOption, field::Symbol, typ::DataType) = set_gtk_property!(option.combobox, field, typ)

function addToGrid!(grid::Gtk.GtkGrid, option::DictionaryChoiceOption, col::Integer=0, row::Integer=0)
    grid[col, row] = option.label
    grid[col + 1, row] = option.combobox

    @guarded signal_connect(option.combobox, "changed") do args...
        text = Gtk.bytestring(GAccessor.active_text(option.combobox))

        if haskey(option.options, text)
            option.value = get(option.options, text, text)

        else
            try
                option.value = Ahorn.convertString(option.as, text)

            catch
                option.value = first(option.options)[2]
            end
        end
    end
end

struct CheckBoxOption <: Option
    name::String
    dataName::String

    checkbox::Gtk.GtkCheckButton

    function CheckBoxOption(name::String; dataName::String=name, tooltip::String="", value::Bool=false)
        checkbox = CheckButton(name, active=value, tooltip_text=tooltip)

        return new(name, dataName, checkbox)
    end
end

Base.size(option::CheckBoxOption) = (1, 1)
getValue(option::CheckBoxOption) = GAccessor.active(option.checkbox)
setValue!(option::CheckBoxOption, value::Bool) = GAccessor.active(option.checkbox, value)
getGroup(option::CheckBoxOption) = 2
setGtkProperty!(option::CheckBoxOption, field::Symbol, value::Any) = set_gtk_property!(option.checkbox, field, value)
getGtkProperty(option::CheckBoxOption, field::Symbol, typ::DataType) = set_gtk_property!(option.checkbox, field, typ)

function addToGrid!(grid::Gtk.GtkGrid, option::CheckBoxOption, col::Integer=0, row::Integer=0)
    grid[col, row] = option.checkbox
end

struct SpinOption <: Option
    name::String
    dataName::String

    label::Gtk.GtkLabel
    spinButton::Gtk.GtkSpinButton

    function SpinOption(name::String, range::UnitRange{T}; dataName::String=name, tooltip::String="", value::T=first(range)) where T
        label = Label(name, xalign=0.0, margin_start=8, tooltip_text=tooltip)
        spinButton = SpinButton(range)

        GAccessor.value(spinButton, value)

        return new(name, dataName, label, spinButton)
    end
end

Base.size(option::SpinOption) = (1, 1)
getValue(option::SpinOption) = GAccessor.value(option.spinButton)
setValue!(option::SpinOption, value::Integer) = GAccessor.value(option.spinButton, value)
getGroup(option::SpinOption) = 0
setGtkProperty!(option::SpinOption, field::Symbol, value::Any) = set_gtk_property!(option.spinButton, field, value)
getGtkProperty(option::SpinOption, field::Symbol, typ::DataType) = set_gtk_property!(option.spinButton, field, typ)

function addToGrid!(grid::Gtk.GtkGrid, option::SpinOption, col::Integer=0, row::Integer=0)
    grid[col, row] = option.label
    grid[col + 1, row] = option.spinButton
end

# TODO - Handle columns kwarg
struct ListOption <: Option
    name::String
    dataName::String

    container::Ahorn.ListContainer
    scrollable::Gtk.GtkScrolledWindow

    addButton::Gtk.GtkButton
    removeButton::Gtk.GtkButton

    editable::Bool

    minRows::Integer
    maxRows::Integer

    columns::Integer

    function ListOption(name::String, data::Array{T, 1}; dataName::String=name, tooltip::String="", editable::Bool=true, minRows::Integer=-1, maxRows::Integer=-1, columns::Integer=4) where T <: Tuple
        @Ahorn.catchall begin
            headers = tuple(string.(split(name, ';'))...)
            container = Ahorn.generateTreeView(headers, data, sortable=false, editable=fill(editable, length(headers)), callbacks=Array{Ahorn.listViewCallbackUnion, 1}(fill(columnEditCallback, length(headers))))
            scrollableWindow = ScrolledWindow(vexpand=true, hscrollbar_policy=Gtk.GtkPolicyType.NEVER)
            push!(scrollableWindow, container.tree)

            addButton = Button("Add row")
            removeButton = Button("Delete row")

            return new(name, dataName, container, scrollableWindow, addButton, removeButton, editable, minRows, maxRows, columns)
        end
    end
end

Base.size(option::ListOption) = (4, option.editable ? 2 : 1)
getValue(option::ListOption) = Ahorn.getListData(option.container)
setValue!(option::ListOption, value::Array{T, 1}) where T <: Tuple = Ahorn.updateTreeView(option.container, value)
getGroup(option::ListOption) = 4
setGtkProperty!(option::ListOption, field::Symbol, value::Any) = set_gtk_property!(option.tree, field, value)
getGtkProperty(option::ListOption, field::Symbol, typ::DataType) = set_gtk_property!(option.tree, field, typ)

function addToGrid!(grid::Gtk.GtkGrid, option::ListOption, col::Integer=0, row::Integer=0)
    rows, cols = Base.size(option.container.store)

    grid[col:col + 4, row] = option.scrollable

    if rows > 0
        Ahorn.selectRow!(option.container, 1)
    end

    if option.editable
        grid[col:col + 1, row + 1] = option.addButton
        grid[col + 2:col + 4, row + 1] = option.removeButton

        signal_connect(w -> addListRow(option), option.addButton, "clicked")
        signal_connect(w -> deleteListRow(option), option.removeButton, "clicked")
    end

    setListButtonSensitivity(option)
end

# Section Functions

function groupOptions(options::Array{Option, 1})
    groups = Dict{Integer, Array{Option, 1}}()

    for option in options
        group = getGroup(option)

        if !haskey(groups, group)
            groups[group] = Option[]
        end

        push!(groups[group], option)
    end

    return groups
end

function getOptionByDataName(options::Array{Option, 1}, name::String)
    index = findfirst(isequal(name), getfield.(options, :dataName))

    return index === nothing ? nothing : options[index]
end

function getOptionsData(options::Array{Option, 1})
    data = Dict{String, Any}()
    incorrectOptions = Option[]

    for option in options
        try
            data[option.dataName] = getValue(option)

        catch e
            push!(incorrectOptions, option)
        end
    end

    if !isempty(incorrectOptions)
        return nothing, incorrectOptions
    end

    return data, incorrectOptions
end

function setOptionsData!(options::Array{Option, 1}, data::Dict{String, Any})
    for option in options
        if haskey(data, option.dataName)
            setValue!(option, data[option.dataName])
        end
    end
end

function getSectionData(section::Section)
    return getOptionsData(section.options)
end

function setSectionData!(section::Section, data::Dict{String, Any})
    return setOptionsData!(section.options, data)
end

function getSectionsData(sections::Array{Section, 1}, packIfSingleSection::Bool=false)
    if length(sections) == 1 && !packIfSingleSection
        return getSectionData(sections[1])
    end

    sectionsData = Dict{String, Dict{String, Any}}()
    incorrectSectionsOptions = Dict{Section, Array{Option, 1}}()
    success = true

    for section in sections
        data, incorrectOptions = getSectionData(section)

        if data !== nothing
            sectionsData[section.dataName] = data

        else
            incorrectSectionsOptions[section] = incorrectOptions
            success = false
        end
    end

    if !success
        return nothing, incorrectSectionsOptions
    end

    return sectionsData, incorrectSectionsOptions
end

function setSectionsData!(sections::Array{Section, 1}, data::Dict{String, Any})
    if length(sections) == 1
        if length(data) == 1 && isa(first(data)[2], Dict)
            setSectionData!(sections[1], first(data)[2])

        else
            setSectionData!(sections[1], data)
        end
    end

    data = Dict{String, Dict{String, Any}}()

    for section in sections
        setSectionData!(section, data[section.dataName])
    end
end

function addOption!(grid::Gtk.GtkGrid, option::Option, col::Integer, row::Integer, columns::Integer=4)
    needsCols, needsRows = size(option)

    if col + needsCols > columns
        row += 1
        col = 0
    end

    addToGrid!(grid, option, col, row)

    col += needsCols
    row += needsRows - 1

    return col, row
end

function generateSectionGrid(section::Section; columns::Integer=4, separateGroups::Bool=isempty(section.fieldOrder))
    grid = Grid()

    col = 0
    row = 0

    orderedOptions = filter((option) -> option.dataName in section.fieldOrder, section.options)
    sort!(orderedOptions, by=(option) -> findfirst(isequal(option.dataName), section.fieldOrder))

    for option in orderedOptions
        col, row = addOption!(grid, option, col, row, columns)
    end

    groups = groupOptions(section.options)

    for (group, options) in sort(collect(groups), by=g -> g[1])
        # Adding in front as this also covers orderedOptions, Grid won't display this as empty row
        if separateGroups
            row += 1
            col = 0
        end

        for option in sort(options, by=o -> o.name)
            # Already handled
            if option.dataName in section.fieldOrder
                continue
            end

            col, row = addOption!(grid, option, col, row, columns)
        end
    end

    return grid
end

function generateSectionsNotebook(sections::Array{Section, 1}; columns::Integer=4, separateGroups::Bool=true, gridIfSingleSection::Bool=true)
    useNotebook = length(sections) > 1
    notebook = Notebook()

    if gridIfSingleSection && !useNotebook
        return generateSectionGrid(sections[1], columns=columns, separateGroups=separateGroups)
    end

    for section in sections
        sectionGrid = generateSectionGrid(section, columns=columns, separateGroups=separateGroups)
        sectionBox = Box(:v)

        push!(sectionBox, sectionGrid)
        push!(notebook, sectionBox, section.name)
    end

    return notebook
end

function getIncorrectOptionsMessage(incorrectOptions::Dict{Section, Array{Option, 1}})
    lines = String[]
    flattenedOptions = Option[]

    for (section, options) in incorrectOptions
        append!(lines, [
            length(options) == 1 ? "The following field is invalid in section '$(section.name)':" : "The following fields are invalid in section '$(section.name)':"
            join(String[option.name for option in options], ", ")
        ])
    end

    return join(lines, "\n")
end

function getIncorrectOptionsMessage(incorrectOptions::Array{Option, 1})
    lines = String[
        length(incorrectOptions) == 1 ? "The following field is invalid:" : "The following fields are invalid:"
        join(String[option.name for option in incorrectOptions], ", ")
    ]

    return join(lines, "\n")
end

# TODO - Respect validation function
function createFormWindow(title::String, sections::Union{Array{Section, 1}, Section}; columns::Integer=4, separateGroups::Bool=true, gridIfSingleSection::Bool=true, buttonText::String="Update", callback::Function=(data) -> println(data), parent::Gtk.GtkWindow=Ahorn.window, icon::Pixbuf=Ahorn.windowIcon, canResize::Bool=false)
    sections = isa(sections, Section) ? Section[sections] : sections
    updateButton = Button(buttonText)

    content = generateSectionsNotebook(sections, columns=columns, separateGroups=separateGroups, gridIfSingleSection=gridIfSingleSection)

    window = Window(title, -1, -1, canResize, icon=icon) |> (Frame() |> (box = Box(:v)))

    push!(box, content)
    push!(box, updateButton)

    @guarded signal_connect(updateButton, "clicked") do args...
        data, incorrectOptions = getSectionsData(sections)

        if data !== nothing
            callback(data)

        else
            Ahorn.topMostInfoDialog(getIncorrectOptionsMessage(incorrectOptions), window)
        end
    end

    return window
end

function createFormWindow(title::String, options::Array{Option, 1}; columns::Integer=4, separateGroups::Bool=true, gridIfSingleSection::Bool=true, buttonText::String="Update", callback::Function=(data) -> println(data), parent::Gtk.GtkWindow=Ahorn.window)
    return createFormWindow(title, Section("section", options), columns=columns, fieldOrder=fieldOrder, separateGroups=separateGroups, gridIfSingleSection=gridIfSingleSection, buttonText=buttonText, callback=callback, parent=parent)
end

function suggestOption(displayName::String, value::Any; tooltip::String="", dataName::String=displayName, choices::Union{Array, Dict, Nothing}=nothing, editable::Bool=false, sortChoices::Bool=false, preferredType::Union{Type, Nothing}=nothing)
    if isa(choices, Array)
        return TextChoiceOption(displayName, choices, dataName=dataName, tooltip=tooltip, value=value, editable=editable, sortChoices=sortChoices)

    elseif isa(choices, Dict)
        return DictionaryChoiceOption(displayName, choices, dataName=dataName, tooltip=tooltip, value=value, editable=editable, sortChoices=sortChoices)

    elseif isa(value, Bool)
        return CheckBoxOption(displayName, value=value, dataName=dataName, tooltip=tooltip)

    elseif isa(value, Char) || isa(value, String) || isa(value, Number)
        fieldType = something(preferredType, typeof(value))

        # Default to Int and Float32 for abstract number types
        # Celeste doesn't use any floats besides Float32
        # Never any Float besides Float32, causes visual "decimal rounding" issues
        if fieldType <: Integer
            fieldType = Int

        elseif fieldType <: Number
            fieldType = Float32
        end

        return TextEntryOption(displayName, value=string(convert(fieldType, value)), dataName=dataName, as=fieldType, tooltip=tooltip)

    elseif isa(value, Array) && !isempty(displayName)
        return ListOption(displayName, value, dataName=dataName, tooltip=tooltip, editable=editable)
    end
end

end

max_stars_repo_name: CelestialCartographers/Ahorn

max_issues_repo_name: CelestialCartographers/Ahorn

max_forks_repo_name: CelestialCartographers/Ahorn


text:
###
# UDY Oracle
###
# At the moment only supports cutting planes
import JuMP.UnsetSolver

using JuMP, Mosek  #Right now SDPs require Mosek

export UDYOracle

type UDYOracle <: AbstractOracle
    eps_::Float64
    gamma1::Float64
    gamma2::Float64
    muhat::Vector{Float64}
    covhat::Array{Float64, 2}
    C::Array{Float64, 2}  #C'C = inv(covhat)

    # Cutting plane algorithm
    cut_tol::Float64  ##defaults to 1e-6
    cut_model::Model
    cut_vars::Vector{Variable}

    # Other options
    debug_printcut::Bool
end

#preferred interface
function UDYOracle(data, eps_, alpha; numBoots=10000, cut_tol=1e-6, debug_printcut=false)
    muhat  = vec(mean(data, 1))
    covhat = cov(data)
    C = chol(inv(covhat))
    gamma1 = bootDY_mu(data, alpha/2, numBoots)
    gamma2 = bootDY_sigma(data, alpha/2, numBoots)
    UDYOracle(eps_, gamma1, gamma2, muhat, covhat, C,
                cut_tol, Model(solver=MosekSolver()), Variable[], debug_printcut)
end

UDYOracle(muhat, covhat, gamma1, gamma2, eps_; cut_tol=1e-6, debug_printcut=false) =
    UDYOracle(eps_, gamma1, gamma2, muhat, covhat, chol(inv(covhat)),
                cut_tol, Model(solver=MosekSolver()), Variable[], debug_printcut)

# JuMPeR alerting us that we're handling this contraint
registerConstraint(w::UDYOracle, rm::Model, ind::Int, prefs) =
    ! get(prefs, :prefer_cuts, true) && error("Only cutting plane supported")

#VG Make setup fcns consistent in terms of checks...
function setup(wDY::UDYOracle, rm::Model, prefs)
    wDY.cut_model = Model(solver=MosekSolver())

    # Extract preferences we care about
    wDY.debug_printcut = get(prefs, :debug_printcut, false)
    wDY.cut_tol        = get(prefs, :cut_tol, wDY.cut_tol)

    rd = JuMPeR.getRobust(rm)
    d = size(wDY.covhat, 1)
    @assert (rd.numUncs == d) "Num Uncertainties doesn't match columns in data"
    @assert (length(rd.uncertaintyset) == 0) #does not support additional cnsts on unctertainties for now

    #Assumes bounds were set directly on the robust contraints
    #build the cut solver
    @defMatrixVar(wDY.cut_model, u[d])
    @defMatrixVar(wDY.cut_model, w[d])
    @defMatrixVar(wDY.cut_model, m[d])
    @defVar(wDY.cut_model, 0 <= lam <=1/wDY.eps_)
    @defSDPVar(wDY.cut_model, A[d, d])
    @defSDPVar(wDY.cut_model, Ahat[d, d])

    addConstraint(wDY.cut_model, norm(wDY.C * w) <= sqrt(wDY.gamma1)*lam)

    for i = 1:d
        @addConstraint(wDY.cut_model, rd.uncLower[i] <= u[i] <= rd.uncUpper[i])
        @addConstraint(wDY.cut_model, rd.uncLower[i]*(lam-1) <= m[i] )
        @addConstraint(wDY.cut_model, m[i] <= rd.uncUpper[i]*(lam-1))

        @addConstraint(wDY.cut_model,
            wDY.muhat[i]*lam  == u[i] + m[i] + w[i])
    end

    Z1 = [lam-1  m';
          m      A]
    addConstraint(wDY.cut_model, Z1 >= 0)

    Z2 = [1. u';
         u  Ahat ]
    addConstraint(wDY.cut_model, Z2 >= 0)

    #add an auxiliary variable for the big matrix expression
    @defMatrixVar(wDY.cut_model, Z3[d, d])
    sigbar = (wDY.gamma2 * wDY.covhat + wDY.muhat * wDY.muhat')
    for i = 1:d
        for j = 1:d
            @addConstraint(wDY.cut_model, Z3[i, j] == sigbar[i,j] * lam)
        end
    end

    temp = wDY.muhat * w'
    addConstraint(wDY.cut_model, Z3-A-Ahat-temp-temp' >=0)

    #VG Duplicate to make other things easier
    @defVar(wDY.cut_model, u2[1:d])
    for i = 1:d
        addConstraint(wDY.cut_model, u2[i] == u[i])
    end
    wDY.cut_vars = u2[:]
end

function generateCut(w::UDYOracle, m::Model, rm::Model, inds::Vector{Int}, active=false)
    rd = JuMPeR.getRobust(rm)
    master_sol = m.colVal
    new_cons = {}

    for ind in inds
        # Update the cutting plane problem's objective
        con = JuMPeR.get_uncertain_constraint(rm, ind)
        #VG This should move to sparse notation
        cut_sense, unc_obj_coeffs, lhs_const =
                JuMPeR.build_cut_objective(rm, con, master_sol)

        setObjective(w.cut_model, cut_sense, sum([unc_obj_coeffs[ix] * w.cut_vars[ix] for ix=1:length(w.cut_vars)]))
        cut_solve_status = solve(w.cut_model, suppress_warnings=true)
        cut_solve_status != :Optimal && error("Cutting plane problem failed: $cut_solve_status")
        ustar = getValue(w.cut_vars)
        zstar = getObjectiveValue(w.cut_model)

        lhs_of_cut = zstar + lhs_const

        # SUBJECT TO CHANGE: active cut detection
        if active
            push!(rd.activecuts[ind],
                JuMPeR.cut_to_scen(ustar,
                    JuMPeR.check_cut_status(con, lhs_of_cut, w.cut_tol) == :Active))
            continue
        end

        # Check violation
        if JuMPeR.check_cut_status(con, lhs_of_cut, w.cut_tol) != :Violate
            w.debug_printcut && JuMPeR.debug_printcut(rm ,m,w,lhs_of_cut,con,nothing)
            continue  # No violation, no new cut
        end

        # Create and add the new constraint
        new_con = JuMPeR.build_certain_constraint(m, con, ustar)
        w.debug_printcut && JuMPeR.debug_printcut(rm, m, w, lhs_of_cut, con, new_con)
        push!(new_cons, new_con)
    end
    return new_cons
end

generateReform(w::UDYOracle, m::Model, rm::Model, inds::Vector{Int}) = 0

ext: jl


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

JLLWrappers.@generate_main_file_header("LibCURL")
JLLWrappers.@generate_main_file("LibCURL", UUID("deac9b47-8bc7-5906-a0fe-35ac56dc84c0"))
end  # module LibCURL_jll


text:
module Metrics

using StaticArrays

import ..GPs: sethyperparameters!, gethyperparameters, hyperparametercount

"""
An abstract metric type for those kernels that support the notion of a metric,
e.g., squared exponential kernels. Every child type should implement the
following functions.

    `sqrdistance(metric, x, y)`
    `sethyperparameters!(metric, θ)`
    `gethyperparameters(metric)`
    `hyperparametercount(metric)`
"""
abstract type Metric{T <: Number} end

export Metric

"""
`sqrdistance(metric, x1, x2)` returns the squared distance between `(x1, x2)`
according to the specified metric.
"""
function sqrdistance end

export sqrdistance

"""
`sqrdistance_dθ(metric, x1, x2)` returns the derivative w.r.t. the
hyperparameters of the squared distance between `(x1, x2)` according to the
specified metric.
"""
function sqrdistance_dθ end

export sqrdistance_dθ

"""
`sqrdistance_dx1(metric, x1, x2)` returns the derivative w.r.t. x1 of the
squared distance between `(x1, x2)` according to the specified metric.
"""
function sqrdistance_dx1 end

export sqrdistance_dx1

"""
`sqrdistance_dx2(metric, x1, x2)` returns the derivative w.r.t. x2 of the
squared distance between `(x1, x2)` according to the specified metric.
"""
function sqrdistance_dx2 end

export sqrdistance_dx2

#--- Standard Euclidean metric

"""
Euclidean{N, T}

The standard Euclidean metric in `N` dimensions. (Note that `N=1` corresponds to
an isotropic metric when the input data is multi-dimensional.)
"""
struct Euclidean{N, T} <: Metric{T}
    ℓ::MVector{N, T}
end
Euclidean(ℓ::Union{SVector, Number}) = Euclidean(MVector(ℓ))

export Euclidean

sqrdistance(metric::Euclidean, x1, x2) = sum(((x1 .- x2)./metric.ℓ).^2)

function sqrdistance_dθ(metric::Euclidean{1, T} where {T}, x1, x2)
    return SVector(-2*sum((x1 .- x2).^2)/metric.ℓ[1]^3)
end

function sqrdistance_dθ(metric::Euclidean{N, T} where {N, T}, x1, x2)
    return -2*(x1 .- x2).^2./metric.ℓ.^3
end

function sqrdistance_dx1(metric::Euclidean, x1, x2)
    return 2*(x1 .- x2)./metric.ℓ.^2
end

function sqrdistance_dx2(metric::Euclidean, x1, x2)
    return -2*(x1 .- x2)./metric.ℓ.^2
end

sethyperparameters!(metric::Euclidean, θ) = (metric.ℓ .= θ)
gethyperparameters(metric::Euclidean) = metric.ℓ
hyperparametercount(metric::Euclidean) = length(metric.ℓ)

end


text:
struct NexSliderCursor <: AbstractNexObject
    _nid::NexID

    widthable::IWidthable
    heightable::IHeightable
end

function setproperty!(obj::NexSliderCursor, property::Symbol, new_val)
    _setcommonproperty!(obj, property, new_val)
end

"""
    NexSlider(nexSerial, name; pid=pid, cid=cid)

A Nextion Slider UI control.

On Nextion side, it have the following properties:
- bco
- pco
- wid
- hig
- val
- maxval
- minval
"""
struct NexSlider <: AbstractNexObject
    _nid::NexID

    viewable::IViewable
    numericalvalued::INumericalValued
    colourable::IColourable
    touchable::ITouchable

    function NexSlider(nexSerial::T, name::Name; pid=PageID(), cid=ComponentID()) where {T <: AbstractNexSerial}
        nid = NexID(nexSerial, name, pid, cid)
        #new(nid, IViewable(nid), INumericalValued(nid, Int16, RangeNumber{Int16, 0:typemax(UInt16)}), IColourable(nid), ITouchable(nid))
        new(nid, IViewable(nid), INumericalValued(nid, UInt16, RangeNumber{UInt16, 0:typemax(UInt16)}), IColourable(nid), ITouchable(nid))
    end
end

function setproperty!(obj::NexSlider, property::Symbol, new_val)
    _setcommonproperty!(obj, property, new_val)
end

function getproperty(obj::NexSlider, property::Symbol)
    if property == :cursor
        nid = NexID(obj)
        NexSliderCursor(nid, IWidthable(nid), IHeightable(nid))
    else
        _getcommonproperty(obj, property)
    end
end


text:
using OrdinaryDiffEq, DiffEqDevTools, DiffEqBase, Test
using DiffEqProblemLibrary.ODEProblemLibrary: importodeproblems; importodeproblems()
import DiffEqProblemLibrary.ODEProblemLibrary: prob_ode_linear, prob_ode_2Dlinear

probArr = Vector{ODEProblem}(undef, 2)

probArr[1] = prob_ode_linear
probArr[2] = prob_ode_2Dlinear

function fixed_step_ϕstar(k)
  ∇ = Vector{typeof(k[end][1])}(undef, 3)
  ∇[1] = k[end][1]
  ∇[2] = ∇[1] - k[end-1][1]
  ∇[3] = ∇[2] - k[end-1][1] + k[end-2][1]
  return ∇
end

for i = 1:2
  prob = probArr[i]
  dt = 1//256
  integrator = init(prob,VCAB3(),dt=dt,adaptive=false)
  for i = 1:3
    step!(integrator)
  end
  @test integrator.cache.g == [1,1/2,5/12]*dt
  step!(integrator)
  @test integrator.cache.g == [1,1/2,5/12]*dt
  step!(integrator)
  @test integrator.cache.g == [1,1/2,5/12]*dt
end

for i = 1:2
  prob = probArr[i]
  # VCAB3
  integrator = init(prob,VCAB3(),dt=1//256,adaptive=false)
  for i = 1:3
    step!(integrator)
  end
  # in perform_step, after swaping array using pointer, ϕstar_nm1 points to ϕstar_n
  @test integrator.cache.ϕstar_nm1 == fixed_step_ϕstar(integrator.sol.k)
  step!(integrator)
  @test integrator.cache.ϕstar_nm1 == fixed_step_ϕstar(integrator.sol.k)
  step!(integrator)
  @test integrator.cache.ϕstar_nm1 == fixed_step_ϕstar(integrator.sol.k)

  # VCAB4
  sol1 = solve(prob,VCAB4(),dt=1//256,adaptive=false)
  sol2 = solve(prob,AB4(),dt=1//256)
  @test sol1.u ≈ sol2.u

  # VCAB5
  sol1 = solve(prob,VCAB5(),dt=1//256,adaptive=false)
  sol2 = solve(prob,AB5(),dt=1//256)
  @test sol1.u ≈ sol2.u
end


text:
import .Soss
import .Soss.TransformVariables as TV
import .Soss.SimpleGraphs as SG

export SossMuseProblem

struct SossMuseProblem{A<:AD.AbstractBackend, M<:Soss.AbstractModel, MP<:Soss.AbstractModel} <: AbstractMuseProblem
    autodiff :: A
    model :: M
    model_for_prior :: MP
    xform_z
    xform_θ
    x
    observed_vars
    latent_vars
    hyper_vars
end

@doc doc"""
    SossMuseProblem(model; params, autodiff = ForwardDiffBackend())

Specify a MUSE problem with a
[Soss](https://github.com/cscherrer/Soss.jl) model.

The Soss model should be conditioned on the variables which comprise
the "data", and all other variables should be unconditioned. By
default, any parameter which doesn't depend on another parameter will
be estimated by MUSE, but this can be overridden by passing `params`
as a list of symbols. All other non-conditioned and non-`params`
variables will be considered the latent space.

The `autodiff` parameter should be either
`MuseInference.ForwardDiffBackend()` or
`MuseInference.ZygoteBackend()`, specifying which library to use for
automatic differenation.

## Example

```julia
# toy hierarchical problem, using MeasureTheory distributions
funnel = Soss.@model (σ) begin
    θ ~ MeasureTheory.Normal(0, σ)
    z ~ MeasureTheory.Normal(0, exp(θ/2)) ^ 512
    x ~ For(z) do zᵢ
        MeasureTheory.Normal(zᵢ, 1)
    end
end

sim = rand(funnel(3))
model = funnel(3) | (;sim.x)
prob = SossMuseProblem(model)

# get solution
result = muse(prob, (θ=0,))
```

!!! note

    You can also call [`muse`](@ref), etc... directly on the model, e.g.
    `muse(model, (θ=0,))`, in which case the parameter names `params`
    will be read from the keys of provided the starting point.
"""
function SossMuseProblem(
    model::Soss.ConditionalModel;
    params = leaf_params(model),
    autodiff = ForwardDiffBackend()
)
    x = model.obs
    !isempty(x) || error("Model must be conditioned on observed data.")
    sim = rand(model)
    observed_vars = keys(x)
    hyper_vars = (((k, eltype(sim[k])) for k in keys(sim) if k in params)...,)
    latent_vars = keys(delete(sim, (observed_vars..., first.(hyper_vars)...)))
    model_for_prior = Soss.likelihood(Soss.Model(model), first.(hyper_vars)...)(Soss.argvals(model))
    xform_z = Soss.xform(model | select(sim, first.(hyper_vars)))
    xform_θ = Soss.xform(model | select(sim, latent_vars))
    SossMuseProblem(
        autodiff,
        model,
        model_for_prior,
        xform_z,
        xform_θ,
        x,
        observed_vars,
        latent_vars,
        hyper_vars
    )
end

function leaf_params(model)
    D = Soss.digraph(Soss.Model(model))
    [x for x in SG.vlist(D) if isempty(filter(!in(keys(Soss.argvals(model))), SG.in_neighbors(D, x)))]
end

function transform_θ(prob::SossMuseProblem, θ)
    TV.inverse(prob.xform_θ, _namedtuple(θ))
end

function inv_transform_θ(prob::SossMuseProblem, θ)
    ComponentVector(TV.transform(prob.xform_θ, θ))
end

function logPriorθ(prob::SossMuseProblem, θ::ComponentVector, ::UnTransformedθ)
    Soss.logdensity(prob.model_for_prior(_namedtuple(θ)))
end
function logPriorθ(prob::SossMuseProblem, θ::AbstractVector, ::Transformedθ)
    logPriorθ(prob, inv_transform_θ(prob, θ), UnTransformedθ())
end

function ∇θ_logLike(prob::SossMuseProblem, x, z::AbstractVector, θ::ComponentVector, ::UnTransformedθ)
    like = prob.model | (;x..., TV.transform(prob.xform_z, z)...)
    first(AD.gradient(prob.autodiff, θ -> Soss.logdensity(like, _namedtuple(θ)), θ))
end
function ∇θ_logLike(prob::SossMuseProblem, x, z::AbstractVector, θ::AbstractVector, ::Transformedθ)
    like = prob.model | (;x..., TV.transform(prob.xform_z, z)...)
    first(AD.gradient(prob.autodiff, θ -> Soss.logdensity(like, _namedtuple(inv_transform_θ(prob, θ))), θ))
end

function logLike_and_∇z_logLike(prob::SossMuseProblem, x, z, θ)
    first.(AD.value_and_gradient(prob.autodiff, z -> Soss.logdensity(prob.model | (;x..., _namedtuple(θ)...), TV.transform(prob.xform_z, z)), z))
end

function sample_x_z(prob::SossMuseProblem, rng::AbstractRNG, θ)
    sim = Soss.predict(rng, prob.model, _namedtuple(θ))
    x = select(sim, prob.observed_vars)
    z = TV.inverse(prob.xform_z, select(sim, prob.latent_vars))
    (;x, z)
end

# all user-provided θ pass through this first so get it right type /
# order so the user can be lazy
function standardizeθ(prob::SossMuseProblem, θ::NamedTuple)
    Set(keys(θ)) == Set(first.(prob.hyper_vars)) || error("Expected θ to have keys: $(first.(prob.hyper_vars)).")
    ComponentVector((;(k => T.(θ[k]) for (k,T) in prob.hyper_vars)...))
end
function standardizeθ(prob::SossMuseProblem, θ::Number)
    standardizeθ(prob, (;(prob.hyper_vars[1][1] => θ)))
end

_params_from_θ₀(θ₀::Number) = (:θ,)
_params_from_θ₀(θ₀) = keys(θ₀)

function muse!(result::MuseResult, model::Soss.ConditionalModel, θ₀ = result.θ; kwargs...)
    muse!(result, SossMuseProblem(model, params=_params_from_θ₀(θ₀)), θ₀; kwargs...)
end
function get_J!(result::MuseResult, model::Soss.ConditionalModel, θ₀ = result.θ; kwargs...)
    get_J!(result, SossMuseProblem(model, params=_params_from_θ₀(θ₀)), θ₀; kwargs...)
end
function get_H!(result::MuseResult, model::Soss.ConditionalModel, θ₀ = result.θ; kwargs...)
    get_H!(result, SossMuseProblem(model, params=_params_from_θ₀(θ₀)), θ₀; kwargs...)
end


text:
function MP.differentiate(m::Monomial{C}, x::PolyVar{C}) where C
    i = findfirst(isequal(x), _vars(m))
    if (i == nothing || i == 0) || m.z[i] == 0
        zeroterm(m)
    else
        z = copy(m.z)
        z[i] -= 1
        m.z[i] * Monomial(_vars(m), z)
    end
end

function MP.differentiate(p::Polynomial{C, T}, x::PolyVar{C}) where {C, T}
    # grlex order preserved
    i = something(findfirst(isequal(x), _vars(p)), 0)
    S = typeof(zero(T) * 0)
    if iszero(i)
        zero(Polynomial{C, S})
    else
        keep = findall(z -> z[i] > 0, p.x.Z)
        Z = copy.(p.x.Z[keep])
        a = Vector{S}(undef, length(keep))
        for j in 1:length(Z)
            a[j] = p.a[keep[j]] * Z[j][i]
            Z[j][i] -= 1
        end
        Polynomial(a, MonomialVector(_vars(p), Z))
    end
end


text:
# Assumes that the ptr does not alias
# Represents this via not actually loading/storing value to the ptr until some finalization step.
struct Reference{T,V}
    ptr::Ptr{T}
    value::V
end

Base.:(+)(r::Reference, x) = Reference(r.ptr, vadd(r.value, x))
Base.:(-)(r::Reference, x) = Reference(r.ptr, vsub(r.value, x))
Base.:(*)(r::Reference, x) = Reference(r.ptr, vmul(r.value, x))
Base.:(/)(r::Reference, x) = Reference(r.ptr, vdiv(r.value, x))
Base.:(+)(x, r::Reference) = Reference(r.ptr, vadd(x, r.value))
Base.:(-)(x, r::Reference) = Reference(r.ptr, vsub(x, r.value))
Base.:(*)(x, r::Reference) = Reference(r.ptr, vmul(x, r.value))
Base.:(/)(x, r::Reference) = Reference(r.ptr, vdiv(x, r.value))
Base.exp(r::Reference) = Reference(r.ptr, exp(r.value))
Base.log(r::Reference) = Reference(r.ptr, log(r.value))


text:
"""
   extract(x, points; order, atol)

Extracts the value of `Raster` or `RasterStack` at given points, returning
a vector of `NamedTuple` with columns for the point dimensions and layer
value/s.

Note that if objects have more dimensions than the length of the point tuples,
sliced arrays or stacks will be returned instead of single values.

# Arguments

- `x`: a `Raster` or `RasterStack` to extract values from.
- `points`: multiple `Vector`s of point values, a `Vector{Tuple}`,
    or a single `Tuple` or `Vector`. `points` can also be a Tables.jl compatible
    table, in which case `order` may need to specify the keys.

# Keywords

- `order`: a tuple of `Dimension` connecting the order of the points to the array
    axes, such as `(X, Y)`, with a defaut `(XDim, YDim, ZDim)` order.
    If `points` is a table, `order` should be a `Tuple` of `Dimension`/`Symbol` pairs
    like `(X => :xcol, Y => :ycol)`. This will be automatically detected wherever
    possible, assuming the keys match the dimensions of the object `x`.
- `atol`: a tolorerance for floating point lookup values for when the `LookupArray`
    contains `Points`. `atol` is ignored for `Intervals`.

Note: extracting polygons in a `GeoInterface.AbstractGeometry` is not yet supported,
but will be in future.

# Example

Here we extact points matching the occurrence of the Mountain Pygmy Possum,
_Burramis parvus_. This could be used to fit a species distribution lookupl.

```jldoctest
using Rasters, GBIF, CSV

# Get a stack of BioClim layers, and replace missing values with `missing`
st = RasterStack(WorldClim{BioClim}, (1, 3, 5, 7, 12))[Band(1)] |> replace_missing

# Download some occurrence data
obs = GBIF.occurrences("scientificName" => "Burramys parvus", "limit" => 5)

# use `extract` to get values for all layers at each observation point.
points = map(o -> (o.longitude, o.latitude), obs)
vals = extract(st, points)

# output
5-element Vector{NamedTuple{(:X, :Y, :bio1, :bio3, :bio5, :bio7, :bio12)}}:
 (X = missing, Y = missing, bio1 = missing, bio3 = missing, bio5 = missing, bio7 = missing, bio12 = missing)
 (X = 147.096394, Y = -36.935687, bio1 = 9.408354f0, bio3 = 40.790546f0, bio5 = 22.39425f0, bio7 = 23.0895f0, bio12 = 1292.0f0)
 (X = 148.450743, Y = -35.999643, bio1 = 8.269542f0, bio3 = 41.030262f0, bio5 = 21.4485f0, bio7 = 23.858f0, bio12 = 1440.0f0)
 (X = 148.461854, Y = -36.009001, bio1 = 6.928167f0, bio3 = 41.78015f0, bio5 = 20.18025f0, bio7 = 23.69975f0, bio12 = 1647.0f0)
 (X = 148.459452, Y = -36.002648, bio1 = 6.928167f0, bio3 = 41.78015f0, bio5 = 20.18025f0, bio7 = 23.69975f0, bio12 = 1647.0f0)

```
"""
function extract(A::RasterStackOrArray, points::NTuple{<:Any,<:AbstractVector}; kw...)
    extract(A, zip(points...); kw...)
end
function extract(A::RasterStackOrArray, points::GI.AbstractGeometry; kw...)
    extract(A, _flat_nodes(GI.coordinates(points)); kw...)
end
function extract(A::RasterStackOrArray, points::AbstractVector{<:Tuple}; kw...)
    extract.(Ref(A), points; kw...)
end
function extract(A::RasterStackOrArray, points::AbstractVector{Union{<:AbstractVector{<:Union{Real,Missing}},Missing}}; kw...)
    extract.(Ref(A), points; kw...)
end
function extract(A::RasterStackOrArray, data; order=nothing, kw...)
    if Tables.istable(data)
        order = isnothing(order) ? _auto_dim_columns(dims(A), data) : order
        rows = Tables.rows(data)
        point_dims = map(p -> DD.basetypeof(p[1])(p[2]), order)
        point_keys = map(val, point_dims)
        map(rows) do row
            point_vals = map(pk -> row[pk], point_keys)
            extract(A, point_vals; order=map(first, order), point_keys, kw...)
        end
    else
        order = isnothing(order) ? DEFAULT_POINT_ORDER : order
        map(data) do point
            extract(A, point; order, kw...)
        end
    end
end
extract(A::RasterStackOrArray, points::Missing; kw...) = missing
function extract(
    x::RasterStackOrArray, point::Union{Tuple,AbstractVector{<:Union{Missing,<:Real}}};
    order=(XDim, YDim, ZDim),
    point_keys=map(DD.dim2key, dims(x, order)),
    layer_keys=_layer_keys(x),
    atol=nothing
)
    # Get the actual dimensions available in the object
    # Usually this will be `X` and `Y`, but `Z` as well if it exists.
    ordered_dims = dims(x, order)
    length(point) == length(ordered_dims) || throw(ArgumentError("Length of `point` does not match dims. Pass `order` dims manually"))
    point = ntuple(i -> point[i], length(ordered_dims))
    dimtypes = map(DD.basetypeof, ordered_dims)

    # Extract the values
    if any(map(ismissing, point))
        point_vals = map(_ -> missing, ordered_dims)
        layer_vals = map(_ -> missing, layer_keys)
    else
        selectors = map((d, x) -> _at_or_contains(d, x, atol), ordered_dims, point)
        point_vals = map(val ∘ val, selectors)
        layer_vals = if DD.hasselection(x, selectors)
            x isa Raster ? (x[selectors...],) : x[selectors...]
        else
            map(_ -> missing, layer_keys)
        end
    end
    return NamedTuple{(point_keys..., layer_keys...)}((point_vals..., layer_vals...))
end

_layer_keys(A::AbstractRaster) = cleankeys(name(A))
_layer_keys(A::AbstractRasterStack) = cleankeys(keys(A))


text:
module TestTools
using DataFrames, CSV

using Test

include("../src/Tools.jl")
import .Tools
export runtests

function test_make_Y_bus()

    VN = 12500
    Sb = 1e6
    Yb = Sb / (VN^2)

    target_pu = [
        20-40im -20+40im 0+0im
        -20+40im 48.5714-97.1429im -28.5714+57.1429im
        0+0im -28.5714+57.1429im 28.5714-56.8429im
    ]
    target = target_pu * Yb

    data = DataFrame(CSV.File("case3_dist.csv"))

    result = Tools.make_Y_bus(data, VN)
    @testset "make_Y_bus" begin
        @testset "Y[$i, $j]" for i = 1:3, j = 1:3
            @test target[i, j] ≈ result[i, j] atol = 1e-4
        end
    end
end

function runtests()
    @testset "Tools" begin
        test_make_Y_bus()
    end
end

end

TestTools.runtests()


text:
RO_data = Load_RO_data()
weatherts_F,GPP_data_F = create_weather_struct_RO(RO_data;stand_type="Fertilized")
weatherts_C,GPP_data_C = create_weather_struct_RO(RO_data;stand_type="Control")

data_F,data_C = Create_RoData_C_F(GPP_data_F,GPP_data_C,weatherts_F,weatherts_C)

treepar = TreePar()
cons = CCPH.Constants()

A_c_F = Float64[]
A_c_C = Float64[]
A_c_F_calc = Float64[]
GPP_F_calc = Float64[]
E_c_F = Calctranpiration.(weatherts_F.VPD.*weatherts_F.daylight/(7*24*3600),CCPH.CalcSₑ.(weatherts_F.θₛ))
E_c_C = Calctranpiration.(weatherts_C.VPD.*weatherts_C.daylight/(7*24*3600),CCPH.CalcSₑ.(weatherts_C.θₛ))

g_c = 4 #mm s⁻¹
J_max_opt = 184.0e-6
J_max_opt = 400.0e-6

E_c_Calc = g_c*weatherts_C.VPD/10^5.0.*weatherts_F.daylight/7*(cons.ρ_vapor/cons.ρ_H2O)
for i = 1:length(weatherts_F.date)

    photorates = CCPH.PhotoKineticRates()
    env = CCPH.EnvironmentStruct(weatherts_F,i)
    photopar = PhotoPar(photorates,env.Tₐ)
    Jₘₐₓ = J_max_opt*photopar.b_Jmax
    Iᵢ = CCPH.Calc_Iᵢ(env.I₀,treepar)
    gₜ = CCPH.Calc_gₜ(data_F.gₛ[i],treepar)

    A_c = CCPH.Farquhar(gₜ,Iᵢ,Jₘₐₓ,photopar,env)[1]

    push!(A_c_F_calc,A_c)

    ind = Find_data_ind(i)

    LAI_F = data_F.Wf[ind]*data_F.N[ind]/treepar.LMA
    N_F = data_F.N[ind]
    push!(A_c_F,Est_C_assimilation(GPP_data_F[i],LAI_F,N_F,weatherts_F,i))

    LAI_F = data_C.Wf[ind]*data_C.N[ind]/treepar.LMA
    N_F = data_C.N[ind]
    push!(A_c_C,Est_C_assimilation(GPP_data_C[i],LAI_F,N_F,weatherts_C,i)[1])

    Xₜ = weatherts_F.acclimation_fac[i]
    daylight = weatherts_F.daylight[i]/7
    GPP = A_c*(1-exp(-treepar.k*LAI_F))/treepar.k*Xₜ*cons.M_C*daylight*10^3 #gC m⁻² ground area day⁻¹
    push!(GPP_F_calc,GPP)
end

plot(weatherts_F.date,E_c_F,xlabel="Date",ylabel="Canopy tranpiraiton (mm day⁻¹)",label="F")
plot!(weatherts_C.date,E_c_C,label="C")
plot!(weatherts_F.date,E_c_Calc,label="g_c = 4 mm s⁻¹")
savefig("plots/Ec.svg")

plot(weatherts_F.date,data_F.gₛ,xlabel="Date",
ylabel="Stomatal conductance (mol C m⁻² leaf area s⁻¹)",label="F")
plot!(weatherts_C.date,data_C.gₛ,label="C")
savefig("plots/gs.svg")

plot(weatherts_F.date,data_F.GPP,xlabel="Date",
ylabel="GPP (g C m⁻² ground area day⁻¹)",label="F data")
plot!(weatherts_C.date,data_C.GPP,label="C data")
plot!(weatherts_F.date,GPP_F_calc,label="F calc")
savefig("plots/GPP.svg")

pl1 = plot(weatherts_F.date,data_F.GPP,xlabel="Date",ylabel="GPP",label="F data")
plot!(weatherts_F.date,GPP_F_calc)
pl2 = plot(weatherts_F.date,E_c_F,ylabel="E_c",label="F")
pl3 = plot(weatherts_F.date,data_F.gₛ,ylabel="gₛ",label="F")
pl4 = plot(weatherts_F.date,weatherts_F.VPD,ylabel="VPD",label="F")
pl5 = plot(weatherts_F.date,weatherts_F.θₛ,ylabel="θₛ",label="F")
pl6 = plot(weatherts_F.date,A_c_F_calc,ylabel="A_c",label="F")
pl7 = plot(weatherts_F.date,GPP_F_calc./A_c_F_calc,ylabel="Γ")

plot(pl1,pl6,pl7,pl3,layout=(4,1),legends=false)
savefig("plots/fig1.svg")
plot(pl2,pl3,pl4,pl5,layout=(4,1),legends=false)
savefig("plots/fig2.svg")


text:
# Autogenerated wrapper script for MongoC_jll for i686-linux-gnu
export libbson, libmongoc

using OpenSSL_jll
using Zlib_jll
using Zstd_jll
using snappy_jll
JLLWrappers.@generate_wrapper_header("MongoC")
JLLWrappers.@declare_library_product(libbson, "libbson-1.0.so.0")
JLLWrappers.@declare_library_product(libmongoc, "libmongoc-1.0.so.0")
function __init__()
    JLLWrappers.@generate_init_header(OpenSSL_jll, Zlib_jll, Zstd_jll, snappy_jll)
    JLLWrappers.@init_library_product(
        libbson,
        "lib/libbson-1.0.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libmongoc,
        "lib/libmongoc-1.0.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
### Day 7 - Problem 2
## Internet Protocol Version 7
## Authors: Thanasis Georgiou, Marina Grigoropoulou

# Looks for an ABA identifier. Compatible with BAB
function findABAs(str)
    foundABAs = []

    len = length(str)
    if len >= 3
        for i = 1:len - 2
            if str[i] == str[i + 2]
                push!(foundABAs, str[i:i + 2])
            end
        end
    end

    return foundABAs
end

# Are the ABA and BAB identifiers a pair
function isPair(aba, bab)
    return collect(aba) == [bab[2], bab[1], bab[2]]
end

# Read input
f = open("input.txt")
input = readlines(f)
close(f)

# How many support TLS
count = 0

# For each IP
for ip in input
    # Split in string that are contained in brackets
    # and ones that are not
    bracketArr = matchall(r"\[\w*\]", ip)
    map!(str -> str[2:end-1], bracketArr)

    inBracket = join(bracketArr, '-')
    outBracket = replace(ip, r"\[\w*\]", "-")

    # Find ABAs and BABs
    inABAs = findABAs(inBracket)
    outABAs = findABAs(outBracket)

    # Check if any pair is valid
    found = false
    for i = 1:length(inABAs)
        for j = 1:length(outABAs)
            if isPair(inABAs[i], outABAs[j])
                found = true
                break
            end
        end

        if found
            break
        end
    end

    # This IP was valid
    if found
        count += 1
    end
end

println("Valid IPs: $count")


text:
# External packages
using Revise, DelimitedFiles, SymPy

# Add Exasim to Julia search path
cdir = pwd(); ii = findlast("Exasim", cdir);
include(cdir[1:ii[end]] * "/Installation/setpath.jl");

# Exasim packages
using Preprocessing, Mesh, Gencode, Postprocessing

# create pde structure and mesh structure
pde, mesh = Preprocessing.initializeexasim();

# Define PDE model: governing equations, initial solutions, and boundary conditions
pde.model = "ModelD";            # ModelC, ModelD, ModelW
include("pdemodel.jl");          # include the PDE model file

# Set discretization parameters, physical parameters, and solver parameters
pde.porder = 3;               # polynomial degree
pde.torder = 2;               # time-stepping order of accuracy
pde.nstage = 2;               # time-stepping number of stages
pde.dt = 0.002*ones(10000);   # time step sizes
pde.soltime = collect(1:length(pde.dt)); # steps at which solution are collected
pde.visdt = pde.dt[1]; # visualization timestep size

gam = 1.4;                      # specific heat ratio
Re = 10000;                     # Reynolds number
Pr = 0.71;                      # Prandtl number
Minf = 0.4;                     # Mach number
alpha = 5*pi/180;               # angle of attack
rinf = 1.0;                     # freestream density
ruinf = cos(alpha);             # freestream horizontal velocity
rvinf = sin(alpha);             # freestream vertical velocity
pinf = 1/(gam*Minf^2);          # freestream pressure
rEinf = 0.5+pinf/(gam-1);       # freestream energy
pde.physicsparam = [gam Re Pr Minf rinf ruinf rvinf rEinf];
pde.tau = [2.0];          # DG stabilization parameter
pde.GMRESrestart=30;            # number of GMRES restarts
pde.linearsolvertol=0.0001;     # GMRES tolerance
pde.linearsolveriter=31;        # number of GMRES iterations
pde.precMatrixType=2;           # preconditioning type
pde.NLtol = 1e-7;               # Newton tolerance
pde.NLiter=3;                   # Newton iterations

# Choose computing platform and set number of processors
#pde.platform = "gpu";           # choose this option if NVIDIA GPUs are available
pde.mpiprocs = 4;                # number of MPI processors

# read a grid from a file
mesh.p,mesh.t = Mesh.readmesh("grid.bin",0);
# expressions for domain boundaries
mesh.boundaryexpr = [p -> (sqrt.((p[1,:].-0.5).*(p[1,:].-0.5)+(p[2,:].*p[2,:])) .< 3), p -> (p[1,:] .< 20)];
mesh.boundarycondition = [1 2];
# expressions for curved boundaries
mesh.curvedboundary = [1 0];
mesh.curvedboundaryexpr = [p -> p[2,:].^2-(5*0.01*12*(0.29690*sqrt.(abs.(p[1,:]))-0.12600*p[1,:]-0.35160*p[1,:].^2+0.28430*p[1,:].^3-0.10150*p[1,:].^4)).^2, p -> 0];

# call exasim to generate and run C++ code to solve the PDE model
sol, pde, mesh,~,~,~,~  = Postprocessing.exasim(pde,mesh);

# visualize the numerical solution of the PDE model using Paraview
pde.visscalars = ["density", 1, "energy", 4];  # list of scalar fields for visualization
pde.visvectors = ["momentum", [2, 3]]; # list of vector fields for visualization
Postprocessing.vis(sol,pde,mesh); # visualize the numerical solution
print("Done!");


text:
using Reproject: parse_input_data, parse_output_projection
@testset "input parser" begin
    fname = tempname() * ".fits"
    f = FITS(fname, "w")
    inhdr = FITSHeader(["FLTKEY", "INTKEY", "BOOLKEY", "STRKEY", "COMMENT",
                        "HISTORY"],
                       [1.0, 1, true, "string value", nothing, nothing],
                       ["floating point keyword",
                        "",
                        "boolean keyword",
                        "string value",
                        "this is a comment",
                        "this is a history"])

    indata = reshape(Float32[1:100;], 5, 20)
    write(f, indata; header=inhdr)

    @testset "ImageHDU type" begin
        result = parse_input_data(f[1])
        @test result[1] isa Array
        @test result[2] isa WCSTransform
    end

    @testset "data matrix and WCSTransform tuple" begin
        wcs = WCSTransform(2;
                          cdelt = [-0.066667, 0.066667],
                          ctype = ["RA---AIR", "DEC--AIR"],
                          crpix = [-234.75, 8.3393],
                          crval = [0., -90],
                          pv    = [(2, 1, 45.0)])
        result = parse_input_data((indata, wcs))
        @test result[1] isa Array
        @test result[2] isa WCSTransform
    end

    @testset "Single HDU FITS file" begin
        result = parse_input_data(f, 1)
        @test result[1] isa Array
        @test result[2] isa WCSTransform
    end
    close(f)

    @testset "String filename input" begin
        result = parse_input_data(fname, 1)
        @test result[1] isa Array
        @test result[2] isa WCSTransform
    end

    f = FITS(fname, "w")
    write(f, indata; header=inhdr)
    write(f, indata; header=inhdr)

    @testset "Multiple HDU FITS file" begin
        result = parse_input_data(f, 2)
        @test result[1] isa Array
        @test result[2] isa WCSTransform

        close(f)
        result = parse_input_data(fname, 1)
        @test result[1] isa Array
        @test result[2] isa WCSTransform
    end
end

@testset "output parser" begin
    fname = tempname() * ".fits"
    f = FITS(fname, "w")
    inhdr = FITSHeader(["FLTKEY", "INTKEY", "BOOLKEY", "STRKEY", "COMMENT",
                        "HISTORY"],
                       [1.0, 1, true, "string value", nothing, nothing],
                       ["floating point keyword",
                        "",
                        "boolean keyword",
                        "string value",
                        "this is a comment",
                        "this is a history"])

    indata = reshape(Float32[1:100;], 5, 20)
    write(f, indata; header=inhdr)

    @testset "ImageHDU type" begin
        result = parse_output_projection(f[1], (12,12))
        @test result[1] isa WCSTransform
        @test result[2] isa Tuple
        @test_throws DomainError parse_output_projection(f[1], ())
    end
    close(f)

    @testset "String filename" begin
        result = parse_output_projection(fname, 1)
        @test result[1] isa WCSTransform
        @test result[2] isa Tuple
    end

    wcs = WCSTransform(2;
                          cdelt = [-0.066667, 0.066667],
                          ctype = ["RA---AIR", "DEC--AIR"],
                          crpix = [-234.75, 8.3393],
                          crval = [0., -90],
                          pv    = [(2, 1, 45.0)])

    @testset "WCSTransform input" begin
        result = parse_output_projection(wcs, (12,12))
        @test result[1] isa WCSTransform
        @test result[2] isa Tuple
        @test_throws DomainError parse_output_projection(wcs, ())
    end
end


text:
using ConstrainedDynamics

# Parameters
joint_axis = [1.0;0.0;0.0]

length1 = 0.5
width, depth = 0.5, 0.5
box1 = Box(width, depth, length1, 1., color = RGBA(1., 1., 0.))
b1 = Box(0.1, 0.1, .1, .1, color = RGBA(1., 0., 0.))

corners = [
    [[length1 / 2;length1 / 2;-length1 / 2]]
    [[length1 / 2;-length1 / 2;-length1 / 2]]
    [[-length1 / 2;length1 / 2;-length1 / 2]]
    [[-length1 / 2;-length1 / 2;-length1 / 2]]
    [[length1 / 2;length1 / 2;length1 / 2]]
    [[length1 / 2;-length1 / 2;length1 / 2]]
    [[-length1 / 2;length1 / 2;length1 / 2]]
    [[-length1 / 2;-length1 / 2;length1 / 2]]
]

# Links
origin = Origin{Float64}()

link1 = Body(box1)
links = [link1;[Body(b1) for i = 1:8]]

# Constraints
ineqcs = [InequalityConstraint(Impact(links[i + 1], [0;0;1.0])) for i = 1:8]

joint0to1 = EqualityConstraint(OriginConnection(origin, link1))
eqcs = [joint0to1;[EqualityConstraint(Fixed(link1, links[i + 1], corners[i], zeros(3))) for i = 1:8]]

shapes = [box1;b1]

mech = Mechanism(origin, links, eqcs, ineqcs, shapes = shapes)
setPosition!(mech,link1,x = [0.;-2;1.5])
for i = 1:8
    setPosition!(mech, link1, links[i + 1], p1 = corners[i])
end

setVelocity!(mech,link1,v = [0;3;7.],ω = (rand(3) .- 0.5) * 100)
for i = 1:8
    setVelocity!(mech, link1, links[i + 1], p1 = corners[i])
end

simulate!(mech,save = true)
visualize!(mech)


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

using Markdown
using InteractiveUtils

# ╔═╡ a2850860-f73b-11ea-159a-ad5bc1704c57
using Images, FileIO

# ╔═╡ 7302a090-f739-11ea-3fd7-2786faf0a7bd
md"## Data utils
Contains helper functions to deal with data extraction, download etc.

**TODO**
Functions for downloading, saving and other actions are pending"

# ╔═╡ cad848b0-f739-11ea-0343-177efd090dc6
md"Many functions from `Datadeps` package can be extended for data handeling.
One such function is `Unpack`.

##### unpack
`Unpack` takes the archive name and uncompresses it. Optionally a destination path can be provided or by default the uncompression is done in the present working directory. The archive file gets deleted by default once the files are extracted. However, this can be changed."

# ╔═╡ 98ef39c0-f73a-11ea-0426-dda23207aeb9
begin
import DataDeps: unpack_cmd, unpack

function unpack(f; filepath=pwd(), keep_originals=false)
        run(unpack_cmd(f, filepath, last(splitext(f)), last(splitext(first(splitext(f))))))
        rm("pax_global_header"; force=true)
        !keep_originals && rm(f)
end

end

# ╔═╡ c9e2a8a0-f73a-11ea-11a1-55fa9f327373
md"This can then be used to extract a compressed file."

# ╔═╡ d8a383a0-f73a-11ea-3fbd-cf9783b8ad3d
#unpack(<path to the archive file, filepath=<destination path>)

# ╔═╡ 2985b0e0-f73b-11ea-0a83-c52b1117dece
md"##### load_image
Reads an image file and loads it to memory.

It is a thin wrapper over the Images.load functions and is more of a convenience function.
"

# ╔═╡ a01c93e0-f73b-11ea-3e76-4fe742ac8301
function load_image(path::AbstractString)
	if isfile(path)
		return load(path)
	end
end

# ╔═╡ d99a2330-f73b-11ea-0fa1-fbd0524c58ca
md"##### image2array

A thin wrapper over `Images.channelview` function. convenience function to get the fixed poitn array of an image"

# ╔═╡ 129a4520-f73c-11ea-1272-cbe8a3ca2f08
function image2array(image)
	channelview(image)
end

# ╔═╡ Cell order:
# ╟─7302a090-f739-11ea-3fd7-2786faf0a7bd
# ╟─cad848b0-f739-11ea-0343-177efd090dc6
# ╠═98ef39c0-f73a-11ea-0426-dda23207aeb9
# ╟─c9e2a8a0-f73a-11ea-11a1-55fa9f327373
# ╠═d8a383a0-f73a-11ea-3fbd-cf9783b8ad3d
# ╟─2985b0e0-f73b-11ea-0a83-c52b1117dece
# ╠═a2850860-f73b-11ea-159a-ad5bc1704c57
# ╠═a01c93e0-f73b-11ea-3e76-4fe742ac8301
# ╟─d99a2330-f73b-11ea-0fa1-fbd0524c58ca
# ╠═129a4520-f73c-11ea-1272-cbe8a3ca2f08


text:
# MIT License
#
# Copyright (c) 2020 Joshua E Gentry

# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.
using Blake3Hash

@testset "Hash - standard" begin
    @testset "Length: ($(test["input_len"]))" for test in CASES
        local case = Case(test)

        local hasher = Blake3Hash.SAImpl.Blake3Ctx()
        local result = Vector{UInt8}(undef, length(case.hash))

        Blake3Hash.SAImpl.update!(hasher, case.input)
        Blake3Hash.SAImpl.digest(hasher, result)

        @test case.hash == result
    end
end

@testset "Hash - keyed" begin
    @testset "Length: ($(test["input_len"]))" for test in CASES
        local case = Case(test)

        local hasher = Blake3Hash.SAImpl.Blake3Ctx(KEY)
        local result = Vector{UInt8}(undef, length(case.hash))

        Blake3Hash.SAImpl.update!(hasher, case.input)
        Blake3Hash.SAImpl.digest(hasher, result)

        @test case.keyed == result
    end
end

@testset "Hash - derive" begin
    @testset "Length: ($(test["input_len"]))" for test in CASES
        local case = Case(test)

        local hasher = Blake3Hash.SAImpl.Blake3Ctx(CONTEXT)
        local result = Vector{UInt8}(undef, length(case.hash))

        Blake3Hash.SAImpl.update!(hasher, case.input)
        Blake3Hash.SAImpl.digest(hasher, result)

        @test case.derive == result
    end
end

@testset "Hash - 256" begin
    @testset "Length: ($(test["input_len"]))" for test in CASES
        local case = Case(test)

        local hasher = Blake3Hash.SAImpl.Blake3Ctx()

        Blake3Hash.SAImpl.update!(hasher, case.input)

        local result = Blake3Hash.SAImpl.digest(hasher)

        @test case.hash[1:32] == result
    end
end

@testset "Hash - digest" begin
    @testset "Length: ($(test["input_len"]))" for test in CASES
        local case = Case(test)

        local hasher = Blake3Hash.SAImpl.Blake3Ctx()
        Blake3Hash.SAImpl.update!(hasher, case.input)

        @test Blake3Hash.SAImpl.digest(hasher, 1) == case.hash[1:1]
        @test Blake3Hash.SAImpl.digest(hasher, 2) == case.hash[1:2]
        @test Blake3Hash.SAImpl.digest(hasher, 3) == case.hash[1:3]
        @test Blake3Hash.SAImpl.digest(hasher, 4) == case.hash[1:4]

        @test Blake3Hash.SAImpl.digest(hasher, 30) == case.hash[1:30]
        @test Blake3Hash.SAImpl.digest(hasher, 31) == case.hash[1:31]
        @test Blake3Hash.SAImpl.digest(hasher, 32) == case.hash[1:32]
        @test Blake3Hash.SAImpl.digest(hasher, 33) == case.hash[1:33]
        @test Blake3Hash.SAImpl.digest(hasher, 34) == case.hash[1:34]
    end
end

@testset "Hash - update" begin
    @testset "Length: ($(test["input_len"]))" for test in CASES
        local case = Case(test)

        local hasher = Blake3Hash.SAImpl.Blake3Ctx()
        local result = Vector{UInt8}(undef, length(case.hash))

        local data = case.input

        while length(data) > 0
            local take = min(rand(1:4), length(data))
            Blake3Hash.SAImpl.update!(hasher, data[1:take])
            data = data[(take+1):end]
        end

        Blake3Hash.SAImpl.digest(hasher, result)

        @test case.hash == result
    end
end

@testset "Hash - stress" begin
    data   = rand(UInt8, TOTAL_SIZE)
    hasher = Blake3Hash.SAImpl.Blake3Ctx()
    Blake3Hash.SAImpl.update!(hasher, data)
    expected = Blake3Hash.SAImpl.digest(hasher)

    for chunks in StressIter(data)
        hasher = Blake3Hash.SAImpl.Blake3Ctx()

        for chunk in chunks
            Blake3Hash.SAImpl.update!(hasher, chunk)
        end

        @test Blake3Hash.SAImpl.digest(hasher) == expected
    end
end


text:
"""
    GMMModel

Abstract type for GMM models.
"""
abstract type GMMModel end

"""
    get_gi(model::GMMModel)

Returns a function `gi(θ)` where the moment condition for a GMM model is

``E[g_i(\\theta)] = 0``

`gi(θ)` should be a number of observations by number of moment conditions matrix.
"""
function get_gi(model::GMMModel)
  error("get_gi not defined for model type $(typeof(model))")
end

"""
     number_parameters(model::GMMModel)

Number of parameters (dimension of θ) for a GMMModel.
"""
function number_parameters(model::GMMModel)
  error("number_parameters not defined for model type $(typeof(model))")
end

"""
     number_observations(model::GMMModel)

Number of observations (rows of gi(θ)) for a GMMModel
"""
function number_observations(model::GMMModel)
  θ = ones(number_parameters(model))
  gi = get_gi(model)
  return(size(gi(θ),1))
end

"""
    number_moments(model::GMModel)

Number of moments (columns of gi(θ)) for a GMMModel
"""
function number_moments(model::GMMModel)
  θ = ones(number_parameters(model))
  gi = get_gi(model)
  return(size(gi(θ),2))
end

"""
    gmm_constraints(model::GMMModel)

Returns constraint as function or parameters, `θ`, where `c(θ) = 0`.
"""
function gmm_constraints(model::GMMModel)
  return nothing
end

"""
     gmm_objective(gi::Function, W=I)

Returns ``Q(θ) = n (1/n \\sum_i g_i(θ)) W (1/n \\sum_i g_i(θ))' ``
"""
function gmm_objective(gi::Function, W=I)
  function(θ)
    g = gi(θ)
    m = mean(g,dims=1)
    (size(g)[1]*( m*W*m')[1])
  end
end

"""
     cue_objective(gi::Function)

Returns the CUE objective function for moment functions gi.

``Q(θ) = n (1/n \\sum_i g_i(θ)) \\widehat{cov}(g_i(θ))^{-1} (1/n \\sum_i g_i(θ))' ``

Calculates cov(gi(θ)) assuming observations are independent.
"""
function cue_objective(gi::Function)
  function(θ)
    g = gi(θ)
    m = mean(g,dims=1)
    W = inv(cov(g))
    (size(g)[1]*( m*W*m')[1])
  end
end

"""
    gmm_objective(model::GMMModel, W=I)

Returns the GMM objective function with weighting matrix `W` for `model`.
"""
function gmm_objective(model::GMMModel, W=I)
  gmm_objective(get_gi(model), W)
end

"""
    cue_objective(model::GMMModel)

Returns the CUE objective function for `model`.

Calculated weighting matrix assuming observations are independent.
"""
function cue_objective(model::GMMModel)
  cue_objective(get_gi(model))
end

"""
    gel_jump_problem(model::GMMModel, h::Function=log)

Returns JuMP problem for generalized empirical likelihood estimation of `model`.
`h` is a generalized likelihood function.
"""
function gel_jump_problem(model::GMMModel, h::Function=log)
  error("gel_jump_problem not defined for model type $(typeof(model))")
end

"""
    gel_nlp_problem(model::GMMModel, h::Function=log)

Returns NLPModel problem for generalized empirical likelihood estimation of `model`.
`h` is a generalized likelihood function.
"""
function gel_nlp_problem(model::GMMModel, h::Function=log)
  gi = get_gi(model)
  n = number_observations(model)
  d = number_parameters(model)
  k = number_moments(model)
  el(x) = -sum(h.(x[1:n]))
  c = gmm_constraints(model)
  nc = isnothing(c) ? 0 : length(c(zeros(d)))
  con = function(x)
    p = x[1:n]
    θ = x[(n+1):end]
    if isnothing(c)
      out=[sum(p) sum(gi(θ).*p, dims=1)]
    else
      out=[sum(p) sum(gi(θ).*p, dims=1) c(θ)]
    end
    return(vec(out))
  end
  lcon = [1.0, zeros(k+nc)...]
  ucon = [1.0, zeros(k+nc)...]
  lvar = [zeros(n)... , fill(-1e10,d)...]
  uvar = [ones(n)... ,  fill(1e10,d)...]
  name = "GEL for $(typeof(model))"
  x0 = [fill(1/n,n)... , ones(d)...]
  return(ADNLPModel(el,x0, lvar, uvar,
                    con, lcon, ucon,
                    name=name))
end

"""
    gmm_nlp_problem(model::GMMModel, obj=gmm_objective)

Constructs NLPModel for GMMModel.
"""
function gmm_nlp_problem(model::GMMModel, obj=gmm_objective)
  con = gmm_constraints(model)
  if con===nothing
    return(ADNLPModel(obj(model), ones(number_parameters(model)),
                      name="$(typeof(model))"))
  else
    θ = zeros(number_parameters(model))
    nc = length(con(θ))
    return(ADNLPModel(obj(model), ones(number_parameters(model)),
                      c=con, lcon=zeros(nc), ucon=zeros(nc),
                      name="$(typeof(model))"))
  end
end

"""
    gmm_jump_problem(model::GMMModel, obj=gmm_objective)

Constructs JuMP problem for GMMModel.
"""
function gmm_jump_problem(model::GMMModel, obj=gmm_objective)
  con = gmm_constraints(model)
  m = Model()
  d = number_parameters(model)
  @variable(m, θ[1:d])
  if con!==nothing
    error("gmm_jump_problem not implemented with constraints for $(typeof(GMMModel))")
  end
  f = obj(model)
  JuMP.register(m, :obj, d, (θ...)->f([θ...]), autodiff=true)
  @NLobjective(m, Min, obj(θ...))
  return(m)
end

"""
    gel_optim_args(model::GMMModel, h::Function=log)

Return tuple, `out` for calling `optimize(out..., IPNewton())` afterward.

It seems that IPNewton() works better if the constraint on p is 0 ≤
sum(p) ≤ 1 instead of sum(p)=1, and you begin the optimizer from a
point with sum(p) < 1.
"""
function gel_optim_args(model::GMMModel, gel::Function=log)
  gi = get_gi(model)
  n = number_observations(model)
  d = number_parameters(model)
  k = number_moments(model)
  el(x) = -sum(gel.(x[1:n]))
  ∇el!(g,x) = (g.= 0; g[1:n] .= -(p->ForwardDiff.derivative(gel,p)).(x[1:n]); g)
  ∇²el! = function (h,x)
    h .= 0
    for i in 1:n
      h[i,i] = -ForwardDiff.derivative(p->ForwardDiff.derivative(gel,p),x[i])
    end
    h
  end
  cfn = gmm_constraints(model)
  nc = cfn===nothing ? 0 : length(cfn(zeros(d)))
  con_c! = function(c,x)
    p = x[1:n]
    θ = x[(n+1):end]
    c[1] = sum(p)
    c[2:(k+1)] .= sum(gi(θ).*p, dims=1)[:]
    if cfn!==nothing
      c[(k+2):(k+nc+1)] .= cfn(θ)
    end
    c
  end
  con_jac! = function(J,x)
    p = x[1:n]
    θ = x[(n+1):end]
    J[1,1:n] .= 1
    J[1,(n+1):end] .= 0
    J[2:(k+1),1:n] .= gi(θ)'
    J[2:(k+1),(n+1):end] .= ForwardDiff.jacobian(t->sum(gi(t).*p, dims=1), θ)
    if (cfn!==nothing)
      J[(k+2):(k+nc+1),1:n].=0
      J[(k+2):(k+nc+1),(n+1):end].=ForwardDiff.jacobian(cfn, θ)
    end
    J
  end
  con_hess! = function(h,x,λ)
    p = x[1:n]
    θ = x[(n+1):end]
    λDgi = ForwardDiff.jacobian(x->gi(x)*λ[2:(1+k)], θ)
    for i in 1:n
      for m in 1:d
        h[i,(n+m)] += λDgi[i,m]
        h[(n+m),i] += λDgi[i,m]
      end
    end
    h[(n+1):(n+d),(n+1):(n+d)] .+=
      ForwardDiff.hessian(t->dot(λ[2:(1+k)],sum(gi(t).*p, dims=1)), θ)
    if (cfn!==nothing)
      h[(n+1):(n+d),(n+1):(n+d)] .+= ForwardDiff.hessian(t->dot(λ[(2+k):end],cfn(t)), θ)
    end
    h
  end

  lcon = [0.0, zeros(k+nc)...]
  ucon = [1.0, zeros(k+nc)...]
  lvar = [zeros(n)... , fill(-Inf,d)...]
  uvar = [ones(n)... ,  fill(Inf,d)...]
  x0 = [fill(1/(1.1*n),n)... , zeros(d)...]
  objective = TwiceDifferentiable(el, ∇el!, ∇²el!, x0)
  constraints = TwiceDifferentiableConstraints(con_c!, con_jac!, con_hess!,
                                               lvar, uvar, lcon, ucon)
  return((objective=objective,constraints=constraints,x0=x0))
end


text:
# Copyright (c) 2020 Idiap Research Institute, http://www.idiap.ch/
#  Niccolò Antonello <nantonel@idiap.ch>

# testing underflow not occuring single precision
# full
T = Float32
Nt = 21
Ns = 10
a = rand(T,Ns)   # initial state probability
a ./= norm(a,1)
A = rand(T,Ns,Ns)
for i in 1:Ns A[i,:] ./= norm(A[i,:],1) end
@assert HMMGradients.is_transition_matrix(A)

y = rand(T,Nt,Ns)
alpha, c = forward(Nt,a,A,y)
beta = backward(Nt,A,c,y)
gamma = posterior(Nt,a,A,y)
cost = nlogML(Nt,a,A,y)

grada = zeros(T,Ns)
gradA = zeros(T,Ns,Ns)
grady = zeros(T,Nt,Ns)
betat1 = ones(T,Ns)
betat0 = zeros(T,Ns)
grada, gradA, grady = HMMGradients.get_grad!(grada,gradA,grady,betat0,betat1,Nt,a,A,y,alpha,c)
grada2, gradA2, grady2 = HMMGradients.get_grad(Nt,a,A,y,alpha,c)
A, a, y = Float64.(A), Float64.(a), Float64.(y) # needs double precision
grada_fd = gradient_fd(a -> nlogML(Nt,a,A,y),a)
gradA_fd = gradient_fd(A -> nlogML(Nt,a,A,y),A)
grady_fd = gradient_fd(y -> nlogML(Nt,a,A,y),y)
grada_an = -1/(Nt*c[1]) .* (Diagonal(y[1,:]) * beta[1,:])  # analyatic formula
grada_an2 = -1/(Nt) .* gamma[1,:] ./ a  # analyatic formula
gradA_an = -1/(Nt) * sum( alpha[t,:]*(Diagonal(1/c[t+1] .* y[t+1,:])*beta[t+1,:])' for t = 1:Nt-1)  # analyatic formula
gradA_an2 = -1/(Nt) * sum( alpha[t,:]*(A\beta[t,:])' for t = 1:Nt-1)  # analyatic formula
grady_an = -1/Nt .* gamma ./ y # analyatic formula
grady_an2 = similar(grady_an)
grady_an2[1,:] = -1/(c[1]*Nt) .* (Diagonal(a)*beta[1,:])
for t = 2:Nt
  grady_an2[t,:] .= -1/(c[t]*Nt) .* Diagonal(A'*alpha[t-1,:])*beta[t,:] # analyatic formula
end
fy,p = ChainRulesCore.rrule(nlogML, Nt, a, A, y)
_,_,grada3,gradA3,grady3 = p(1.0)

if iseven(Nt)
  @test norm(beta[1,:]-betat0) < 1e-5
else
  @test norm(beta[1,:]-betat1) < 1e-5
end
@test norm(grada  - grada_fd, Inf)/norm(grada_fd) < 1e-4
@test norm(gradA  - gradA_fd, Inf)/norm(gradA_fd) < 1e-4
@test norm(grady  - grady_fd, Inf)/norm(grady_fd) < 1e-4
@test norm(grada2  - grada_fd, Inf)/norm(grada_fd) < 1e-4
@test norm(gradA2  - gradA_fd, Inf)/norm(gradA_fd) < 1e-4
@test norm(grady2  - grady_fd, Inf)/norm(grady_fd) < 1e-4
@test norm(grada3  - grada_fd, Inf)/norm(grada_fd) < 1e-4
@test norm(gradA3  - gradA_fd, Inf)/norm(gradA_fd) < 1e-4
@test norm(grady3  - grady_fd, Inf)/norm(grady_fd) < 1e-4
@test norm(grada_an   - grada_fd, Inf)/norm(grady_fd) < 1e-4
@test norm(grada_an2  - grada_fd, Inf)/norm(grady_fd) < 1e-4
@test norm(gradA_an   - gradA_fd, Inf)/norm(grady_fd) < 1e-4
@test norm(gradA_an2  - gradA_fd, Inf)/norm(grady_fd) < 1e-4
@test norm(grady_an   - grady_fd, Inf)/norm(grady_fd) < 1e-4
@test norm(grady_an2  - grady_fd, Inf)/norm(grady_fd) < 1e-4

# test with zeropadded
Nz = 20
y2 = [y;zeros(eltype(y),Nz,size(y,2))]
alpha, c = forward(Nt,a,A,y2)
grada,gradA,grady = HMMGradients.get_grad(Nt,a,A,y2,alpha,c)
gradA_fd = gradient_fd(A->nlogML(Nt,a,A,y),Float64.(A))
grady_fd = gradient_fd(y->nlogML(Nt,a,A,y),Float64.(y2))

@test size(grady) == size(y2)
@test norm(grada  - grada_fd, Inf)/norm(grada_fd) < 1e-4
@test norm(gradA  - gradA_fd, Inf)/norm(gradA_fd) < 1e-4
@test norm( grady  - grady_fd, Inf)/norm(grady_fd) < 1e-5
@test sum(grady[end-(Nz-1):end,:]) == 0.0

## test with batch Nt × Ns × Nb
T = Float32
Nb = 3
Ns = 10
Nts = [10,15,20]
A  = rand(T,Ns,Ns)
As = zeros(T,Ns,Ns,Nb)
for i = 1:Nb As[:,:,i] .= A end
as = rand(T,Ns,Nb)
ys = rand(T,maximum(Nts),Ns,Nb)
alphas,cs = forward(Nts,as,As,ys)
alphas_cs = [forward(Nts[i],as[:,i],As[:,:,i],ys[:,:,i]) for i=1:Nb]
@test all([norm(alphas[i]-alphas_cs[i][1])<1e-5 for i=1:Nb] )
@test all([norm(cs[i]-alphas_cs[i][2])<1e-5 for i=1:Nb] )
@test all(length(cs) .== Nb)
@test all(size.(alphas,2) .== Ns)
@test all(size.(alphas,1) .== Nts)
grada, gradA, grady = HMMGradients.get_grad(Nts,as,As,ys,alphas,cs)
grada_gradA_grady = [HMMGradients.get_grad(Nts[i],as[:,i],As[:,:,i],ys[:,:,i],alphas[i],cs[i]) for i=1:Nb]
fys = [nlogML(Nts[i],as[:,i],As[:,:,i],ys[:,:,i]) for i=1:Nb]
fy2 = nlogML(Nts,as,As,ys)
fy,p = ChainRulesCore.rrule(nlogML, Nts, as, As, ys)
_,_,grada2,gradA2,grady2 = p(1.0)

@test sum(fys) ≈ fy
@test fy2 == fy
@test all([norm(grada[:,i]-grada_gradA_grady[i][1]) for i = 1:Nb] .< 1e-5)
@test all([norm(gradA[:,:,i]-grada_gradA_grady[i][2]) for i = 1:Nb] .< 1e-5)
@test all([norm(grady[:,:,i]-grada_gradA_grady[i][3]) for i = 1:Nb] .< 1e-5)
@test all([norm(grada2[:,i]-grada_gradA_grady[i][1]) for i = 1:Nb] .< 1e-5)
@test all([norm(gradA2[:,:,i]-grada_gradA_grady[i][2]) for i = 1:Nb] .< 1e-5)
@test all([norm(grady2[:,:,i]-grada_gradA_grady[i][3]) for i = 1:Nb] .< 1e-5)

# constrained path
T = Float32
A = T[0.0 1.0 0.0; 0.0 0.5 0.5; 1.0 0.0 0.0]
t2tr = Dict{Int,Array{Int,1}}[Dict(1 => [2]), Dict(2 => [3, 2]), Dict(2 => [3, 2],3 => [1]), Dict(2 => [3, 2],3 => [1],1 => [2]), Dict(2 => [3, 2],3 => [1],1 => [2]), Dict(2 => [3, 2],3 => [1],1 => [2]), Dict(2 => [2],3 => [1],1 => [2]), Dict(2 => [2],1 => [2]), Dict(2 => [3])]
t2IJ= HMMGradients.t2tr2t2IJ(t2tr)

Nt = length(t2tr)+1
Ns = size(A,1)
y = rand(T,Nt,Ns)

alpha, c = forward(Nt,t2IJ,A,y)
@test all( sum(alpha,dims=2) .≈ 1.0)

beta = backward(Nt,A,c,t2IJ,y)

gamma = posterior(Nt,t2IJ,A,y)
@test all( sum(gamma,dims=2) .≈ 1.0 )

grada = zero(T)
gradA = zero(T)
grady = zeros(T,Nt,Ns)
betat1 = ones(T,Ns)
betat0 = zeros(T,Ns)
_,_,grady  = HMMGradients.get_grad!(grada,gradA,grady,betat0,betat1,Nt,t2IJ,A,y,alpha,c)
_,_,grady2 = HMMGradients.get_grad(Nt,t2IJ,A,y,alpha,c)

fy,p = ChainRulesCore.rrule(nlogML, Nt, t2IJ, A, y)
_,_,grada3,gradA3,grady3 = p(1.0)

A,y = Float64.(A), Float64.(y) # needs double precision
grady_fd = gradient_fd(y->nlogML(Nt,t2IJ,A,y),y)

@test norm(beta[1,:]-betat0) < 1e-5
@test norm(grady  - grady_fd, Inf)/norm(grady_fd) < 1e-5
@test norm(grady2  - grady_fd, Inf)/norm(grady_fd) < 1e-5
@test norm(grady3  - grady_fd, Inf)/norm(grady_fd) < 1e-5
@test gradA3 == ChainRulesCore.NoTangent()
@test grada3 == ChainRulesCore.NoTangent()

# test with zeropadded
Nz = 20
y2 = [y;zeros(eltype(y),Nz,size(y,2))]
alpha, c = forward(Nt,t2IJ,A,y2)
_,_,grady = HMMGradients.get_grad(Nt,t2IJ,A,y2,alpha,c)
grady_fd = gradient_fd(y->nlogML(Nt,t2IJ,A,y),Float64.(y2)) # needs double precision

@test size(grady) == size(y2)
@test norm( grady  - grady_fd, Inf)/norm(grady_fd) < 1e-5
@test sum(grady[end-(Nz-1):end,:]) == 0.0

## test with batch Nt × Ns × Nb
T = Float32
Nb = 3
A = T[0.0 1.0 0.0; 0.0 0.5 0.5; 1.0 0.0 0.0]
t2IJs = [t2IJ for i=1:Nb]
Nts = length.(t2IJs).+1
y = rand(T,maximum(Nts),Ns,Nb)
alphas,cs = forward(Nts,t2IJs,A,y)
alphas_cs = [forward(Nts[i],t2IJs[i],A,y[:,:,i]) for i=1:Nb]
@test all([norm(alphas[i]-alphas_cs[i][1])<1e-5 for i=1:Nb] )
@test all([norm(cs[i]-alphas_cs[i][2])<1e-5 for i=1:Nb] )
@test all(length(cs) .== Nb)
@test all(size.(alphas,2) .== Ns)
@test all(size.(alphas,1) .== Nts)
grada, gradA, grady = HMMGradients.get_grad(Nts,t2IJs,A,y,alphas,cs)
grada_gradA_grady = [HMMGradients.get_grad(Nts[i],t2IJs[i],A,y[:,:,i],alphas[i],cs[i]) for i=1:Nb]
fy,p = ChainRulesCore.rrule(nlogML, Nts, t2IJs, A, y)
fy2 = nlogML(Nts,t2IJs,A,y)
_,_,grada2,gradA2,grady2 = p(1.0)

@test fy == fy2
@test all([norm(grady[:,:,i]-grada_gradA_grady[i][3]) for i = 1:Nb] .< 1e-5)
@test all([norm(grady2[:,:,i]-grada_gradA_grady[i][3]) for i = 1:Nb] .< 1e-5)
@test gradA2 == ChainRulesCore.NoTangent()
@test grada2 == ChainRulesCore.NoTangent()


text:
# Creates `PGFPlotsX.Coordinates` from hits
function ZTCoordinates(hits::Vector{CalibratedHit}; offset=0)
    PGFPlotsX.Coordinates(map(h->h.t - offset, hits), map(h->h.pos.z, hits))
end

"""
    function ztplot(hits::Vector{CalibratedHit}; du=nothing, t₀=0) -> PGFPlotsX.Axis

Create a zt-plot from a set of hits using PGFPlotsX. If no `du` is provided,
the brightest DU will be determined.
"""
function ztplot(hits::Vector{CalibratedHit}; du=nothing, t₀=0)
    if du == nothing
        du = NeRCA.most_frequent(h -> h.du, hits)
        hits = filter(h -> h.du == du, hits)
    end
    sort!(hits, by=h->h.t)
    triggered_hits = filter(h -> h.triggered, hits)
    ax = @pgf Axis({xlabel=raw"time [\si{\ns}]", ylabel=raw"z [\si{\m}]"})
    push!(ax,  @pgf PlotInc({only_marks, mark="*", mark_size=".8", thick}, ZTCoordinates(hits, offset=t₀)))
    push!(ax, LegendEntry("hits"))
    push!(ax,  @pgf PlotInc({only_marks, mark="+", mark_size="5", thick}, ZTCoordinates(triggered_hits, offset=t₀)))
    push!(ax, LegendEntry("triggered"))
    ax
end

"""
    function ztplot(hits::Vector{CalibratedHit}, fit::NeRCA.ROyFit; t₀=0) -> PGFPlotsX.Axis

Create a zt-plot from hits and a ROyFit using PGFPlotsX.
"""
function ztplot(hits::Vector{CalibratedHit}, fit::NeRCA.ROyFit; t₀=0)
    ax = ztplot(hits; t₀=t₀)
    zs = range(0, maximum(h->h.pos.z, hits) * 1.1, length=200)
    dᵧ, ccalc = NeRCA.make_cherenkov_calculator(fit.sdp)
    push!(ax, @pgf PlotInc({ultra_thick}, Coordinates(ccalc.(zs) .- t₀, zs)))
    push!(ax, LegendEntry("fit"))
    push!(ax, @pgf PlotInc(
        {thick, only_marks, mark="o", mark_size="4"},
        ZTCoordinates(fit.selected_hits, offset=t₀))
    )
    push!(ax, LegendEntry("selected"))
    ax
end

"""
    f(dom_positions, track::NeRCA.Track)

Plot recipe for z-t-plots.
"""
@recipe function f(dom_positions, track::NeRCA.Track)
    # layout := @layout [a b]
    layout := (2, 2)

    seriestype := :scatter

    xlabel := "time [ns]"
    ylabel := "z [m]"

    @series begin
        subplot := 1
        ccalc = NeRCA.make_cherenkov_calculator(track)
        ccalc.(dom_positions), [p.z for p in dom_positions]
    end

    xlabel := "x [m]"
    ylabel := "y [m]"

    @series begin
        subplot := 2
        label := "DOMs"
        [p.x for p in dom_positions], [p.y for p in dom_positions]
        [track.pos.x], [track.pos.y]
    end

    @series begin
        subplot := 2
        label := "Track"
        [track.pos.x], [track.pos.y]
    end

    @series begin
        subplot := 3
    end
end

"""
    f(hits::Vector{CalibratedHit}; label="hits", highlight_triggered=false, multiplicities=false, Δt=20, du=0, t₀=0)

Plot recipe to plot simple z-t-plots.
"""
@recipe function f(hits::Vector{CalibratedHit}; label="hits", markersize=4, highlight_triggered=false, multiplicities=false, Δt=20, du=0, t₀=0)
    seriestype := :scatter

    xlabel := "time [ns]"
    ylabel := "z [m]"
    markerstrokewidth := 0

    if du > 0
        hits = filter(h -> h.du == du, hits)
    end

    thits = filter(h -> h.triggered, hits)

    @series begin
        label := label
        if multiplicities
            markersize := count_multiplicities(hits)[1]
            markeralpha := 0.8
        end
        [h.t - t₀ for h in hits], [h.pos.z for h in hits]
    end

    if highlight_triggered
        @series begin
            label := "triggered"
            markersize := 1
            marker := :x
            markerstrokewidth := 1
            [h.t - t₀ for h in thits], [h.pos.z for h in thits]
        end
    end
end

@recipe function f(hits::Vector{CalibratedHit}, track::Track; max_z=nothing)
    seriestype := :scatter

    xlabel := "time [ns]"
    ylabel := "z [m]"
    background_color_legend := PlotThemes.RGBA{Float64}(1.0,1.0,1.0,0.4)
    markerstrokewidth := 0

    thits = filter(h -> h.triggered, hits)
    if max_z == nothing
        max_z = maximum(map(h->h.pos.z, hits))
    end

    dus = sort(unique(map(h->h.du, thits)))
    ccalc = make_cherenkov_calculator(track, v=norm(track.dir)*1e9)

    colours = palette(:default)
    for (idx, du) in enumerate(dus)
        du_hits = filter(h -> h.du == du, hits)
        @series begin
            label := "DU $(du)"
            linewidth := 3
            markersize := 5
            markercolor := colours[idx]
            du_hits
        end
        x = du_hits[1].pos.x
        y = du_hits[1].pos.y
        zs = range(0, max_z, length=200)
        @series begin
            linewidth := 10
            label := ""
            markersize := 1
            markercolor := colours[idx]
            [ccalc(Position(x, y, z)) for z in zs], zs
        end
    end
end

@recipe function f(hits::Vector{CalibratedHit}, fit::ROyFit; label="", max_z=nothing)
    seriestype := :scatter

    xlabel := "time [ns]"
    ylabel := "z [m]"
    background_color_legend := PlotThemes.RGBA{Float64}(1.0,1.0,1.0,0.4)
    markerstrokewidth := 0

    triggered_hits = filter(h -> h.triggered, hits)
    if max_z == nothing
        max_z = maximum(map(h->h.pos.z, hits))
    end

    dᵧ, ccalc = make_cherenkov_calculator(fit.sdp)

    @series begin
        markersize := 2
        marker := :circle
        markeralpha := 0.7
        label := label
        hits
    end
    @series begin
        markersize := 5
        markeralpha := 0.4
        fillcolor := nothing
        marker := :circle
        label := ""
        triggered_hits
    end
    @series begin
        markersize := 10
        marker := :cross
        label := ""
        fit.selected_hits
    end
    zs = range(0, max_z, length=200)
    @series begin
        linewidth := 10
        label := ""
        markersize := 1
        label := ""
        ccalc.(zs), zs
    end
end


text:
# ------------------------------------------------------------------
# Licensed under the ISC License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    LearningProblem(sdata, tdata, task)

A spatial learning problem with source data `sdata`,
target data `tdata`, and learning `task`.

## Examples

Create a clustering problem based on a set of soil features:

```julia
julia> LearningProblem(sdata, tdata,
                       ClusteringTask((:moisture,:mineral)))
```
"""
struct LearningProblem{DΩₛ<:AbstractData,
                       DΩₜ<:AbstractData,
                       T<:AbstractLearningTask} <: AbstractProblem
  sdata::DΩₛ
  tdata::DΩₜ
  task::T

  function LearningProblem{DΩₛ,DΩₜ,T}(sdata, tdata, task) where {DΩₛ<:AbstractData,
                                                                 DΩₜ<:AbstractData,
                                                                 T<:AbstractLearningTask}
    sourcevars = keys(variables(sdata))
    targetvars = keys(variables(tdata))

    # assert spatial configuration
    @assert ndims(sdata) == ndims(tdata) "source and target data must have the same number of dimensions"
    @assert coordtype(sdata) == coordtype(tdata) "source and target data must have the same coordinate type"

    # assert task is compatible with the data
    if iscomposite(task)
      # TODO
    else
      @assert features(task) ⊆ sourcevars "features must be present in source data"
      @assert features(task) ⊆ targetvars "features must be present in target data"
      if issupervised(task)
        @assert label(task) ∈ sourcevars "label must be present in source data"
      end
    end

    new(sdata, tdata, task)
  end
end

LearningProblem(sdata::DΩₛ, tdata::DΩₜ, task::T) where {DΩₛ<:AbstractData,
                                                        DΩₜ<:AbstractData,
                                                        T<:AbstractLearningTask} =
  LearningProblem{DΩₛ,DΩₜ,T}(sdata, tdata, task)

"""
    sourcedata(problem)

Return the source data of the learning `problem`.
"""
sourcedata(problem::LearningProblem) = problem.sdata

"""
    targetdata(problem)

Return the target data of the learning `problem`.
"""
targetdata(problem::LearningProblem) = problem.tdata

"""
    task(problem)

Return the learning task of the learning `problem`.
"""
task(problem::LearningProblem) = problem.task

# ------------
# IO methods
# ------------
function Base.show(io::IO, problem::LearningProblem)
  dim = ndims(problem.sdata)
  print(io, "$(dim)D LearningProblem")
end

function Base.show(io::IO, ::MIME"text/plain", problem::LearningProblem)
  println(io, problem)
  println(io, "  source: ", problem.sdata)
  println(io, "  target: ", problem.tdata)
  print(  io, "  task:   ", problem.task)
end


text:
using ECON627_2020
using Test

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


text:
# const priΛ = "Dir"
# const priΛ = "SDM"
# const priΛ = "SBM"

# const priλ = "Dir"
# const priλ = "SDM"
# const priλ = "SBM"
# const priλ = "SBMSDM" # SBM on single lag, SDM on higher orders

# const priQ = "Dir"
# const priQ = "SDM"
# const priQ = "SBM"

# const srnd = 22219
# const TT = 50
# const L = 6
# const R = 1
#
# whichdat = "sim1_sets"
# whichdat = "sim2_sets"
# whichdat = "pinkSalmon_first30_K4"

### if ARGS used

const priΛ = ARGS[1]
const priλ = ARGS[2]
const priQ = ARGS[3]
const srnd = parse(Int64, ARGS[4])
const TT = parse(Int64, ARGS[5])
const L = parse(Int64, ARGS[6])
const R = parse(Int64, ARGS[7])
const whichdat = ARGS[8]

### after user selections:

using BSON
using StatsBase

# pkg> add "https://github.com/mheiner/SparseProbVec.jl"
using SparseProbVec

# pkg> add "https://github.com/mheiner/MTD.jl"
using MTD

using Random
using Dates

# srnd = 10091
Random.seed!(srnd)

### Read in data
simdata = occursin(r"sim\d+_sets", whichdat) # run validation at end?

if simdata
    BSON.@load "data/$(whichdat).bson" y X P
    const S = deepcopy(y[1])
else
    BSON.@load "data/$(whichdat).bson" S
end

### Setup
const SS = S[1:TT]
println(sort(unique(SS[(L+1):TT])))
println(counts(SS))
const K = maximum(S)
const λ_indx = build_λ_indx(L, R)

const n_keep = Int(2e3)
const n_burn = Int(200e3)
const thin = 200

const mesg = "$(whichdat)_MMTD_T$(TT)_L$(L)_R$(R)_K$(K)_$(priΛ)Lam_$(priλ)lam_$(priQ)Qmarg_$(srnd)"

### Prior

fieldnames(PriorMMTD)
symm_prior = symmetricDirPrior_mmtd(1.0, 1.0, 1.0, 1.0, L, R, K, λ_indx)

η_Λ = 1.0e3
π_small_Λ = vcat(0.0, fill(0.25, R-1))
# π_small_Λ = vcat(0.0, fill(1.0/float(R+1), R-1))
π_large_Λ = π_small_Λ .* 1.0
α_Λ = [ 5.0 / 2.0^float(r-1) for r = 1:(R+1) ] # decreasing alpha when mimicking Dirichlet; currently only used in initialization
# α_Λ = vcat(float(R), fill(1.0/float(R), R))
# γ_Λ, δ_Λ = shape_Dir2GenDir(α_Λ)
γ_Λ = 1.0
δ_Λ = 1.0
# δc_Λ = del_correctionSBM(π_small_Λ[2], π_large_Λ[2], L)
δc_Λ = 1.0

# β_λ = log(float(TT))
β_λ = sqrt(float(TT))
# β_λ = float(TT) / 10.0
# β_λ = 1.001

η_λ = 1.0e3
π_small_λ = 0.75
π_large_λ = 0.10
γ_λ, δ_λ = shape_Dir2GenDir(symm_prior[2][1]) # this works for SBMSDM only (will need one for each set otherwise)
δc_λ = del_correctionSBM(π_small_λ, π_large_λ, L)
# γ_λ = 1.5
# δ_λ = 1.5
# δc_λ = 1.0

# β_Q = float(TT) / 8.0 # for SDM on Q
β_Q = float(K) / 1.0

η_Q = 1.0e3
π_small_Q = 1.0 / float(K)
π_large_Q = π_small_Q / 5.0
γ_Q, δ_Q = shape_Dir2GenDir(symm_prior[4][1][:,1])
δc_Q = del_correctionSBM(π_small_Q, π_large_Q, K)

if priΛ == "Dir"
        prior_Λ = deepcopy(symm_prior[1]) # Dirichlet
    elseif priΛ == "DirLowM"
        prior_Λ = deepcopy(α_Λ) # Dirichlet with decreasing α
    elseif priΛ == "SDM"
        prior_Λ = SparseDirMix(symm_prior[1], β_Λ)
    elseif priΛ == "SBM"
        prior_Λ = SBMprior(R+1, η_Λ, π_small_Λ, π_large_Λ, γ_Λ, δc_Λ.*δ_Λ)
end

if priλ == "Dir"
        prior_λ = deepcopy(symm_prior[2])
    elseif priλ == "SDM"
        prior_λ = [ SparseDirMix(symm_prior[2][r], β_λ) for r in 1:R ]
    elseif priλ == "SBM"
        prior_λ = [ SBMprior(λ_indx.lens[r], η_λ, π_small_λ, π_large_λ, γ_λ, δc_λ.*δ_λ) for r in 1:R ]
    elseif priλ == "SBMSDM"
        prior_λ = Vector{Union{SparseDirMix, SBMprior}}(undef, R)
        prior_λ[1] = SBMprior(λ_indx.lens[1], η_λ, π_small_λ, π_large_λ, γ_λ, δc_λ.*δ_λ)
        for r in 2:R
            prior_λ[r] = SparseDirMix(symm_prior[2][r], β_λ)
        end
end

prior_Q0 = deepcopy(symm_prior[3])

if priQ == "Dir"
        prior_Q = deepcopy(symm_prior[4]) # for Dirichlet
    elseif priQ == "SDM"
        prior_Q = [ [ SparseDirMix(symm_prior[4][1][:,1], β_Q) for kk in 1:(K^r) ] for r in 1:R ] # SDM
        prior_Q = [ reshape(prior_Q[r], fill(K, r)... ) for r in 1:R ]
    elseif priQ == "SBM"
        prior_Q = [ [ SBMprior(K, η_Q, π_small_Q, π_large_Q, γ_Q, δc_Q.*δ_Q) for kk in 1:(K^r) ] for r in 1:R ] # SBM
        prior_Q = [ reshape(prior_Q[r], fill(K, r)... ) for r in 1:R ]
end

bfact = bfact_MC(SS, L, R, K, prior_Q0, prior_Q)

using Printf
bfact_file = open("BayesFactors/bfact_$(mesg).txt", "a+")
write(bfact_file, "order 0, lags 0: $(@sprintf("%.2e", bfact[2][1,3])), llik: $(@sprintf("%.2e", bfact[2][1,4])) \n")
for i in 1:bfact[1].nZζ
    Znow, ζnow = bfact[1].Zζindx[i,:]
    write(bfact_file, "order $(Znow), lags $(bfact[1].indxs[Znow][ζnow]): $(@sprintf("%.2e", bfact[2][i+1,3])), llik: $(@sprintf("%.2e", bfact[2][i+1,4])) \n")
end
close(bfact_file)

### Inits
fieldnames(ParamsMMTD)
inits = ParamsMMTD(
  SparseProbVec.rDirichlet(α_Λ, logout=true), # Λ
  [ SparseProbVec.rDirichlet(symm_prior[2][r], logout=true) for r in 1:R], # λ
  [ StatsBase.sample(1:(λ_indx.nZζ+1)) - 1 for i in 1:(TT-L) ], # Zζ
  SparseProbVec.rDirichlet(prior_Q0, logout=true),
  # [ log(fill(1.0/K, fill(K, r+1)...)) for r in 1:R ] # lQ (const)
  [ reshape(vcat([ SparseProbVec.rDirichlet(ones(K), logout=true) for k in 1:K^(r) ]...), fill(K, r+1)... ) for r in 1:R ] # lQ (rand)
  # [ log( transTens_MLE( count_trans_L(SS, K, r) + 1.0 ) ) for r in 1:R ] # lQ (mle)
  )

### Create model
fieldnames(ModMMTD)
model = ModMMTD(L, R, K, TT, SS,
  PriorMMTD(prior_Λ, prior_λ, prior_Q0, prior_Q),
  inits, λ_indx)

logfilename = "postsim_progress/out_prog_$(mesg).txt"
timestart = Dates.now()

# timing for benchmarks
timemod!(5, model, 100, logfilename)

## Burn in
mcmc!(model, Int(0.15*n_burn)-500, save=false, report_filename=logfilename,
      thin=1, jmpstart_iter=5, report_freq=1000)
mcmc!(model, Int(0.85*n_burn), save=false, report_filename=logfilename,
      thin=1, jmpstart_iter=10, report_freq=10000)

## Estimate time remaining
etr(timestart, n_keep, thin, logfilename)

### MCMC
sims = mcmc!(model, n_keep, save=true, report_filename=logfilename,
             thin=thin, jmpstart_iter=10)

### Output
bson("./postsim/mcmc_$(mesg).bson", sims=deepcopy(sims),
    model=deepcopy(model), whichdat=deepcopy(whichdat),
    mesg=deepcopy(mesg))

addDecomp_postsims!(sims, model) # saved before adding decompositions just in case it encounters an error

bson("./postsim/mcmc_$(mesg).bson", sims=deepcopy(sims),
    model=deepcopy(model), whichdat=deepcopy(whichdat),
    mesg=deepcopy(mesg))

### Send simulations to R
using RCall

nsim = length(sims)

R"rm(list=ls())"
priorinfo = "$(model.prior)"
sims_llik = [ sims[i][:llik] for i=1:nsim ]

lamindx = model.λ_indx.indxs
lamlens = model.λ_indx.lens
lamZzetaindx = model.λ_indx.Zζindx
lamNZzeta = model.λ_indx.nZζ

sims_Lam = permutedims(hcat([ exp.(sims[i][:lΛ]) for i=1:nsim ]...))
sims_lam = [ permutedims(hcat([ exp.(sims[i][:lλ][r]) for i=1:nsim]...)) for r in 1:R ]
sims_Q0 = permutedims(hcat([ exp.(sims[i][:lQ0]) for i=1:nsim ]...))
sims_Q = [ permutedims(hcat([ exp.(vec(sims[i][:lQ][r])) for i=1:nsim]...)) for r in 1:R ]

sims_Dlevelwgt = permutedims(hcat([ sims[ii][:Dlevelweight] for ii in 1:nsim ]...)) # decomposition level weights
sims_Dlaginvolv = permutedims(hcat([ sims[ii][:Dlaginvolv] for ii in 1:nsim ]...)) # decomposition lag weights

@rput sims_Lam sims_lam sims_Q0 sims_Q sims_llik priorinfo K L R TT lamindx lamlens lamZzetaindx lamNZzeta mesg sims_Dlevelwgt sims_Dlaginvolv

R"ls()"
R"save.image(file=paste0('postsim/mcmc_', $(mesg), '.rda'))"

### GoF
if simdata

    using Random
    using Dates
    using Statistics

    modeltype = "MMTD"
    report_filename = "forecast_validation/gof_$(mesg).txt"
    test_indx = 2

    report_file = open(report_filename, "a+")

    nprime = length(y[test_indx])
    niter = length(sims)
    nsim = 2000

    Random.seed!(srnd)
    simind = sort(StatsBase.sample(1:niter, nsim, replace=false))

    (FL1, MC, NLL, SQE, PL1, FOREC) = meanForecLoss(y[test_indx], X[test_indx][:,1:model.L],
        P[test_indx], nprime, lossL1, sims,
        model.TT, model.L, K,
        simind, λ_indx=model.λ_indx, modeltype=modeltype)

    PMforec = reshape(mean(FOREC, dims=1), (nprime, K) )
    PML1P = mean(abs.(PMforec .- P[test_indx])) # lower because of Jensen's inequality

    write(report_file, "L1 loss: $(round(100.0*mean(FL1), digits=2))\n" )
    write(report_file, "Misclass rate: $(round(100.0*mean(MC), digits=2))\n" )
    write(report_file, "NLL: $(round(mean(sum(NLL, dims=2)), digits=2))\n" )
    write(report_file, "SqErr: $(round(100.0*mean(SQE), digits=2))\n" )
    write(report_file, "postmean L1 loss, Ptrue: $(round(100.0*mean(PL1), digits=2))\n" )
    write(report_file, "L1 loss, Ptrue postmean forec: $(round(100.0*PML1P, digits=2))\n" )

    close(report_file)
end


text:
# syntax: proto3
using ProtoBuf
import ProtoBuf.meta
import ._ProtoBuf_Top_.tensorflow

mutable struct TensorMetadata <: ProtoType
    _type::Int32
    shape::tensorflow.TensorShapeProto
    TensorMetadata(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct TensorMetadata

mutable struct HostTransferMetadata <: ProtoType
    key::AbstractString
    metadata::Base.Vector{TensorMetadata}
    HostTransferMetadata(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct HostTransferMetadata

mutable struct HostComputeMetadata <: ProtoType
    device_to_host::Base.Vector{HostTransferMetadata}
    host_to_device::Base.Vector{HostTransferMetadata}
    HostComputeMetadata(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct HostComputeMetadata

export TensorMetadata, HostTransferMetadata, HostComputeMetadata


text:
mass(U, q, w, v, x, t, μ, η) = [1.0, 1.0, 1.0]
source(U, q, w, v, x, t, μ, η) = [0.0, 0.0, 0.0]
ubou(U, q, w, v, x, t, μ, η, û, n, τ) = [0.0, 0.0, 0.0]

function flux(U, q, w, v, x, t, μ, η)
    h, hu, hv = U

    h⁻¹ = 1/h

    uv = hu * h⁻¹ # u
    vv = hv * h⁻¹ # v

    g = μ[1] # gravitational acceleration

    p = g * h^2 / 2 # g * h^2 / 2

    f = [hu, hu * uv + p, hv * uv, hv, hu * vv, hv * vv + p]

    f = reshape(f, (3, 2))

    return f
end

function fbou(U, q, w, v, X, t, μ, η, Û, n, τ)
    f = flux(U, q, w, v, X, t, μ, η)
    f̂ = f[:, 1] * n[1] + f[:, 2] * n[2] + τ[1] * (U - Û)
    return f̂
end

function initu(X, μ, η)
    ϵ = 0.1 # Perturbation magnitude
    l = 0.5 # Gaussian width
    k = 0.5 # Sinusoidal wavenumber

    x = X[1]
    y = X[2]

    # The Bickley jet
    U = 1 / cosh(y)^2

    # Slightly off-center vortical perturbations
    ψ′ = exp(-(y + l/10)^2 / 2l^2) * cos(k * x) * cos(k * y)

    # Vortical velocity fields (ũ, ṽ) = (-∂_y, +∂_x) ψ̃
    u′ = + ψ′ * (k * tan(k * y) + y / l^2)
    v′ = - ψ′ * k * tan(k * x)

    Uᵢ₁ = 1.0        # h
    Uᵢ₂ = U + ϵ * u′ # h * u
    Uᵢ₃ = ϵ * v′     # h * v

    return [Uᵢ₁, Uᵢ₂, Uᵢ₃]
end

max_stars_repo_name: glwagner/Exasim

max_issues_repo_name: glwagner/Exasim

max_forks_repo_name: glwagner/Exasim


text:
"""
Description

      The last() method returns the last n number of observations in a TimeArray

Usage

      using MarketData
      last(cl)     # returns the last observation
      last(cl, 10) # returns the last 10 observations

Method Signature(s)

      last(ts::TimeSeries.TimeArray, n::Int=1)

Details

      The last function subsets a TimeSeries.Timearray by the last n observations. When n is not supplied, it
      defaults to n=1.

See Also

      last() which takes the last n observations from a TimeSeries.TimeArray

"""

function last(ts::TimeSeries.TimeArray, n::Int=1)
   ts[length(ts)-n+1:length(ts)]
end


text:
module LuckyImaging

using DataDeps
using Statistics
using SubpixelRegistration

export lucky_image,
       lucky_image!,
       testcube

include("util.jl")
include("classic.jl")
include("data.jl")
# include("fourier.jl")

function __init__()
    DataDeps.register(testcube_datadep)
end

end