text:
"""
OGRStyleMgr factory.

### Parameters
* `styletable`: OGRStyleTable or NULL if not working with a style table.

### Returns
an handle to the new style manager object.
"""
unsafe_createstylemanager(styletable = GDALStyleTable(C_NULL)) =
    StyleManager(GDAL.sm_create(styletable))

"""
Destroy Style Manager.

### Parameters
* `stylemanager`: handle to the style manager to destroy.
"""
destroy(sm::StyleManager) = (GDAL.destroy(sm.ptr); sm.ptr = C_NULL)

"""
Initialize style manager from the style string of a feature.

### Parameters
* `stylemanager`: handle to the style manager.
* `feature`: handle to the new feature from which to read the style.

### Returns
the style string read from the feature, or NULL in case of error.
"""
initialize!(stylemanager::StyleManager, feature::Feature) =
    GDAL.initfromfeature(stylemanager.ptr, feature.ptr)

"""
Initialize style manager from the style string.

### Parameters
* `stylemanager`: handle to the style manager.
* `stylestring`: the style string to use (can be NULL).
### Returns
TRUE on success, FALSE on errors.
"""
initialize!(stylemanager::StyleManager, stylestring = C_NULL) =
    Bool(GDAL.initstylestring(stylemanager.ptr, stylestring))

"""
Get the number of parts in a style.

### Parameters
* `stylemanager`: handle to the style manager.
* `stylestring`: (optional) the style string on which to operate. If NULL then
                 the current style string stored in the style manager is used.

### Returns
the number of parts (style tools) in the style.
"""
npart(stylemanager::StyleManager) = GDAL.getpartcount(stylemanager.ptr, C_NULL)

npart(stylemanager::StyleManager, stylestring::AbstractString) =
    GDAL.getpartcount(stylemanager.ptr, stylestring)

"""
Fetch a part (style tool) from the current style.

### Parameters
* `stylemanager`: handle to the style manager.
* `id`: the part number (0-based index).
* `stylestring`: (optional) the style string on which to operate. If not
    provided, then the current style string stored in the style manager is used.

### Returns
OGRStyleToolH of the requested part (style tools) or NULL on error.
"""
unsafe_getpart(stylemanager::StyleManager, id::Integer, stylestring = C_NULL) =
    StyleTool(GDAL.getpart(stylemanager.ptr, id, stylestring))

"""
Add a part (style tool) to the current style.

### Parameters
* `stylemanager`: handle to the style manager.
* `styletool`: the style tool defining the part to add.

### Returns
TRUE on success, FALSE on errors.
"""
addpart!(stylemanager::StyleManager, styletool::StyleTool) =
    Bool(GDAL.addpart(stylemanager.ptr, styletool.ptr))

"""
Add a style to the current style table.

### Parameters
* `stylemanager`: handle to the style manager.
* `stylename`: the name of the style to add.
* `stylestring`: (optional) the style string to use, or (if not provided) to use
    the style stored in the manager.

### Returns
TRUE on success, FALSE on errors.
"""
function addstyle!(
        stylemanager::StyleManager,
        stylename::AbstractString,
        stylestring::AbstractString
    )
    Bool(GDAL.addstyle(stylemanager.ptr, stylename, stylestring))
end

addstyle!(stylemanager::StyleManager, stylename::AbstractString) =
    Bool(GDAL.addstyle(stylemanager.ptr, stylename, C_NULL))

"""
OGRStyleTool factory.

### Parameters
* `classid`: subclass of style tool to create. One of OGRSTCPen (1),
             OGRSTCBrush (2), OGRSTCSymbol (3) or OGRSTCLabel (4).

### Returns
an handle to the new style tool object or NULL if the creation failed.
"""
unsafe_createstyletool(classid::OGRSTClassId) =
    StyleTool(GDAL.st_create(classid))

"""
Destroy Style Tool.

### Parameters
* `styletool`: handle to the style tool to destroy.
"""
destroy(styletool::StyleTool) =
    (GDAL.destroy(styletool.ptr); styletool.ptr = C_NULL)

"""
Determine type of Style Tool.

### Parameters
* `styletool`: handle to the style tool.

### Returns
the style tool type, one of OGRSTCPen (1), OGRSTCBrush (2), OGRSTCSymbol (3) or
OGRSTCLabel (4). Returns OGRSTCNone (0) if the OGRStyleToolH is invalid.
"""
gettype(styletool::StyleTool) = GDAL.gettype(styletool.ptr)

"""
Get Style Tool units.

### Parameters
* `styletool`: handle to the style tool.

### Returns
the style tool units.
"""
getunit(styletool::StyleTool) = GDAL.getunit(styletool.ptr)

"""
    OGR_ST_SetUnit(OGRStyleToolH styletool,
                   OGRSTUnitId eUnit,
                   double scale) -> void
Set Style Tool units.
### Parameters
* `styletool`: handle to the style tool.
* `newunit`: the new unit.
* `scale`: ground to paper scale factor.
"""
setunit!(styletool::StyleTool, newunit::OGRSTUnitId, scale::Real) =
    GDAL.setunit(styletool.ptr, newunit, scale)

"""
Get Style Tool parameter value as a string.

### Parameters
* `styletool`: handle to the style tool.
* `id`: the parameter id from the enumeration corresponding to the type of this
        style tool (one of the OGRSTPenParam, OGRSTBrushParam, OGRSTSymbolParam
        or OGRSTLabelParam enumerations)
* `nullflag`: pointer to an integer that will be set to TRUE or FALSE to
        indicate whether the parameter value is NULL.

### Returns
the parameter value as a string and sets `nullflag`.
"""
asstring(styletool::StyleTool, id::Integer, nullflag::Ref{Cint}) =
    GDAL.getparamstr(styletool.ptr, id, nullflag)

asstring(styletool::StyleTool, id::Integer) =
    asstring(styletool, id, Ref{Cint}(0))

"""
Get Style Tool parameter value as an integer.

### Parameters
* `styletool`: handle to the style tool.
* `id`: the parameter id from the enumeration corresponding to the type of this
        style tool (one of the OGRSTPenParam, OGRSTBrushParam, OGRSTSymbolParam
        or OGRSTLabelParam enumerations)
* `nullflag`: pointer to an integer that will be set to TRUE or FALSE to
        indicate whether the parameter value is NULL.

### Returns
the parameter value as an integer and sets `nullflag`.
"""
asint(styletool::StyleTool, id::Integer, nullflag::Ref{Cint}) =
    GDAL.getparamnum(styletool.ptr, id, nullflag)

asint(styletool::StyleTool, id::Integer) =
    asint(styletool, id, Ref{Cint}(0))

"""
Get Style Tool parameter value as a double.

### Parameters
* `styletool`: handle to the style tool.
* `id`: the parameter id from the enumeration corresponding to the type of this
        style tool (one of the OGRSTPenParam, OGRSTBrushParam, OGRSTSymbolParam
        or OGRSTLabelParam enumerations)
* `nullflag`: pointer to an integer that will be set to TRUE or FALSE to
        indicate whether the parameter value is NULL.

### Returns
the parameter value as a double and sets `nullflag`.
"""
asdouble(styletool::StyleTool, id::Integer, nullflag::Ref{Cint}) =
    GDAL.getparamdbl(styletool.ptr, id, nullflag)

asdouble(styletool::StyleTool, id::Integer) =
    asdouble(styletool, id, Ref{Cint}(0))

"""
Set Style Tool parameter value from a string.

### Parameters
* `styletool`: handle to the style tool.
* `id`: the parameter id from the enumeration corresponding to the type of this
        style tool (one of the OGRSTPenParam, OGRSTBrushParam, OGRSTSymbolParam
        or OGRSTLabelParam enumerations)
* `value`: the new parameter value
"""
setparam!(styletool::StyleTool, id::Integer, value::AbstractString) =
    GDAL.setparamstr(styletool.ptr, id, value)

"""
Set Style Tool parameter value from an integer.

### Parameters
* `styletool`: handle to the style tool.
* `id`: the parameter id from the enumeration corresponding to the type of this
        style tool (one of the OGRSTPenParam, OGRSTBrushParam, OGRSTSymbolParam
        or OGRSTLabelParam enumerations)
* `value`: the new parameter value
"""
setparam!(styletool::StyleTool, id::Integer, value::Integer) =
    GDAL.setparamnum(styletool.ptr, id, value)

"""
Set Style Tool parameter value from a double.

### Parameters
* `styletool`: handle to the style tool.
* `id`: the parameter id from the enumeration corresponding to the type of this
        style tool (one of the OGRSTPenParam, OGRSTBrushParam, OGRSTSymbolParam
        or OGRSTLabelParam enumerations)
* `value`: the new parameter value
"""
setparam!(styletool::StyleTool, id::Integer, value::Float64) =
    GDAL.setparamdbl(styletool.ptr, id, value)

"""
Get the style string for this Style Tool.

### Parameters
* `styletool`: handle to the style tool.

### Returns
the style string for this style tool or "" if the styletool is invalid.
"""
getstylestring(styletool::StyleTool) = GDAL.getstylestring(styletool.ptr)

"""
Return the r,g,b,a components of a color encoded in #RRGGBB[AA] format.

### Parameters
* `styletool`: handle to the style tool.
* `pszColor`: the color to parse

### Returns
(R,G,B,A) tuple of Cints.
"""
function getrgba(styletool::StyleTool, color::AbstractString)
    red = Ref{Cint}(0); green = Ref{Cint}(0)
    blue = Ref{Cint}(0); alpha = Ref{Cint}(0)
    result = Bool(GDAL.getrgbfromstring(styletool.ptr, color,
        red, green,
        blue, alpha
    ))
    result || error("Error in getting RGBA from Styletool")
    (red[], green[], blue[], alpha[])
end

"""
OGRStyleTable factory.

### Returns
an handle to the new style table object.
"""
unsafe_createstyletable() = StyleTable(GDAL.stbl_create())

"""
Destroy Style Table.

### Parameters
* `styletable`: handle to the style table to destroy.
"""
destroy(st::StyleTable) = (GDAL.destroy(st.ptr); st.ptr = C_NULL)

"""
Add a new style in the table.

### Parameters
* `styletable`: handle to the style table.
* `name`: the name the style to add.
* `stylestring`: the style string to add.

### Returns
TRUE on success, FALSE on error
"""
function addstyle!(
        styletable::StyleTable,
        stylename::AbstractString,
        stylestring::AbstractString
    )
    Bool(GDAL.addstyle(styletable.ptr, stylename, stylestring))
end

"""
Save a style table to a file.
### Parameters
* `styletable`: handle to the style table.
* `filename`: the name of the file to save to.
### Returns
TRUE on success, FALSE on error
"""
savestyletable(styletable::StyleTable, filename::AbstractString) =
    Bool(GDAL.savestyletable(styletable.ptr, filename))

"""
Load a style table from a file.

### Parameters
* `styletable`: handle to the style table.
* `filename`: the name of the file to load from.

### Returns
TRUE on success, FALSE on error
"""
loadstyletable!(styletable::StyleTable, filename::AbstractString) =
    Bool(GDAL.loadstyletable(styletable.ptr, filename))

"""
Get a style string by name.

### Parameters
* `styletable`: handle to the style table.
* `name`: the name of the style string to find.

### Returns
the style string matching the name or NULL if not found or error.
"""
find(styletable::StyleTable, name::AbstractString) =
    GDAL.find(styletable.ptr, name)

"""
Reset the next style pointer to 0.

### Parameters
* `styletable`: handle to the style table.
"""
resetreading!(styletable::StyleTable) =
    GDAL.resetstylestringreading(styletable.ptr)

"""
Get the next style string from the table.

### Parameters
* `styletable`: handle to the style table.

### Returns
the next style string or NULL on error.
"""
nextstyle(styletable::StyleTable) = GDAL.getnextstyle(styletable.ptr)

"""
Get the style name of the last style string fetched with OGR_STBL_GetNextStyle.

### Parameters
* `styletable`: handle to the style table.

### Returns
the Name of the last style string or NULL on error.
"""
laststyle(styletable::StyleTable) = GDAL.getlaststylename(styletable.ptr)


text:
module CommunalHelperDreamJellyfish

using ..Ahorn, Maple

@mapdef Entity "CommunalHelper/DreamJellyfish" DreamJellyfish(
    x::Integer,
    y::Integer,
    bubble::Bool=false,
    tutorial::Bool=false,
)

const placements = Ahorn.PlacementDict(
    "Dream Jellyfish (Communal Helper)" => Ahorn.EntityPlacement(
        DreamJellyfish,
    ),
    "Dream Jellyfish (Floating) (Communal Helper)" => Ahorn.EntityPlacement(
        DreamJellyfish,
        "point",
        Dict{String, Any}(
            "bubble" => true,
        ),
    )
)

const sprite = "objects/CommunalHelper/dreamJellyfish/jello"

function Ahorn.selection(entity::DreamJellyfish)
    x, y = Ahorn.position(entity)

    return Ahorn.Rectangle(x - 14, y - 15, 30, 19)
end

function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::DreamJellyfish, room::Maple.Room)
    Ahorn.drawSprite(ctx, sprite, 0, 0)

    if get(entity, "bubble", false)
        curve = Ahorn.SimpleCurve((-7, -1), (7, -1), (0, -6))
        Ahorn.drawSimpleCurve(ctx, curve, (1.0, 1.0, 1.0, 1.0), thickness=1)
    end
end

end

ext: jl


text:
function Newton_Stop(
  nlp::AbstractNLPModel,
  x₀::AbstractVector;
  τ₀::Float64 = 0.0005,
  stp::NLPStopping = NLPStopping(nlp, NLPAtX(x₀)),
  optimality_check = unconstrained_check,
  atol = 1e-6,
  max_iter = 200,
)
  x = copy(x₀)
  f, g = obj(nlp, x), grad(nlp, x)

  OK = update_and_start!(stp, x = x, fx = f, gx = g)

  while !OK
    Hx = hess(nlp, x)
    H = Matrix(Symmetric(Hx, :L))
    Δ, O = eigen(H)

    # Boost negative values of Δ to 1e-8
    D = Δ .+ max.((1e-8 .- Δ), 0.0)

    d = -O * diagm(1.0 ./ D) * O' * g

    # Armijo bactracking
    t, f, g = Armijo(nlp, f, g, x, d, τ₀)

    x += t * d
    OK = update_and_stop!(stp, x = x, gx = g, Hx = H)
  end

  if !stp.meta.optimal
    @warn "Optimality not reached"
    @info status(stp, list = true)
  end

  return x, f, norm(g), stp
end


text:
 __precompile__()

module Abstractions

export AnyInteger,
       AbsVector,   AbsMatrix,   AbsArray,
       NumVector,   NumMatrix,   NumArray,
       RealVector,  RealMatrix,  RealArray,
       FloatVector, FloatMatrix, FloatArray,
       IntVector,   IntMatrix,   IntArray,
       TimeVector,  TimeMatrix,  TimeArray

export AbstractPlace,    AbstractTime,
       AbstractSequence, AbstractMultiSequence,
       AbstractTimeSeries, AbstractTimeMultiSeries
using Base.Dates
import Base.Dates:AbstractTime
import Base.TimeSeries: AbstractTimeSeries

# from StatsBase.jl/common.jl
# common utilities

## convenient type alias
#
#  These types signficantly reduces the need of using
#  type parameters in functions (which are often just
#  for the purpose of restricting the arrays to real)
#
# These could be removed when the Base supports
# covariant type notation, i.e. AbstractVector{<:Real}
#

const RealArray{T<:Real,N} = AbstractArray{T,N}
const RealVector{T<:Real} = AbstractArray{T,1}
const RealMatrix{T<:Real} = AbstractArray{T,2}

const IntegerArray{T<:Integer,N} = AbstractArray{T,N}
const IntegerVector{T<:Integer} = AbstractArray{T,1}
const IntegerMatrix{T<:Integer} = AbstractArray{T,2}

const RealFP = Union{Float32, Float64}

## conversion from real to fp types

fptype{T<:Union{Float32,Bool,Int8,UInt8,Int16,UInt16}}(::Type{T}) = Float32
fptype{T<:Union{Float64,Int32,UInt32,Int64,UInt64,Int128,UInt128}}(::Type{T}) = Float64
fptype(::Type{Complex64}) = Complex64
fptype(::Type{Complex128}) = Complex128

# end from StatsBase

import Base.Dates: AbstractTime, TimeType
import TimeSeries: AbstractTimeSeries

const AnyInteger  = Union{Signed, Unsigned}

const AbsVector   = AbstractVector{T}  where T
const AbsMatrix   = AbstractMatrix{T}  where T
const AbsArray    = AbstractArray{T,N} where T where N

const NumVector   = AbstractVector{T}  where T<:Number
const NumMatrix   = AbstractMatrix{T}  where T<:Number
const NumArray    = AbstractArray{T,N} where T<:Number where N

const RealVector  = AbstractVector{T}  where T<:Real
const RealMatrix  = AbstractMatrix{T}  where T<:Real
const RealArray   = AbstractArray{T,N} where T<:Real where N

const FloatVector = AbstractVector{T}  where T<:AbstractFloat
const FloatMatrix = AbstractMatrix{T}  where T<:AbstractFloat
const FloatArray  = AbstractArray{T,N} where T<:AbstractFloat where N

const IntVector   = AbstractVector{T}  where T<:AnyInteger
const IntMatrix   = AbstractMatrix{T}  where T<:AnyInteger
const IntArray    = AbstractArray{T,N} where T<:AnyInteger where N

const TimeVector  = AbstractVector{T}  where T<:AbstractTime
const TimeMatrix  = AbstractMatrix{T}  where T<:AbstractTime
const TimeArray   = AbstractArray{T,N} where T<:AbstractTime where N

abstract type AbstractPlace    end

abstract type AbstractSequence end
abstract type AbstractMultiSequence <: AbstractSequence end

abstract type AbstractSeries   end
abstract type AbstractMultiSeries <: AbstractSeries   end

abstract type AbstractMultiTimeSeries <: AbstractTimeSeries end

abstract type AbstractSequence end

abstract type AbstractOrder end
abstract type PartialOrder <: AbstractOrder end
abstract type LatticeOrder <: PartialOrder  end
abstract type TotalOrder   <: LatticeOrder  end

export Abstraction,
    AnyVec, AnyMat, AnyArray,
    AnyTimeVec, AnyTimeMat, AnyTimeArray,
    AnyNumericVec, AnyNumericArray,
    AnyRealVec, AnyRealArray,
    AnyFloatVec, AnyFloatArray,
    SysFloat, SysInt, SysUInt

#       AbstractTime, AbstractSpace, AbstractInterval,
#       AbstractMeasure, AbstractValue, AbstractOperator

abstract type Abstraction end

abstract type AbstractType <: Abstraction end
abstract type AbstractKind <: Abstraction end

abstract type AbstractQuality  <: Abstraction end
abstract type AbstractQuantity <: Abstraction end

abstract type AbstractNeighborhood <: Abstraction end
abstract type AbstractRelationship <: Abstraction end

abstract type AbstractConnection <: Abstraction end
abstract type AbstractDirection  <: Abstraction end

const AnyTimeVec = AbstractVector{T} where T<:TimeType
const AnyNumericVec = AbstractVector{T} where T<:Number

const AnyTimeSeq = Union{AnyTimeVec, StepRange{P,Q}} where P<:TimeType where Q

end # module


text:
function ueberpruefeZeile(Zeile)
if round.(Zeile;digits=2) == [0 0 0.13 0 0.13 1.29 0.13 1.29 0]
    println("Eure Bestimmung der Zeile ist korrekt!")
else
    println("Überprüft noch einmal eure Rechnung und die Bedeutung der einzelnen Matrixeinträge.")
end
end


text:
# MIT License
#
# Copyright (c) 2018 Martin Biel
#
# 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.

"""
    SerialExecution

Functor object for using serial execution in a quasi-gradient algorithm. Create by supplying a [`Serial`](@ref) object through `execution` in `QuasiGradient.Optimizer` or by setting the [`Execution`](@ref) attribute.

"""
struct SerialExecution{T <: AbstractFloat, S <: AbstractSubProblem{T}} <: AbstractQuasiGradientExecution
    subproblems::Vector{S}
    decisions::DecisionMap

    function SerialExecution(structure::StageDecompositionStructure{2, 1, <:Tuple{ScenarioProblems}}, x::AbstractVector, subproblems::Unaltered, ::Type{T}) where T <: AbstractFloat
        execution = new{T,SubProblem{T}}(Vector{SubProblem{T}}(), structure.decisions[2])
        initialize_subproblems!(execution, scenarioproblems(structure, 2))
        return execution
    end

    function SerialExecution(structure::StageDecompositionStructure{2, 1, <:Tuple{ScenarioProblems}}, x::AbstractVector, subproblems::Smoothed, ::Type{T}) where T <: AbstractFloat
        execution = new{T,SmoothSubProblem{T}}(Vector{SmoothSubProblem{T}}(), structure.decisions[2])
        initialize_subproblems!(execution, scenarioproblems(structure, 2), x; type2dict(subproblems.parameters)...)
        return execution
    end
end

function initialize_subproblems!(execution::SerialExecution{T,SubProblem{T}}, scenarioproblems::ScenarioProblems) where T <: AbstractFloat
    for i in 1:num_subproblems(scenarioproblems)
        push!(execution.subproblems, SubProblem(
            subproblem(scenarioproblems, i),
            i,
            T(probability(scenario(scenarioproblems, i)))))
    end
    return nothing
end

function initialize_subproblems!(execution::SerialExecution{T,SmoothSubProblem{T}},
                                 scenarioproblems::ScenarioProblems,
                                 x::AbstractVector;
                                 kw...) where T <: AbstractFloat
    for vi in all_known_decisions(execution.decisions)
        # Unfix first-stage decisions
        execution.decisions[vi].state = NotTaken
    end
    # Load smooth subproblems (through Moreau envelope)
    for i in 1:num_subproblems(scenarioproblems)
        push!(execution.subproblems, SmoothSubProblem(
            subproblem(scenarioproblems, i),
            i,
            T(probability(scenario(scenarioproblems, i))),
            x; kw...))
    end
    return nothing
end

function resolve_subproblems!(quasigradient::AbstractQuasiGradient, execution::SerialExecution{T,SubProblem{T}}) where T <: AbstractFloat
    # Update subproblems
    update_known_decisions!(execution.decisions, quasigradient.x)
    # Initialize subgradient
    quasigradient.gradient .= quasigradient.c
    Q = zero(T)
    # Update and solve subproblems
    for subproblem in execution.subproblems
        update_subproblem!(subproblem)
        subgradient = subproblem(quasigradient.x)
        if isinf(subgradient.Q)
            return subgradient.Q
        end
        quasigradient.gradient .-= subgradient.δQ
        Q += subgradient.Q
    end
    # Return current objective value and subgradient
    return current_objective_value(quasigradient, Q)
end

function resolve_subproblems!(quasigradient::AbstractQuasiGradient, execution::SerialExecution{T,SmoothSubProblem{T}}) where T <: AbstractFloat
    # Update subproblems
    update_known_decisions!(execution.decisions, quasigradient.x)
    # Initialize subgradient
    quasigradient.gradient .= quasigradient.c
    Q = zero(T)
    # Update and solve subproblems
    for subproblem in execution.subproblems
        update_subproblem!(subproblem)
        gradient = subproblem(quasigradient.x)
        if isinf(gradient.Q)
            return gradient.Q
        end
        quasigradient.gradient .+= gradient.δQ
        Q += gradient.Q
    end
    # Return current objective value and subgradient
    return current_objective_value(quasigradient, Q)
end

function restore_subproblems!(::AbstractQuasiGradient, execution::SerialExecution{T,SubProblem{T}}) where T <: AbstractFloat
    return nothing
end

function restore_subproblems!(::AbstractQuasiGradient, execution::SerialExecution{T,SmoothSubProblem{T}}) where T <: AbstractFloat
    for subproblem in execution.subproblems
        restore_subproblem!(subproblem)
    end
    for vi in all_known_decisions(execution.decisions)
        # Remove common projection targets
        remove_decision!(execution.decisions, vi)
    end
    for vi in all_decisions(execution.decisions)
        # Re-fix first-stage decisions
        execution.decisions[vi].state = Known
    end
    return nothing
end

# API
# ------------------------------------------------------------
function (execution::Serial)(structure::StageDecompositionStructure{2, 1, <:Tuple{ScenarioProblems}},
                             x::AbstractVector,
                             subproblems::AbstractSubProblemState,
                             ::Type{T}) where T <: AbstractFloat
    return SerialExecution(structure, x, subproblems, T)
end

function str(::Serial)
    return ""
end


text:
# ------------------------------------------------------------------------------------------
# # Julia for Data Science - Algorithms
#
# ### Data processing
# In what's next, we will see how to use some of the standard algorithms for data analysis
# implemented in Julia. In particular, we'll look at
#
# 1. Kmeans clustering
# 1. Nearest neighbors with a KDTree
# 1. PCA for dimensionality reduction
# ------------------------------------------------------------------------------------------

using DataFrames, Statistics

# ------------------------------------------------------------------------------------------
# We'll be using the same data for all three of these examples -- the Sacramento real estate
# transactions file that we download next. This is a list of 985 real estate transactions in
# the Sacramento area reported over a five-day period.
# ------------------------------------------------------------------------------------------

download("http://samplecsvs.s3.amazonaws.com/Sacramentorealestatetransactions.csv","houses.csv")
houses = readtable("houses.csv")

# ------------------------------------------------------------------------------------------
# Let's use `Plots` to plot with the `pyplot` backend and start familiarizing ourselves with
# this data set!
# ------------------------------------------------------------------------------------------

using Plots
pyplot()
plot(size=(500,500),leg=false)

# ------------------------------------------------------------------------------------------
# Now let's create a scatter plot to show the price of a house vs. its square footage,
# ------------------------------------------------------------------------------------------

x = houses[!, :sq__ft]
# x = houses[7] # equivalent, useful if file has no header
y = houses[!, :price]
# y = houses[10] # equivalent
scatter(x,y,markersize=3)

# ------------------------------------------------------------------------------------------
# *Houses with 0 square feet that cost money?*
#
# The square footage seems to not have been recorded in these cases.
#
# Filtering these houses out is easy to do!
# ------------------------------------------------------------------------------------------

filter_houses = houses[houses[!, :sq__ft] .> 0, :]  # dot broadcasting
x = filter_houses[!, :sq__ft]
y = filter_houses[!, :price]
scatter(x,y)

# ------------------------------------------------------------------------------------------
# This makes sense! The higher the square footage, the higher the price.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# We can filter a `DataFrame` by feature value too, using the `by` function.
# ------------------------------------------------------------------------------------------

by(filter_houses,:type,size)

by(filter_houses,:type,filter_houses->mean(filter_houses[!, :price]))

# ------------------------------------------------------------------------------------------
# ### Example 1: Kmeans Clustering
#
# Now let's do some kmeans clustering on this data.
#
# First, we can load the `Clustering` package to do this.
# ------------------------------------------------------------------------------------------

#Pkg.add("Clustering")
using Clustering

# ------------------------------------------------------------------------------------------
# Let's store the features `:latitude` and `:longitude` in an array `X` that we will pass to
# `kmeans`.
# ------------------------------------------------------------------------------------------

X = filter_houses[!, [:latitude,:longitude]]
X = Array{Float64}(X)

# ------------------------------------------------------------------------------------------
# Each feature is stored as a row of `X`, but we can transpose to make these features
# columns of `X`.
# ------------------------------------------------------------------------------------------

X = transpose(X)

# ------------------------------------------------------------------------------------------
# Now let's plot longitudes vs. latitudes for this data!
# ------------------------------------------------------------------------------------------

location_figure = scatter(X[1, :], X[2, :])
# Alternatively....
# location_figure = scatter(filter_houses[!, :latitude], filter_houses[!, :longitude])
xlabel!("Latitude")
ylabel!("Longitude")
title!("Houses plotted by location")
display(location_figure)

# ------------------------------------------------------------------------------------------
# We want to identify clusters in the data above. As a first pass at guessing how many
# clusters we might need, let's use the number of zip codes in our data.
#
# (Try changing this to see how it impacts results!)
# ------------------------------------------------------------------------------------------

k = length(unique(filter_houses[!, :zip]))

# ------------------------------------------------------------------------------------------
# We can use the `kmeans` function to do kmeans clustering!
# ------------------------------------------------------------------------------------------

C = kmeans(X,k) # try changing k

# ------------------------------------------------------------------------------------------
# Now let's create a new data frame, `df`, with all the same data as `filter_houses` that
# also includes a column for the cluster to which each house has been assigned.
# ------------------------------------------------------------------------------------------

df = DataFrame(cluster = C.assignments,city = filter_houses[!,:city],
    latitude = filter_houses[!,:latitude],longitude = filter_houses[!,:longitude],zip = filter_houses[!,:zip])

# ------------------------------------------------------------------------------------------
# Let's plot each cluster as a different color.
# ------------------------------------------------------------------------------------------

clusters_figure = plot()
for i = 1:k
    # filter df to grab all houses in the ith cluster
    clustered_houses = df[df[!,:cluster].== i,:]
    # grab latitudes and longitudes of all houses in the ith cluster
    xvals = clustered_houses[!,:latitude]
    yvals = clustered_houses[!,:longitude]
    # plot latitudes and longitudes of all houses in the ith cluster
    scatter!(clusters_figure,xvals,yvals,markersize=4)
end
xlabel!("Latitude")
ylabel!("Longitude")
title!("Houses color-coded by cluster")
display(clusters_figure)

# ------------------------------------------------------------------------------------------
# And now let's try coloring them by zip code.
# ------------------------------------------------------------------------------------------

unique_zips = unique(filter_houses[!,:zip])
zips_figure = plot()
for uzip in unique_zips
    # filter houses by zipcode
    subs = filter_houses[filter_houses[!,:zip].==uzip,:]
    # grab the latitudes and longitudes of all houses in a given zipcode/subdivision
    x = subs[!, :latitude]
    y = subs[!, :longitude]
    # plot the houses in this zipcode by latitude and longitude!
    scatter!(zips_figure,x,y)
end
xlabel!("Latitude")
ylabel!("Longitude")
title!("Houses color-coded by zip code")
display(zips_figure)

# ------------------------------------------------------------------------------------------
# Let's see the two plots side by side.
# ------------------------------------------------------------------------------------------

plot(clusters_figure,zips_figure,layout=(2, 1))

# ------------------------------------------------------------------------------------------
# Not exactly! but almost... Now we know that ZIP codes are not randomly assigned!
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ### Example 2: Nearest Neighbor with a KDTree
#
# For this example, let's start by loading the `NearestNeighbors` package.
# ------------------------------------------------------------------------------------------

using NearestNeighbors

# ------------------------------------------------------------------------------------------
# With this package, we'll look for the `knearest` neighbors of one of the houses, `point`.
# ------------------------------------------------------------------------------------------

knearest = 10
id = 70 # try changing this
point = X[:,id]

# ------------------------------------------------------------------------------------------
# Now we can build a `KDTree` and use `knn` to look for `point`'s nearest neighbors!
# ------------------------------------------------------------------------------------------

kdtree = KDTree(X)
idxs, dists = knn(kdtree, point, knearest, true)

# ------------------------------------------------------------------------------------------
# We'll first generate a plot with all of the houses in the same color,
# ------------------------------------------------------------------------------------------

x = filter_houses[!,:latitude];
y = filter_houses[!,:longitude];
scatter(x,y);

# ------------------------------------------------------------------------------------------
# and then overlay the data corresponding to the nearest neighbors of `point` in a different
# color.
# ------------------------------------------------------------------------------------------

x = filter_houses[idxs,:latitude];
y = filter_houses[idxs,:longitude];
scatter!(x,y)

# ------------------------------------------------------------------------------------------
# There are those nearest neighbors in red!
#
# We can see the cities of the neighboring houses by using the indices, `idxs`, and the
# feature, `:city`, to index into the `DataFrame` `filter_houses`.
# ------------------------------------------------------------------------------------------

cities = filter_houses[idxs,:city]

# ------------------------------------------------------------------------------------------
# ### Example 3: PCA for dimensionality reduction
#
# Let us try to reduce the dimensions of the price/area data from the houses dataset.
#
# We can start by grabbing the square footage and prices of the houses and storing them in
# an `Array`.
# ------------------------------------------------------------------------------------------

F = filter_houses[[:sq__ft,:price]]
F = convert(Array{Float64,2},F)'

# ------------------------------------------------------------------------------------------
# Recall how the data looks when we plot housing prices against square footage.
# ------------------------------------------------------------------------------------------

scatter(F[1,:],F[2,:])
xlabel!("Square footage")
ylabel!("Housing prices")

# ------------------------------------------------------------------------------------------
# We can use the `MultivariateStats` package to run PCA
# ------------------------------------------------------------------------------------------

# Pkg.add("MultivariateStats")
using MultivariateStats

# ------------------------------------------------------------------------------------------
# Use `fit` to fit the model
# ------------------------------------------------------------------------------------------

M = fit(PCA, F)

# ------------------------------------------------------------------------------------------
# Note that you can choose the maximum dimension of the new space by setting `maxoutdim`,
# and you can change the method to, for example, `:svd` with the following syntax.
#
# ```julia
# fit(PCA, F; maxoutdim = 1,method=:svd)
# ```
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# It seems like we only get one dimension with PCA! Let's use `transform` to map all of our
# 2D data in `F` to `1D` data with our model, `M`.
# ------------------------------------------------------------------------------------------

y = transform(M, F)

# ------------------------------------------------------------------------------------------
# Let's use `reconstruct` to put our now 1D data, `y`, in a form that we can easily overlay
# (`Xr`) with our 2D data in `F` along the principle direction/component.
# ------------------------------------------------------------------------------------------

Xr = reconstruct(M, y)

# ------------------------------------------------------------------------------------------
# And now we create that overlay, where we can see points along the principle component in
# red.
#
# (Each blue point maps uniquely to some red point!)
# ------------------------------------------------------------------------------------------

scatter(F[1,:],F[2,:])
scatter!(Xr[1,:],Xr[2,:])
xlabel!("Square footage")
ylabel!("Housing prices")

max_stars_repo_name: grenkoca/JuliaTutorials

max_issues_repo_name: grenkoca/JuliaTutorials

max_forks_repo_name: grenkoca/JuliaTutorials


text:
# ---------------------------------------------------------------------------- #
#
#   setup.jl
#
#   Completely useless pretty script
#
#   λυτέος
#   Fall 2017
#
#   Max Opgenoord
#
# ---------------------------------------------------------------------------- #

"""
    setup( )

Completely useless pretty script.
"""
function setup( )

  println("   ")
  println("   ")
  println("       ▼▲ ")
  println("        ▼▲ ")
  println("        ▲▼▲    ■   ◣  ◢ ■ ■ ■ ◤  ◢ ■ ◣   ▲▼▲▼▲    ▲▼▲▼▲ ")
  println("       ▲▼ ▼▲   ■   ■      ■      ◥ ◣    ▲▼   ▼▲  ▲▼ ")
  println("      ▲▼   ▼▲  ■   ■      ■      ◢ ◤    ▼▲   ▲▼  ▼▲ ")
  println("     ▲▼     ▼▲ ◥ ■ ◤      ◥ ◤    ◥ ■ ◤   ▼▲▼▲▼    ▼▲▼▲ ")
  println("                                                    ▲▼ ")
  println("  --------------------------------------------------------")
  println("   ")
  println("                         FEM Library ")
  println("   ")

end


text:
@testset "PolygonalAmoeba" begin
    @polyvar x y
    f = x^2*y + y^2 + 3x^2*y^3 + y^4 + x^4*y^4

    A = amoeba(f, accuracy=0.01)
    @test 0 < accuracy(A) < 0.01
end


text:
using Requests

function make_badge_simple(subject::String, status::String, color::String, outfile::String)
    subject = replace(subject, " ", "%20")
    url = "https://img.shields.io/badge/$subject-$status-$color.svg"
    res = get(url)
    open(outfile, "w") do fid
        write(fid, res.data)
    end
end

"""
Make badge based on two numbers
"""
function make_badge(subject::String, a::Int, b::Int, outfile::String)
    limits = Float64[100, 80, 60, 40, 20, 0]
    colors = ["brightgreen", "green", "yellowgreen", "yellow", "orange", "red"]
    picked_color = "lightgray"
    perce = a/b*100
    println("perce: $perce")
    for (i, limit) in enumerate(limits)
        if perce >= limit
            picked_color = colors[i]
            break
        end
    end
    println("color: $picked_color")
    message = "$a/$b"
    make_badge_simple(subject, message, picked_color, outfile)
end

"""
Make badge based on one number
"""
function make_badge(subject::String, a::Int, outfile::ASCIIString)
    limits = Float64[0, 20, 40, 60, 80, 100]
    colors = ["brightgreen", "green", "yellowgreen", "yellow", "orange", "red"]
    picked_color = "red"
    for (i, limit) in enumerate(limits)
        if a <= limit
            picked_color = colors[i]
            break
        end
    end
    println("color: $picked_color")
    a = 100 - a
    if a < 0
        a = 0
    end
    message = "$a"
    make_badge_simple(subject, message, picked_color, outfile)
end


text:
module Reddit

# Inspired by Twitter.jl

using Requests, DataFrames, Codecs, JSON

#import DataFrames: DataFrame
import Requests: get, post
import Requests: Response
import Base: unique
export token!, cred, new, hot, gilded, controversial,
	   top, rising, revoke!, next!, get_all!, Subreddit,
	   @sorting

include("./Types.jl")
include("./Utilities.jl")

end # module Reddit


text:
using Test, CacheServers

struct FakeServer{K, ELT} <: CacheServers.AbstractCacheServer{K, ELT}
end

struct FakeFragment
end

s = FakeServer{Int, FakeFragment}()

# TODO: use custom Exception: InterfaceError
@testset "check apis" begin
    @test_throws MethodError alloc!(s, 1, FakeFragment())
    @test_throws MethodError iscacheable(s, 1)
    @test_throws MethodError iscached(s, 1)
    @test_throws MethodError push!(s, 1, 1)
    @test_throws MethodError update!(FakeFragment(), 1)
    @test_throws MethodError pull(s, 1)
    @test_throws MethodError getindex(s, 1)
    @test_throws MethodError delete!(s, 1)
    @test_throws MethodError clear!(s, 1)
end


text:
# Tests using TaylorModel1 and RTaylorModel1

using TaylorModels
using LinearAlgebra: norm
using Test

const _num_tests = 1000

setformat(:full)

function check_containment(ftest, xx::TaylorModelN{N,T,S}, tma::TaylorModelN{N,T,S}) where {N,T,S}
    xfp = diam.(tma.I) .* (rand(N) .- 0.5) .+ mid(tma.x0)
    xbf = [big(xfp[i]) for i=1:N]
    ib = IntervalBox([@interval(xfp[i]) for i=1:N]...)
    range = evaluate(tma, ib-tma.x0)
    bb = all(ftest(xx(xbf .- mid(tma.x0))) ⊆ range)
    bb || @show(ftest, ib, xbf, ftest(xbf...), range)
    return bb
end

const _order = 2
const _order_max = 2*(_order+1)
set_variables(Interval{Float64}, [:x, :y], order=_order_max)

@testset "Tests for TaylorModelN " begin
    b0 = Interval(0.0) × Interval(0.0)
    ib0 = Interval(-0.5, 0.5) × Interval(-0.5, 0.5)
    b1 =  Interval(0.0) × Interval(1.0)
    ib1 = Interval(-0.5, 0.5) × Interval(0.5, 1.5)

    zi = 0..0
    xT = TaylorN(Interval{Float64}, 1, order=_order)
    yT = TaylorN(Interval{Float64}, 2, order=_order)

    @testset "TaylorModelN constructors" begin
        xm = TaylorModelN{2, Interval{Float64}, Float64}(xT, zi, b0, ib0)
        ym = TaylorModelN{2, Interval{Float64}, Float64}(yT, zi, b0, ib0)
        @test xm == TaylorModelN(xT, zi, b0, ib0)
        @test ym == TaylorModelN(yT, zi, b0, ib0)
        @test xm == TaylorModelN(1, _order, b0, ib0)
        @test ym == TaylorModelN(2, _order, b0, ib0)
        @test TaylorModelN( b1[1], 2, b0, ib0) ==
                TaylorModelN(TaylorN(b1[1], _order), zi, b0, ib0)

        @test isa(xm, AbstractSeries)
        @test TaylorModelN{2, Interval{Float64},Float64} <: AbstractSeries{Interval{Float64}}

        # Test errors in construction
        @test_throws AssertionError TaylorModelN(xT, zi, IntervalBox(1..1), IntervalBox(1..1))
        @test_throws AssertionError TaylorModelN(xT, zi, b0, ib1)
        @test_throws AssertionError TaylorModelN(xT, 1..1, b0, ib0)
        @test_throws BoundsError TaylorModelN(5, _order, b0, ib0) # wrong variable number

        # Tests for get_order and remainder
        @test get_order() == 6
        @test get_order(xm) == 2
        @test remainder(ym) == zi
    end

    @testset "Arithmetic operations" begin
        Δ = interval(-0.25, 0.25)
        xm = TaylorModelN(xT, zi, b1, ib1)
        ym = TaylorModelN(yT, zi, b1, ib1)
        a = TaylorModelN( b1[1]+xT, Δ, b1, ib1)
        @test zero(a) == TaylorModelN(zero(a.pol), 0..0, b1, ib1)
        @test one(a) == TaylorModelN(one(a.pol), 0..0, b1, ib1)
        @test a + a == TaylorModelN(2*(b1[1]+ xT), 2*Δ, b1, ib1)
        @test -a == TaylorModelN( -(b1[1]+xT), -Δ, b1, ib1)
        @test a - a == TaylorModelN(zero(a.pol), 2*Δ, b1, ib1)
        @test b1[2] + ym == TaylorModelN(b1[2] + yT, zi, b1, ib1)
        @test a - b1[1] == TaylorModelN(zero(b1[1])+xT, Δ, b1, ib1)

        @test xm * ym == TaylorModelN( xT * yT, zi, b1, ib1)
        b = a * ym
        @test b == TaylorModelN( xT * yT, Δ*ym.pol(ib1-b1), b1, ib1)
        b = ym * TaylorModelN(xT^2, zi, b1, ib1)
        @test b == TaylorModelN( zero(xT), (ib1[1]-b1[1])^2*(ib1[2]-b1[2]), b1, ib1 )
        b = b1[2] * a
        @test b == TaylorModelN( b1[2]*a.pol, Δ*b1[2], b1, ib1 )
        @test b / b1[2] == a
        @test_throws AssertionError TaylorModelN(TaylorN(1, order=_order_max), zi, b1, ib1) *
            TaylorModelN(TaylorN(2, order=_order_max), zi, b1, ib1)

        remt = remainder(1/(1-TaylorModel1(_order, b1[1], ib1[1])))
        @test remainder(1 / (1-xm)) == remt
        @test remainder(ym / (1-xm)) == Interval(-0.25, 0.25)

        @test remainder(xm^2) == remainder(ym^2)
        @test (xm.I[1]-xm.x0[1])^3 == remainder(xm^3)
        @test (ym.I[2]-ym.x0[2])^4 ⊆ remainder(ym^4)
    end

    @testset "RPAs, functions and remainders" begin
        xm = TaylorModelN(1, _order, b1, ib1)
        ym = TaylorModelN(2, _order, b1, ib1)

        @test rpa(x->5+zero(x), xm) == 5+zero(xm)
        @test rpa(x->5+one(x), ym) == 5+one(ym)
        @test rpa(x->5*x, ym) == 5*ym
        remT = remainder(5*TaylorModel1(2, b1[1], ib1[1])^4)
        @test rpa(x->5*x^4, xm) == TaylorModelN(zero(xT), remT, b1, ib1)
        remT = 5 * (ib1[1]-b1[1])^2 * (2*(ib1[2]-b1[2])+(ib1[2]-b1[2])^2)
        @test rpa(x->5*x^2, xm*ym) == TaylorModelN( 5*xT^2, remT, b1, ib1)

        # Testing remainders of an RPA
        ftest = x -> exp(x)-1
        xx = xm + ym
        tma = rpa(ftest, xx)
        tmb = ftest(xx)
        @test tma == tmb
        # fT, Δ, ξ0 = rpafp(tma)
        # @test interval(ftest(ii.lo)-fT(ii.lo-ξ0),
        #                 ftest(ii.hi)-fT(ii.hi-ξ0)) ⊆ remainder(tma)
        for ind = 1:_num_tests
            @test check_containment(ftest, xx, tma)
        end

        ftest = x -> exp(x)
        xx = xm + 2*ym
        tma = rpa(ftest, xx)
        tmb = ftest(xx)
        @test tma == tmb
        # fT, Δ, ξ0 = rpafp(tma)
        # @test interval(ftest(ii.lo)-fT(ii.lo-ξ0),
        #                 ftest(ii.hi)-fT(ii.hi-ξ0)) ⊆ remainder(tma)
        for ind = 1:_num_tests
            @test check_containment(ftest, xx, tma)
        end

        ftest = x -> sin(x)
        xx = 1- xm^2 + ym
        tma = rpa(ftest, xx)
        tmb = ftest(xx)
        @test tma == tmb
        # fT, Δ, ξ0 = rpafp(tma)
        # @test interval(ftest(ii.lo)-fT(ii.lo-ξ0),
        #                 ftest(ii.hi)-fT(ii.hi-ξ0)) ⊆ remainder(tma)
        for ind = 1:_num_tests
            @test check_containment(ftest, xx, tma)
        end

        ftest = x -> sqrt(x)
        xx = xm^2 + ym
        tma = rpa(ftest, xx)
        tmb = ftest(xx)
        @test tma == tmb
        # fT, Δ, ξ0 = rpafp(tma)
        # @test interval(ftest(ii.lo)-fT(ii.lo-ξ0),
        #                 ftest(ii.hi)-fT(ii.hi-ξ0)) ⊆ remainder(tma)
        for ind = 1:_num_tests
            @test check_containment(ftest, xx, tma)
        end

        ftest = x -> inv(x)
        xx = 1 + xm + ym
        tma = rpa(ftest, xx)
        tmb = ftest(xx)
        @test tma == tmb
        # fT, Δ, ξ0 = rpafp(tma)
        # @test interval(ftest(ii.lo)-fT(ii.lo-ξ0),
        #                 ftest(ii.hi)-fT(ii.hi-ξ0)) ⊆ remainder(tma)
        for ind = 1:_num_tests
            @test check_containment(ftest, xx, tma)
        end
    end

    @testset "Composition of functions and their inverses" begin
        xm = TaylorModelN(xT, zi, b1, ib1)
        ym = TaylorModelN(yT, zi, b1, ib1)

        tma = exp(ym)
        tmb = log(tma)
        @test tmb == log(exp(ym))
        @test tmb.pol == ym.pol

        tma = sin(xm)
        tmb = asin(tma)
        @test tmb == asin(sin(xm))
        @test tmb.pol == xm.pol

        tma = asin(xm)
        tmb = sin(tma)
        @test tmb == sin(asin(xm))
        @test tmb.pol == xm.pol

        tma = acos(ym)
        tmb = cos(tma)
        @test tmb == cos(acos(ym))
        @test sup(norm(tmb.pol - ym.pol, Inf)) < 1.0e-16

        tma = tan(xm)
        tmb = atan(tma)
        @test tmb == atan(tan(xm))
        @test tmb.pol == xm.pol

        tma = atan(xm)
        tmb = tan(tma)
        @test tmb == tan(atan(xm))
        @test tmb.pol == xm.pol

        ####
        tma = log(1+ym)
        tmb = exp(tma)
        @test tmb == exp(log(1+ym))
        @test tmb.pol == (1+ym).pol

        tma = cos(1+ym)
        tmb = acos(tma)
        @test tmb == acos(cos(1+ym))
        @test sup(norm(tmb.pol - (1+ym).pol, Inf)) < 1.0e-15
    end

    @testset "Display" begin
        xm = TaylorModelN(1, _order, b1, ib1)
        ym = TaylorModelN(2, _order, b1, ib1)
        use_show_default(true)
        @test string(xm+ym) == "TaylorModelN{2,Interval{Float64},Float64}" *
            "(TaylorN{Interval{Float64}}" *
            "(HomogeneousPolynomial{Interval{Float64}}" *
            "[HomogeneousPolynomial{Interval{Float64}}" *
            "(Interval{Float64}[Interval(1.0, 1.0)], 0), " *
            "HomogeneousPolynomial{Interval{Float64}}" *
            "(Interval{Float64}[Interval(1.0, 1.0), Interval(1.0, 1.0)], 1), " *
            "HomogeneousPolynomial{Interval{Float64}}" *
            "(Interval{Float64}[Interval(0.0, 0.0), Interval(0.0, 0.0), " *
            "Interval(0.0, 0.0)], 2)], 2), Interval(0.0, 0.0), IntervalBox(Interval(0.0, 0.0), " *
            "Interval(1.0, 1.0)), IntervalBox(Interval(-0.5, 0.5), Interval(0.5, 1.5)))"
        use_show_default(false)
        @test string((xm+ym)^2) == " Interval(1.0, 1.0) + Interval(2.0, 2.0) x + " *
            "Interval(2.0, 2.0) y + Interval(1.0, 1.0) x² + Interval(2.0, 2.0) x y + " *
            "Interval(1.0, 1.0) y² + Interval(0.0, 0.0)"
    end

end


text:
using Weber
using Base.Test
include("find_timing.jl")

loop_events,_,rows = find_timing() do
  @addtrials let i = 0
    @addtrials while i < 3
      addtrial(moment(() -> i+=1),
               moment(() -> record(:a)),
               moment(() -> record(:b)),
               moment(() -> record(:c)))
    end
  end
end
loop_index = map(x -> (x[:trial],x[:offset]),rows)

@testset "Looping Moments" begin
  @test loop_events == [:a,:b,:c,:a,:b,:c,:a,:b,:c]
  @test loop_index == [(1,1),(1,1),(1,1),(2,1),(2,1),(2,1),(3,1),(3,1),(3,1)]
end


text:
## ------------------------------------------------------
mutable struct SimD3
    # Space
    V::Space

    # Chemostat
    X::Float64
    sg::Float64
    cg::Float64
    ϵ::Float64
    Δt::Float64
    D::Float64

    # Simulation
    it::Int
    niters::Int
    dXdt::Float64
    z_av::Float64
    ug_av::Float64
    uo_av::Float64
    cgD_X::Float64
    P::Vector{Float64}

    function SimD3(;
            # Space
            V::Space,
            # Chemostat
            X::Float64,
            cg::Float64,
            sg::Float64,
            ϵ::Float64,
            Δt::Float64,
            D::Float64,
            # Simulation
            it::Int = 1,
            niters::Int = 5000,
            dXdt::Float64 = 0.0,
            z_av::Float64 = 0.0,
            ug_av::Float64 = 0.0,
            uo_av::Float64 = 0.0,
            cgD_X::Float64 = 0.0,
            P::Vector{Float64} = ones(length(V))
        )

        # Space
        @assert length(P) == length(V)
        normalizeP!(P)

        new(V, X, sg, cg, ϵ, Δt, D, it, niters, dXdt, z_av, ug_av, uo_av, cgD_X, P)
    end
end

## ------------------------------------------------------
function Base.show(io::IO, S::SimD3)
    println(io, "SimD3")
    println(io, "niters: ", S.niters)
    println(io, S.V)
end
Base.show(S::SimD3) = show(stdout, S)

## ------------------------------------------------------
function run_simD3!(S::SimD3;
        dobreak::Function,
        tranformation::Function,
        feedback::Function = () -> nothing
    )

    V = S.V
    P = S.P
    Vz = Vi(V, :z)
    Vug = Vi(V, :ug)
    Vuo = Vi(V, :uo)

    normalizeP!(P)

    # globals
    cgD_X = 0.0
    dXidt = similar(P)
    loc = similar(P)
    glob = similar(P)
    drain = similar(P)
    vol = length(V)

    it0 = (S.it - 1)
    for _ in it0:S.niters
        S.it += 1

        # dXidt = loc + glob - drain
        loc .=  (1.0 - S.ϵ) .* Vz .* S.X .* P
        glob .=  S.ϵ * sum(Vz .* S.X .* P) / vol
        drain .=  S.D .* S.X .* P
        dXidt .= (loc .+ glob .- drain) .* S.Δt

        # update P
        P .= (P .+ (dXidt ./ S.X))
        normalizeP!(P)
        S.dXdt = sum(dXidt)

        # update X
        S.X += S.dXdt

        # P transformation
        _ug_av = sum(Vug .* P)
        _cgD_X = S.cg * S.D / S.X
        if (_ug_av > _cgD_X)
            tranformation()
            normalizeP!(P)
        end

        # update S
        S.cgD_X = _cgD_X

        S.z_av = sum(Vz .* P)
        S.ug_av = sum(Vug .* P)
        S.uo_av = sum(Vuo .* P)

        S.sg += (- S.ug_av * S.X + (S.cg - S.sg) * S.D) * S.Δt
        S.sg = max(S.sg, 0.0)

        # feedback
        feedback()

        # stst
        (dobreak() === true) && break

    end # for it

    return S

end


text:
using SafeTestsets

@safetestset "Element tests" begin include("element_test.jl") end

@safetestset "Mesh tests" begin include("mesh_test.jl") end

@safetestset "CartesianMesh tests" begin include("cartesianmesh_test.jl") end

@safetestset "Meshgen tests" begin include("meshgen_test.jl") end

max_stars_repo_name: WaveProp/WaveProp

max_issues_repo_name: WaveProp/WaveProp

max_forks_repo_name: WaveProp/WaveProp


text:
@testset "Vanilla VAE" begin
    Random.seed!(0)
    T = Float32
    xlen = 4
    zlen = 2
    batch = 20
    test_data = randn(T, 4, batch)/100 .+ hcat(ones(T,xlen,Int(batch/2)), -ones(T,xlen,Int(batch/2))) |> gpu

    enc = Chain(Dense(xlen, xlen, relu),
                Dense(xlen, xlen, relu),
                Dense(xlen, xlen, relu),
                SplitLayer(xlen, [zlen,zlen], [identity,softplus]))
    enc_dist = ConditionalMvNormal(enc)

    dec = Chain(Dense(zlen, xlen, relu),
                Dense(xlen, xlen, relu),
                Dense(xlen, xlen, relu),
                SplitLayer(xlen, [xlen,1], [identity,softplus]))
    dec_dist = ConditionalMvNormal(dec)

    model = VAE(zlen, enc_dist, dec_dist) |> gpu

    loss = - elbo(model, test_data)
    ps = Flux.params(model)
    @test length(ps) > 0
    @test isa(loss, Real)

    zs = rand(model.encoder, test_data)
    @test size(zs) == (zlen, batch)
    xs = mean(model.decoder, zs)
    @test size(xs) == (xlen, batch)

    # test training
    params_init = get_params(model)
    opt = ADAM()
    data = [(test_data,) for i in 1:2000]
    lossf(x) = - elbo(model, x, β=1e0)
    Flux.train!(lossf, Flux.params(model), data, opt)

    @test all(param_change(params_init, model)) # did the params change?
    nxs = mean(model.decoder, rand(model.encoder, test_data))
    @debug maximum(test_data - nxs)
    @test Flux.mse(test_data, nxs) < Flux.mse(test_data, xs)
    # this seems to be a stable way to test reconstruction
end

@testset "Wasserstein VAE" begin
    Random.seed!(0)

    T = Float32
    xlen = 4
    zlen = 2
    batch = 20
    test_data = hcat(ones(T,xlen,Int(batch/2)), -ones(T,xlen,Int(batch/2))) |> gpu

    enc = GenerativeModels.stack_layers([xlen, 10, 10, zlen], relu, Dense)
    enc_dist = ConditionalMvNormal(enc)

    dec = GenerativeModels.stack_layers([zlen, 10, 10, xlen], relu, Dense)
    dec_dist = ConditionalMvNormal(dec)

    model = VAE(zlen, enc_dist, dec_dist) |> gpu

    # test training
    params_init = get_params(model)
    opt = ADAM()
    k = IMQKernel(1f0)
    mmd(x) = GenerativeModels.mmd_mean(model, x, k)
    data = (test_data,)
    lossf(x) = Flux.mse(x, mean(model.decoder, mean(model.encoder,x))) + mmd(x)
    GenerativeModels.update_params!(model, data, lossf, opt)
    ps = Flux.params(model)
    @test all(param_change(params_init, model))

    # this works well but has quite a large variance
    data = [(test_data,) for _ in 1:10000]
    Flux.train!(lossf, Flux.params(model), data, opt)
    zs = mean(model.encoder, test_data)
    xs = mean(model.decoder, zs)
    @debug maximum(abs.(test_data - xs))
    @test all(abs.(test_data - xs) .< 0.2)

    msg = sprint(show, model)
    @test occursin("VAE", msg)
    Random.seed!()  # reset the seed
end


text:
export EpitopeMap

struct EpitopeMap
    name::String
    epitopes::Vector{String}
    starts::Vector{Int}
    stops::Vector{Int}
    alleles::Vector{HLAAllele}
    maplength::Int

    function EpitopeMap(name, epitopes, starts, stops, alleles, maplength)
        if !(length(epitopes) == length(starts) == length(stops))
            error("epitope sequence and start and stop positions" *
                  "must have the same length.")
        elseif !all(starts .<= stops)
            error("epitope start must be less than epitope stop.")
        elseif !(all((starts .> 0) .& (stops .> 0)))
            error("start and stop positions must be positive.")
        elseif !all(length.(epitopes) .== (stops .- starts .+ 1))
            error("epitope length must match the given start and stop positions.")
        elseif (length(stops) != 0) && !(maplength >= maximum(stops))
            error("maplength must be >= largest stop position.")
        end

        new(name , epitopes, starts, stops, alleles, maplength)
    end
end

function epitope_map(data::AbstractHLAData; rank_threshold::Real = 100)
    df = epitope_prediction(data, rank_threshold = rank_threshold)

    epitopes = df[!, :peptide]
    starts = df[!, :pos] .+ 1
    stops = starts .+ length.(df[!, :peptide]) .- 1
    alleles = parse_allele.(df[!, :best_allele])
    maplength = sequence_length(data)

    map = EpitopeMap(name(data), epitopes, starts, stops, alleles, maplength)

    return map
end

function Base.in(x::Pair{Replacement, HLAAllele}, y::EpitopeMap)
    for (i, (start, stop)) in enumerate(zip(y.starts, y.stops))
        if (start <= x[1].position <= stop) & (y.alleles[i] == x[2])
            return true
        end
    end

    return false
end

function match_epitope(epitope::String, data::AbstractHLAData;
        max_levenshtein_distance::Int = 0
)
    position = Int[]
    query = ApproximateSearchQuery(LongAminoAcidSeq(epitope))

    for record in records(data)
        sequence = BioSequences.sequence(record)
        sequence_string = replace(string(sequence), 'X' => '-')
        sequence = LongAminoAcidSeq(sequence_string)

        idx = approxsearchindex(sequence, query, max_levenshtein_distance)

        if idx != 0
            push!(position, idx)
        end
    end

    length(unique(position)) == 1 || return nothing

    return position[1]
end


text:
#From http://home.online.no/~pjacklam/notes/invnorm/

function CTΦInv(p::Float64)::Float64

    a1::Float64 = -39.6968302866538
    a2::Float64 = 220.946098424521
    a3::Float64 = -275.928510446969
    a4::Float64 = 138.357751867269
    a5::Float64 = -30.6647980661472
    a6::Float64 = 2.50662827745924

    b1::Float64 = -54.4760987982241
    b2::Float64 = 161.585836858041
    b3::Float64 = -155.698979859887
    b4::Float64 = 66.8013118877197
    b5::Float64 = -13.2806815528857

    c1::Float64 = -7.78489400243029E-03
    c2::Float64 = -0.322396458041136
    c3::Float64 = -2.40075827716184
    c4::Float64 = -2.54973253934373
    c5::Float64 = 4.37466414146497
    c6::Float64 = 2.93816398269878

    d1::Float64 = 7.78469570904146E-03
    d2::Float64 = 0.32246712907004
    d3::Float64 = 2.445134137143
    d4::Float64 = 3.75440866190742

    #Rational approximation for lower region.

    if p < 0.02425
        q::Float64 = √(-2.0 * log(p))
        ans::Float64 = (((((c1 * q + c2) * q + c3) * q + c4) * q + c5) * q + c6) / ((((d1 * q + d2) * q + d3) * q + d4) * q + 1.0)
    elseif p < 0.97575
        q = (p - 0.5) * (p - 0.5)
        ans = (((((a1 * q + a2) * q + a3) * q + a4) * q + a5) * q + a6) * (p - 0.5) / (((((b1 * q + b2) * q + b3) * q + b4) * q + b5) * q + 1)
    else
        q = √(-2 * log(1 - p))
        ans = -(((((c1 * q + c2) * q + c3) * q + c4) * q + c5) * q + c6) / ((((d1 * q + d2) * q + d3) * q + d4) * q + 1)
    end

    return ans
end

CTΦ(z::Real)::Real = 0.5+0.5*erf(z*sqrt2Inv)
CTϕ(z::Real)::Real = sqrt2PiInv*exp(-z^2.0/2.0)

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
using Pkg
Pkg.activate(".")
using RadMod2D

# model
include("./models2D.jl")
r1x = 0.8
r1y = 0.8
r2x = 1.2
r2y = 1.2
elemsize = 0.05
m = model_rect_in_rect(r1x, r1y, r2x, r2y, elemsize)

# view factors
vfmat = zeros(Float64, m.nelem, m.nelem)
@time existing_vf!(m, vfmat)
n = 15
dx, dy = create_tiles(m, n)
@time t_occ = check_tile_occupation(m, dx, dy, n)
@time blocking_vf_with_tiles!(m, vfmat, dx, dy, n, t_occ)
@time calculating_vf!(m, vfmat, normit = true)
vfmatp = vfmat_to_parts(m, vfmat, normit = true)
println("vf mat for parts:")
for i = 1:size(vfmatp,1)
    println(vfmatp[i,:])
end

# net radiation method
epsilon = zeros(m.nelem,1)
set_bc_part!(m, epsilon, 1, 0.3)
set_bc_part!(m, epsilon, 2, 0.6)
temp = zeros(m.nelem,1)
set_bc_part!(m, temp, 1:2, 300)
temp[79:106,1] .= 600 # bc for specific elements
@time Q, G = tempsolver(m, vfmat, temp, epsilon)
area = [m.elem[i].area for i = m.elem2par[1].first:m.elem2par[end].last]
q = Q[:] ./ area[:]

println("calculation finished")

# include("./test/run_example1.jl")


text:
using RollingFunctions
using Statistics
using StatsBase
using LinearAlgebra
using Test

datavec = collect(1.0f0:5.0f0)
weighting = normalize([1.0f0, 2.0f0, 4.0f0])

windowsize = 3

@test rollmean(datavec, windowsize) == Float32[2.0, 3.0, 4.0]

expected = [2.0f0, 3.0f0, 4.0f0]
obtained = rolling(mean, datavec, 3)
@test Float32(eps(Float64)) > abs(sum(expected .- obtained))

expected = [1.236568f0, 1.7457432f0, 2.254918f0]
obtained = rollmean(datavec, windowsize, weighting)
@test Float32(eps(Float32)) > abs(sum(expected .- obtained))

expected = [1.0f0, 1.1180339f0, 1.236568f0, 1.7457432f0, 2.254918f0]
obtained = running(mean, datavec, 3, weighting)
@test Float32(eps(Float32)) > abs(sum(expected .- obtained))

data1=[1,2,3,4,5,6,7];data2=[1,4,6,16,25,36,49];
weighting = normalize([1.0, 2.0, 1.0, 2.0]);

expected = [2.5, 6.0, 9.5, 10.0, 12.0]
obtained = rolling(cov, data1, data2, 3)
@test eps(Float32) > abs(sum(expected .- obtained))

expected = [0.9751641759537001, 0.9283031155706918, 0.9749534026403799, 0.9732162133211492]
obtained = rollcor(data1, data2, 4, weighting)
@test eps(Float32) > abs(sum(expected .- obtained))

expected = [0.0, 1.5, 2.5, 7.833333333333333, 15.0, 20.75, 26.5]
obtained = runcov(data1, data2, 5)
@test eps(Float32) > abs(sum(expected .- obtained))


text:
@testset "Virtual" begin
  include("SequenceControllerTest.jl")
end


text:
coord = TrackCoord(0, 1.0 ,2.0 ,3.0 ,4.0 ,5.0 ,6.0)
@test coord.trackid == 0
@test coord.s == 1.0
@test coord.t == 2.0
@test coord.z == 3.0
@test coord.h == 4.0
@test coord.p == 5.0
@test coord.r == 6.0

@test coord == TrackCoord(0, 1.0 ,2.0 ,3.0 ,4.0 ,5.0 ,6.0)

coordIM = TrackCoordIm(0, 1.0, 2.0)
@test coordIM.trackid == 0
@test coordIM.s == 1.0
@test coordIM.t == 2.0
@test coordIM.z == 0.0
@test coordIM.h == 0.0
@test coordIM.p == 0.0
@test coordIM.r == 0.0

coordIM = TrackCoordIm(0, 1.0 ,2.0 ,3.0 ,4.0 ,5.0 ,6.0)
@test coordIM.s == 1.0
@test coordIM.t == 2.0
@test coordIM.z == 3.0
@test coordIM.h == 4.0
@test coordIM.p == 5.0
@test coordIM.r == 6.0

@test coordIM == TrackCoordIm(0, 1.0, 2.0, 3.0, 4.0 ,5.0, 6.0)

show(coord); print("\n")
println(coord)

@test copy(coord) == TrackCoord(0, 1.0 ,2.0 ,3.0 ,4.0 ,5.0 ,6.0)
@test deepcopy(coord) == TrackCoord(0, 1.0 ,2.0 ,3.0 ,4.0 ,5.0 ,6.0)
@test copy!(TrackCoord(), coord) == TrackCoord(0, 1.0 ,2.0 ,3.0 ,4.0 ,5.0 ,6.0)
@test init!(copy(coord)) == TrackCoord(0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0)

coordA = TrackCoord(0, 1.0 ,2.0 ,3.0 ,4.0 ,5.0 ,6.0)
coordB = TrackCoord(1, 0.5 ,1.5 ,2.5 ,3.5 ,4.5 ,5.5)
@test plus!(coordA, coordB) == TrackCoord(0, 1.5 ,3.5 ,5.5 ,7.5 ,9.5 ,11.5)

coordim = TrackCoordIm()
coordim = TrackCoordIm(0, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5)

@test isa(convert(TrackCoordIm, coordA), TrackCoordIm)
@test isa(convert(TrackCoord, coordA), TrackCoord)
@test isa(convert(TrackCoordIm, coordim), TrackCoordIm)
@test isa(convert(TrackCoord, coordim), TrackCoord)

@test convert(TrackCoord, coordA) == coordA
# @test convert(TrackCoord, coordim) == coordim
# @test convert(TrackCoordIm, coordA) == coordA
@test convert(TrackCoordIm, coordim) == coordim

max_stars_repo_name: sisl/OpenDriveManagerLite

max_issues_repo_name: sisl/OpenDriveManagerLite

max_forks_repo_name: sisl/OpenDriveManagerLite


text:
using Base.Iterators: partition
using Flux
using Flux.Optimise: update!
using Flux: logitbinarycrossentropy, logitcrossentropy
using Images
using MLDatasets
using Statistics
using Parameters: @with_kw
using Random
using Printf
using CUDAapi
using Zygote
if has_cuda()		# Check if CUDA is available
    @info "CUDA is on"
    import CuArrays		# If CUDA is available, import CuArrays
    CuArrays.allowscalar(false)
end

@with_kw struct HyperParams
    batch_size::Int = 128
    latent_dim::Int = 100
    nclasses::Int = 10
    epochs::Int = 30
    verbose_freq::Int = 1500
    output_x::Int = 10
    output_y::Int = 10
    lr_dscr::Float64 = 0.0001
    lr_gen::Float64 = 0.0001
end

mutable struct discriminator
    discr
    aux
    model
end

function discriminator(args)
    model = Chain(Conv((3,3), 1=>128, pad=(1,1), stride = (2,2)),
                  x->leakyrelu.(x, 0.2f0),
                  Dropout(0.4),
                  Conv((3,3), 128=>128, pad= (1,1), stride = (2,2), leakyrelu),
                  x->leakyrelu.(x, 0.2f0),
                  x -> reshape(x, :, size(x, 4)),
                  Dropout(0.4)) |> gpu
    discr = Dense(6272, 1, sigmoid) |> gpu
    aux = Chain(Dense(6272, 10), softmax) |> gpu
    discriminator(discr,aux,model)
end

mutable struct generator
    m1
    m2
    model
end

function generator(args)
    m1 = Chain(Dense(args.nclasses, 49), x-> reshape(x, 7 , 7 , 1 , size(x,2))) |> gpu
    m2 = Chain(Dense(args.latent_dim, 6272), x->leakyrelu.(x, 0.2f0), x-> reshape(x, 7,7,128, size(x,2))) |> gpu
    model = Chain(ConvTranspose((4, 4), 129 => 128; stride = 2, pad = 1),
            BatchNorm(128, leakyrelu),
            Dropout(0.25),
            ConvTranspose((4, 4), 128 => 64; stride = 2, pad = 1),
            BatchNorm(64, leakyrelu),
            Conv((7, 7), 64 => 1, tanh; stride = 1, pad = 3)) |> gpu
    return generator(m1, m2, model)
end

function (m::generator)(x, y)
    t = cat(m.m1(x), m.m2(y), dims = 3)
    m.model(t)
end

function load_data(hparams)
    # Load MNIST dataset
    images, labels = MLDatasets.FashionMNIST.traindata(Float32)
    # Normalize to [-1, 1] and convert it to WHCN
    image_tensor = permutedims(reshape(@.(2f0 * images - 1f0), 28, 28, 1, :), (2, 1, 3, 4))
    y = float.(Flux.onehotbatch(labels, 0:hparams.nclasses-1))
    # Partition into batches
    data = [(image_tensor[:, :, :, r], y[:,r]) |> gpu for r in partition(1:60000, hparams.batch_size)]
    return data
end

# Loss functions
function discr_loss(real_output, fake_output)
    real_loss = mean(logitbinarycrossentropy.(real_output, 1f0))
    fake_loss = mean(logitbinarycrossentropy.(fake_output, 0f0))
    return real_loss + fake_loss
end

function aux_loss(output, labels)
   return logitcrossentropy(output, labels)
end

generator_loss(fake_output) = mean(logitbinarycrossentropy.(fake_output, 1f0))

function train_discr(discr, fake_data, original_data ,labels, opt_discr)
    ps = params(discr.aux, discr.discr,discr.model)
    loss = 0.0
    for i in 1:2
    loss, back = Zygote.pullback(ps) do
           original_common = discr.model(original_data) |> gpu
           fake_common =  discr.model(fake_data) |> gpu
           d1 = discr_loss(discr.discr(original_common), discr.discr(fake_common))
           x = discr.aux(original_common)
           a1 = aux_loss(x, labels |> gpu)
           d1+ a1
           end
    grads = back(1f0)
    update!(opt_discr, ps, grads)
    end
    return loss
end

Zygote.@nograd train_discr

function train_gan(gen, discr, original_data, label, opt_gen, opt_discr, hparams)
    noise = randn!(similar(original_data, (hparams.latent_dim, hparams.batch_size))) |> gpu
    ps = params(gen.m1,gen.m2, gen.model)
    loss = Dict()
    loss["gen"], back = Zygote.pullback(ps) do
            fake = gen(label, noise) |> gpu
            loss["discr"] = train_discr(discr, fake, original_data, label, opt_discr)
            fake_common = discr.model(fake)
            generator_loss(discr.discr(fake_common)) + aux_loss(discr.aux(fake_common), label)
    end
    grads = back(1f0)
    update!(opt_gen, ps, grads)
    return loss
end

function create_output_image(gen, fixed_noise, fixed_labels, hparams)
    @eval Flux.istraining() = false
    fake_images = @. cpu(gen(fixed_labels, fixed_noise))
    @eval Flux.istraining() = true
    image_array = dropdims(reduce(vcat, reduce.(hcat, partition(fake_images, hparams.output_y))); dims=(3, 4))
    image_array = @. Gray(image_array + 1f0) / 2f0
    return image_array
end

function train()
    hparams = HyperParams()

    data = load_data(hparams)

    fixed_noise = [randn(hparams.latent_dim, 1) |> gpu for _=1:hparams.output_x*hparams.output_y]

    fixed_labels = [float.(Flux.onehotbatch(rand(0:hparams.nclasses-1,1), 0:hparams.nclasses-1)) |> gpu for _=1:hparams.output_x*hparams.output_y]

    # Discriminator
    dscr = discriminator(hparams)

    # Generator
    gen =  generator(hparams)

    # Optimizers
    opt_dscr = ADAM(hparams.lr_dscr, (0.5,0.999))
    opt_gen = ADAM(hparams.lr_gen, (0.5,0.999))

    isdir("output")||mkdir("output")

    # Training
    train_steps = 0
    for ep in 1:hparams.epochs
        @info "Epoch $ep"
        for (x,y) in data[1:468]
            # Update discriminator and generator
            loss = train_gan(gen, dscr, x, y, opt_gen, opt_dscr, hparams)

            if train_steps % hparams.verbose_freq == 0
                @info("Train step $(train_steps), Discriminator loss = $(loss["discr"]), Generator loss = $(loss["gen"])")
                # Save generated fake image
                output_image = create_output_image(gen, fixed_noise, fixed_labels, hparams)
                save(@sprintf("output/cgan_steps_%06d.png", train_steps), output_image)
            end
            train_steps += 1
        end
    end

    output_image = create_output_image(gen, fixed_noise, fixed_labels, hparams)
    save(@sprintf("output/cgan_steps_%06d.png", train_steps), output_image)
    return fixed_labels
end

cd(@__DIR__)
fixed_labels = train()


text:
"""
2D traveling wave solution for the Euler equations in a parallelogram domain

     --------
    /      /
   /      /
  /      /
 --------

δx = 1, δy = 0.5, θ = π / 4

"""

using KitBase, FluxReconstruction, LinearAlgebra, OrdinaryDiffEq, Plots
using KitBase.OffsetArrays
using KitBase.ProgressMeter: @showprogress
cd(@__DIR__)
pyplot()

set = Setup(
    case = "dev",
    space = "2d0f0v",
    flux = "hll",
    collision = "nothing",
    interpOrder = 3,
    limiter = "positivity",
    boundary = "euler",
    cfl = 0.1,
    maxTime = 1.0,
)

begin
    lx, ly = 1.0, 0.5
    x0 = 0.0
    x1 = x0 + lx
    nx = 100
    Δx = lx / nx
    y0 = 0.0
    y1 = y0 + ly
    ny = 50
    Δy = ly / ny

    vertices = OffsetArray{Float64}(undef, 0:nx+1, 0:ny+1, 4, 2)
    for j in axes(vertices, 2), i in axes(vertices, 1)
        vertices[i, j, 1, 1] = x0 + 0.5 * (j - 1) / ny + (i - 1) * Δx
        vertices[i, j, 2, 1] = vertices[i, j, 1, 1] + Δx
        vertices[i, j, 3, 1] = vertices[i, j, 2, 1] + Δy
        vertices[i, j, 4, 1] = vertices[i, j, 3, 1] - Δx

        vertices[i, j, 1, 2] = y0 + 0.5 * (j - 1) * Δy
        vertices[i, j, 2, 2] = vertices[i, j, 1, 2]
        vertices[i, j, 3, 2] = vertices[i, j, 1, 2] + Δy
        vertices[i, j, 4, 2] = vertices[i, j, 3, 2]
    end

    x = OffsetArray{Float64}(undef, 0:nx+1, 0:ny+1)
    y = similar(x)
    dx = similar(x)
    dy = similar(y)
    for i = 0:nx+1, j = 0:ny+1
        x[i, j] = (vertices[i, j, 1, 1] + vertices[i, j, 4, 1]) / 2 + Δx / 2
        y[i, j] = vertices[i, j, 1, 2] + 0.5 * Δy
        dx[i, j] = Δx
        dy[i, j] = Δy
    end

    ps = PSpace2D(x0, x1, nx, y0, y1, ny, x, y, dx, dy, vertices)
end

vs = nothing
gas = Gas(Kn = 1e-6, Ma = 0.0, K = 1.0)
fw = function (x...)
    zeros(4)
end
ib = IB(fw, gas)

ks = SolverSet(set, ps, vs, gas, ib)

ctr, a1face, a2face = init_fvm(ks)
for i in axes(ctr, 1), j in axes(ctr, 2)
    ρ = 1.0 + 0.1 * sin(2π * i / ps.nx)
    ctr[i, j].prim .= [ρ, 1.0, 0.0, ρ]
    ctr[i, j].w .= prim_conserve(ctr[i, j].prim, ks.gas.γ)
end

function get_solution(ks, ctr)
    sol = zeros(ks.ps.nx, ks.ps.ny, 4)
    for i = 1:ks.ps.nx, j = 1:ks.ps.ny
        sol[i, j, :] .= conserve_prim(ctr[i, j].w, ks.gas.γ)
        sol[i, j, 4] = 1 / sol[i, j, 4]
    end
    return sol
end

sol0 = get_solution(ks, ctr)

tspan = (0.0, 1.0)
dt = 0.001
nt = tspan[2] ÷ dt |> Int
t = 0.0

@showprogress for iter = 1:nt÷2
    evolve!(
        ks,
        ctr,
        a1face,
        a2face,
        dt;
        mode = :hll,
        bc = [:period, :period, :extra, :extra],
    )
    update!(
        ks,
        ctr,
        a1face,
        a2face,
        dt,
        zeros(4);
        bc = [:period, :period, :extra, :extra],
    )

    t += dt
end

sol = get_solution(ks, ctr)

begin
    idx = 4
    plot(ps.x[1:ps.nx, 1], sol0[1:ps.nx, 1, 1])
    plot!(ps.x[1:ps.nx, 1], sol[1:ps.nx, 1, 1])
end
contourf(ps.x[1:ps.nx, 1:ps.ny], ps.y[1:ps.nx, 1:ps.ny], sol[1:ps.nx, 1:ps.ny, 2])


text:
export NLLob, PY, CY, NY

const PY, CY, NY = 1, 2, 3   ## previous, current, next year

"Line of Business (Non-Life)"
mutable struct NLLob
  "`AbstractString`: Name of the line of business"
  name::Symbol
  "`Int`: Identifying index of the line of business"
  index::Int
  "`Vector{Float64}`: Written gross premium"
  prem_gross_w::Vector{Float64}
  "`Vector{Float64}`: Written premium"
  prem_w::Vector{Float64}
  "`Float64`: Gross premium earned"
  prem_gross_cy::Float64
  "`Vector{Float64}`: Premium earned"
  prem::Vector{Float64}
  "`Vector{Float64}`: Unearned gross prem. reserve"
  upr_gross::Vector{Float64}
  "`Vector{Float64}`: Unearned premium reserve"
  upr::Vector{Float64}
  "`Vector{Float64}`: Prop. reinsurence ceded"
  re_prop_q::Vector{Float64}
  "`Vector{Float64}`: Reserve pattern"
  β::Vector{Float64}
  "`Float64`:  Volume factor premium"
  vol_prem::Float64
  "`Float64`: Var. coeff. net prem risk"
  vol_prem_vc::Float64
  "`Float64`: Volume factor reserve"
  vol_res::Float64
  "`Float64`: Standard deviation prem+res"
  vol_prem_res_sd::Float64
end


text:
"""
Demonstrates optimizing a network with the Metaheuristics package. Other possibilities include using the package: Evolutionary, evolutionary algorithm package.

"""
using Plots
using SpikeNetOpt
SNO = SpikeNetOpt
using SpikingNeuralNetworks
using Evolutionary
using Memoize
SNN = SpikingNeuralNetworks
SNN.@load_units

##
# Ground truths
##
const Ne = 200;
const Ni = 50
const σee = 1.0
const pee = 0.5
const σei = 1.0
const pei = 0.5
MU = 10

const E
const spkd_ground

##
# The memorize macro/decorator caches big data values.
##
@memoize function get_constant_gw()
    ##
    # Load a connectome from disk.
    ##
    try
        filename = string("P1C1_weights.jld")
        P = load(filename,"P")
        C = load(filename,"C")
        return P, C

    catch e
        P, C = SNO.make_net_SNN(Ne, Ni, σee = 0.5, pee = 0.8, σei = 0.5, pei = 0.8)
        filename = string("P1C1_weights.jld")
        save(filename, "P", P)
        save(filename, "C", C)
        return P, C
    end
end

P, C = get_constant_gw()
E, I = P
EE, EI, IE, II = C
SNN.monitor([E, I], [:fire])
sim_length = 1000
@inbounds for t = 1:sim_length*ms
    E.I = vec([11.5 for i = 1:sim_length])#vec(E_stim[t,:])#[i]#3randn(Ne)
    SNN.sim!(P, C, 1ms)

end
spkd_ground = SpikeNetOpt.get_trains(P[1])
display(SNN.raster(P[1]))
SNN.raster(P[1]) |> display

function loss(weight_gain_factor)
    @show(Ne,Ni)
    @show(model)

    P1, C1 = weight_gain_factor[1]
    for c in C1
        c.W = c.W.*weight_gain_factor

    E1, I1 = P1
    SNN.monitor([E1, I1], [:fire])
    sim_length = 1000
    @inbounds for t = 1:sim_length*ms
        E1.I = vec([11.5 for i = 1:sim_length])
        SNN.sim!(P1, C1, 1ms)
    end

    spkd_found = get_trains(P1[1])
    println("Ground Truth \n")
    SNN.raster([E]) |> display
    println("Best Candidate \n")
    SNN.raster([E1]) |> display
    error = SNO.spike_train_difference(spkd_ground, spkd_found)
    error = sum(error)
    @show(error)
    error
end

D = 10
bounds = [3ones(D) 40ones(D)]'
a = view(bounds, 1, 1)
b = view(bounds, 1, 2)
information = Information(f_optimum = 0.0)
options = Options( seed = 1, iterations=10, f_calls_limit =10)

D = size(bounds, 2)
nobjectives=1

options = Options( seed = 1, iterations=10000, f_calls_limit = 25000)
npartitions = nobjectives == 2 ? 100 : 12

methods = [
        SMS_EMOA(N = 5, n_samples=5, options=options),
        NSGA2(options=options)
        ]

for method in methods
    f_calls = 0
    result = optimize(SpikeNetOpt.loss, bounds, method)
    #result = optimize(SpikeNetOpt.loss, bounds)#, method)
    @show(result)

end


text:
"""
StableFluids.jl

See https://github.com/Ceyron/StableFluids.jl#readme

Functions
* `StableFluids.run_2D()`: Run 2D simulation
* `StableFluids.run_3D()`: Run 3D simulation

"""
module StableFluids
    include("TwoDims.jl")
    include("ThreeDims.jl")

    run_2D() = TwoDims.main()
    run_3D() = ThreeDims.main()
end

max_stars_repo_name: Ceyron/stable-fluids-julia

max_issues_repo_name: Ceyron/stable-fluids-julia

max_forks_repo_name: Ceyron/stable-fluids-julia


text:
#############################################################################
#############################################################################
#
# This file contains units tests for polynomial operations
#
#############################################################################
#############################################################################

"""
Test product of polynomials.
"""
function prod_test_poly(;N::Int = 10^3, N_prods::Int = 20, seed::Int = 0)
    Random.seed!(seed)
    for _ in 1:N
        p1 = rand(Polynomial)
        p2 = rand(Polynomial)
        prod = p1*p2
        @assert leading(prod) == leading(p1)*leading(p2)
    end

    for _ in 1:N
        p_base = Polynomial(Term(1,0))
        for _ in 1:N_prods
            p = rand(Polynomial)
            prod = p_base*p
            @assert leading(prod) == leading(p_base)*leading(p)
            p_base = prod
        end
    end
    println("prod_test_poly - PASSED")
end

"""
Test derivative of polynomials (as well as product).
"""
function prod_derivative_test_poly(;N::Int = 10^2,  seed::Int = 0)
    Random.seed!(seed)
    for _ in 1:N
        p1 = rand(Polynomial)
        p2 = rand(Polynomial)
        p1d = derivative(p1)
        p2d = derivative(p2)
        @assert (p1d*p2) + (p1*p2d) == derivative(p1*p2)
    end
    println("prod_derivative_test_poly - PASSED")
end

"""
Test division of polynomials modulo p.
"""
function division_test_poly(;prime::Int = 101, N::Int = 10^4, seed::Int = 0)
    Random.seed!(seed)
    for _ in 1:N
        p1 = rand(Polynomial)
        p2 = rand(Polynomial)
        p_prod = p1*p2
        q, r = Polynomial(), Polynomial()
        try
            q, r = divide(p_prod, p2)(prime)
            if (q, r) == (nothing,nothing)
                println("Unlucky prime: $p1 is reduced to $(p1 % prime) modulo $prime")
                continue
            end
        catch e
            if typeof(e) == DivideError
                @assert mod(p2, prime) == 0
            else
                throw(e)
            end
        end
        @assert iszero( mod(q*p2+r - p_prod, prime) )
    end
    println("division_test_poly - PASSED")
end

"""
Test the extended euclid algorithm for polynomials modulo p.
"""
function ext_euclid_test_poly(;prime::Int=101, N::Int = 10^3, seed::Int = 0)
    Random.seed!(seed)
    for _ in 1:N
        p1 = rand(Polynomial)
        p2 = rand(Polynomial)
        g, s, t = extended_euclid_alg(p1, p2, prime)
        @assert mod(s*p1 + t*p2 - g, prime) == 0
    end
    println("ext_euclid_test_poly - PASSED")
end


text:
module PatModules

export @once, @mainmodule, @auxmodule, @mainmodule!, @auxmodule!

# Adapted from https://github.com/JuliaLang/julia/issues/29966
macro once(m)
    @assert m.head === :module
    name = m.args[2]
    gvname = Symbol("#@once#gv#for#$name")
    quote
        if !isdefined($__module__, $(QuoteNode(gvname)))
            global $(esc(gvname)) = 1
            Core.eval($__module__, $(QuoteNode(m)))
        end
        nothing
    end
end

function _include(modul, to_import)
    for importname ∈ to_import
        if isdir(importname)
            _, filename = splitdir(importname)
            Base.include(modul, joinpath(importname, filename) * ".jl")
        else
            Base.include(modul, importname * ".jl")
        end
    end
end

function _module(name::Symbol, to_import_expr::Expr, block::Expr, __module__, __source__; is_aux::Bool, is_use::Bool)
    to_import::Tuple{Vararg{String}} = map(String, eval(to_import_expr))
    for importname ∈ to_import
        @assert !occursin("/", importname)
        @assert !occursin("\\", importname)
        @assert !occursin(".", importname)
    end
    path = String(__source__.file)
    path = abspath(path)
    path = realpath(path)
    path, _ = splitdir(path)
    to_include = map((importname)->joinpath(path, importname), to_import)

    m = quote
        if $is_aux
            $_include($__module__, $to_include)
        end
        @once module $name
            if !$is_aux
                $_include($name, $to_include)
            end
            for importname ∈ $to_import
                importname_symbol = Symbol(importname)
                if $is_aux
                    if $is_use
                        Core.eval($name, :(using .. $importname_symbol))
                    else
                        Core.eval($name, :(import .. $importname_symbol))
                    end
                else
                    if $is_use
                        Core.eval($name, :(using .$importname_symbol))
                    else
                        Core.eval($name, :(import .$importname_symbol))
                    end
                end
            end
            $(block.args...)
        end
    end
    quote
        using PatModules
        Core.eval($__module__, $(QuoteNode(m)))
    end
end

macro mainmodule(name::Symbol, to_import_expr::Expr, block::Expr)
    _module(name, to_import_expr, block, __module__, __source__; is_aux=false, is_use=false)
end

macro auxmodule(name::Symbol, to_import_expr::Expr, block::Expr)
    _module(name, to_import_expr, block, __module__, __source__; is_aux=true, is_use=false)
end

macro mainmodule!(name::Symbol, to_import_expr::Expr, block::Expr)
    _module(name, to_import_expr, block, __module__, __source__; is_aux=false, is_use=true)
end

macro auxmodule!(name::Symbol, to_import_expr::Expr, block::Expr)
    _module(name, to_import_expr, block, __module__, __source__; is_aux=true, is_use=true)
end

end


text:
export BMSOM
export Online, Batch
export train!, labeling!, predict

abstract UpdateMode

type Online <: UpdateMode
end
type Batch <: UpdateMode
end

type Block
    pos :: Tuple{Int,Int}
    indices :: Tuple
    blksize :: Int
    distance :: Number
    function Block(mapsize::Tuple, pos::Tuple{Int,Int}, blksize::Int)
        indices = (:, vec(Int[sub2ind(mapsize, i,j)
            for i=pos[1]:pos[1]+blksize-1,
                j=pos[2]:pos[2]+blksize-1]))
        new(pos, indices, blksize)
    end
end

type BMSOM{T<:AbstractFloat}
    mapsize :: Int
    refvecs :: Array{T,3}
    labels :: Array{Integer,2}
    metric :: Distances.PreMetric
    labeldict :: Array{Any}
    labelhist :: Array{T,3}

    function BMSOM(mapsize::Int, dim::Int;
                   init::Initializer=UniformInit(), metric::Distances.PreMetric=Euclidean())

        refvecs = initialize(T, (dim, mapsize, mapsize), init)
        labels = Array{Integer,2}(mapsize, mapsize)

        new(mapsize, refvecs, labels, metric)
    end
end

function Base.getindex(som::BMSOM, index...)
    som.refvecs[index...]
end

function Base.ndims(som::BMSOM)
    size(som.refvecs,1)
end

function Base.size(som::BMSOM)
    (som.mapsize, som.mapsize)
end

function calcdist!(blk::Block, iv::AbstractArray, refvecs::AbstractArray,  metric::Distances.PreMetric)
    mv = mean(sub(refvecs, blk.indices), 2)
    blk.distance = evaluate(metric, vec(mv), iv)
end

function update_loop!(refvecs::AbstractArray, drefvecs::AbstractArray, tv::AbstractArray, blk::Block, umode::Batch)
    blkrvs = sub(refvecs, blk.indices)
    dblkrvs = sub(drefvecs, blk.indices)
    d = (tv*ones(1, size(blkrvs,2)) - blkrvs)
    dblkrvs[:,:] += d/blk.blksize^2
end

function update_loop!(refvecs::AbstractArray, drefvecs::AbstractArray, tv::AbstractArray, blk::Block, umode::Online)
    blkrvs = sub(refvecs, blk.indices)
    d = (tv*ones(1, size(blkrvs,2)) - blkrvs)
    blkrvs[:,:] += d/blk.blksize^2
end

function update_after!{T}(refvecs::AbstractArray{T,3}, drefvecs::AbstractArray{T,3}, umode::Batch)
    refvecs[:,:,:] += drefvecs
    drefvecs[:] = zero(T)
end

function update_after!(refvecs::AbstractArray, drefvecs::AbstractArray, umode::Online)
    # do nothing
end

function findwinner(som::BMSOM, tv::AbstractArray, minblksize::Int=2)
    x,y = 1,1
    winblk = Block(size(som), (x,y), som.mapsize)
    calcdist!(winblk, tv, som.refvecs, som.metric)
    for blksize = som.mapsize-1:-1:minblksize
        wcands = Block[Block(size(som), (x+i,y+j), blksize) for i=0:1, j=0:1]
        ds = map(blk->calcdist!(blk, tv, som.refvecs, som.metric), wcands)
        val,idx = findmin(ds)
        if val < winblk.distance
            winblk = wcands[idx]
        end
        x,y = winblk.pos
    end

    return winblk
end

function train!{T<:AbstractFloat}(som::BMSOM{T}, tvs::AbstractArray{T,2}, epoch::Int;
                           minblksize::Int=2, umode::UpdateMode=Batch(), control::Bool=true)
    d_refvecs = zeros(T, size(som.refvecs))
    blksize = minblksize
    for i=1:epoch
        if control
            blksize = Int(floor(som.mapsize-(som.mapsize-(minblksize))/epoch*i))
        end
        for di=1:size(tvs,2)
            winblk = findwinner(som, tvs[:,di], blksize)
            update_loop!(som.refvecs, d_refvecs, tvs[:,di], winblk, umode)
        end
        update_after!(som.refvecs, d_refvecs, umode)
    end
end

function labeling!{T<:AbstractFloat, S<:Any}(som::BMSOM{T}, tvs::AbstractArray{T,2}, tvls::AbstractArray{S,1})
    @assert(size(tvs,2) == size(tvls,1), "training data size error.")

    som.labeldict = unique(tvls)
    nlabels = length(som.labeldict)

    lbl2idx = Dict{S, Integer}()
    for i=1:nlabels
        lbl2idx[som.labeldict[i]] = i
    end

    som.labelhist = zeros(T, (nlabels, som.mapsize, som.mapsize))

    for j=1:som.mapsize, i=1:som.mapsize
        for di=1:size(tvs,2)
            som.labelhist[lbl2idx[tvls[di]],i,j] += evaluate(som.metric, som.refvecs[:,i,j], tvs[:,di])
        end
    end

    for i=1:som.mapsize,j=1:som.mapsize
        val, idx = findmin(som.labelhist[:,i,j])
        som.labels[i,j] = idx
    end
end

function predict{T<:AbstractFloat}(som::BMSOM{T}, iv::AbstractArray{T})
    dist = zeros(T, size(som))
    for j=1:som.mapsize, i=1:som.mapsize
        dist[i,j] = evaluate(som.metric, vec(som.refvecs[:,i,j]), iv)
    end
    val, idx = findmin(dist)

    som.labeldict[som.labels[idx]], val
end


text:
# author: Tommaso Sartor

module riccati

using LinearAlgebra
using MKL

# Riccat recursion

function riccati_trf(N, nx, nu, BAt, RSQ, L, LN, BAtL, M)

	LN = RSQ[nu+1:nu+nx, nu+1:nu+nx];
	LAPACK.potrf!('L', LN);

	for ii in 1:N
		BAtL = copy(BAt);
		if ii==1
#			BLAS.gemm!('N', 'N', 1.0, BAt, LN, 0.0, BAtL);
			BLAS.trmm!('R', 'L', 'N', 'N', 1.0, LN, BAtL);
		else
#			BLAS.gemm!('N', 'N', 1.0, BAt, L[nu+1:nu+nx,nu+1:nu+nx,N+2-ii], 0.0, BAtL);
#			display(L[:,:,N+2-ii])
#			println()
#			BLAS.trmm!('R', 'L', 'N', 'N', 1.0, L[nu+1:nu+nx,nu+1:nu+nx,N+2-ii], BAtL);
			BLAS.trmm!('R', 'L', 'N', 'N', 1.0, view(L, nu+1:nu+nx, nu+1:nu+nx, N+2-ii), BAtL);
		end
#		M = copy(RSQ);
#		BLAS.syrk!('L', 'N', 1.0, BAtL, 1.0, M);
#		LAPACK.potrf!('L', M);
#		L[:,:,N+1-ii] = copy(M);
		L[:,:,N+1-ii] .= RSQ;
		MM = view(L, :, :, N+1-ii);
		BLAS.syrk!('L', 'N', 1.0, BAtL, 1.0, MM);
		LAPACK.potrf!('L', MM);
	end

#display(L)
#println()

end

export riccati_trf

end

max_stars_repo_name: kokizzu/prometeo

max_issues_repo_name: kokizzu/prometeo

max_forks_repo_name: kokizzu/prometeo


text:
@generated function iszeroinf{N,ES}(x::Posit{N,ES})
  IntType = __UT(Posit{N,ES})
  otherbits = (one(IntType) << (N - 1) - one(IntType))
  z = zero(IntType)
  :(__u(x) & $otherbits == $z)
end

@generated function Base.isinf{N,ES}(x::Posit{N,ES})
  IntType = __UT(Posit{N,ES})
  topbit = one(IntType) << (N - 1)
  :(__u(x) == $topbit)
end
#this code may need to change if we activate NaN mode.
Base.isfinite{N,ES}(x::Posit{N,ES}) = !isinf(x)
Base.isnan{N,ES}(x::Posit{N,ES}) = false

Base.isequal{N,ES}(lhs::Posit{N,ES}, rhs::Posit{N,ES}) = __u(lhs) == __u(rhs)
Base.isless{N,ES}(lhs::Posit{N,ES}, rhs::Posit{N,ES}) = __s(lhs) < __s(rhs)

#zero and one constants
@generated function Base.zero{T <: Posit}(::Type{T})
  IntType = __UT(T)
  PositType = __PT(T)
  :(reinterpret($PositType, zero($IntType)))
end
@generated function Base.zero{N,ES}(::Posit{N,ES})
  IntType = __UT(Posit{N,ES})
  PositType = __PT(Posit{N,ES})
  :(reinterpret($PositType, zero($IntType)))
end

@generated function Base.one{T <: Posit}(::Type{T})
  IntType = __UT(T)
  PositType = __PT(T)
  N = sizeof(T) * 8
  :(reinterpret($PositType, one($IntType) << ($N - 2)))
end
@generated function Base.one{N,ES}(::Posit{N,ES})
  IntType = __UT(Posit{N,ES})
  PositType = __PT(Posit{N,ES})
  :(reinterpret($PositType, one($IntType) << (N - 2)))
end

#spoof the sizeof operator to do the right thing here.
Base.sizeof{N,ES}(::Type{Posit{N,ES}}) = div(N,8)
Base.sizeof{N}(::Type{Posit{N}}) = div(N,8)


text:
#!/usr/bin/env bash
    #=
    exec julia --project="$(realpath $(dirname $0))/../" "${BASH_SOURCE[0]}" "$@" -e 'include(popfirst!(ARGS))' \
    "${BASH_SOURCE[0]}" "$@"
    =#

# TODO: select optimal threshold for each feature
# TODO: attentional cascading

using Base.Threads: @threads
using Base.Iterators: partition
using ProgressMeter: @showprogress, Progress, next!

function β(i::T)::T where T
    return @fastmath(T(0.5) * log((one(i) - i) / i))
end

function get_feature_votes(
    positive_path::AbstractString,
    negative_path::AbstractString,
    num_classifiers::Integer=-one(Int32),
    min_feature_width::Integer=one(Int32),
    max_feature_width::Integer=-one(Int32),
    min_feature_height::Integer=one(Int32),
    max_feature_height::Integer=-one(Int32);
    scale::Bool = false,
    scale_to::Tuple = (Int32(200), Int32(200))
    )

    #this transforms everything to maintain type stability
    s₁, s₂ = scale_to
    min_feature_width,
    max_feature_width,
    min_feature_height,
    max_feature_height, s₁, s₂ = promote(
        min_feature_width,
        max_feature_width,
        min_feature_height,
        max_feature_height,
        s₁,
        s₂
    )
    scale_to = (s₁, s₂)
    _Int = typeof(max_feature_width)

    # get number of positive and negative image
    positive_files = filtered_ls(positive_path)
    negative_files = filtered_ls(negative_path)
    image_files = vcat(positive_files, negative_files)
    num_pos = length(positive_files)
    num_neg = length(negative_files)
    num_imgs = num_pos + num_neg

    # get image height and width
    temp_image = load_image(rand(positive_files), scale=scale, scale_to=scale_to)
    img_height, img_width = size(temp_image)
    temp_image = nothing # unload temporary image

    # Maximum feature width and height default to image width and height
    max_feature_height = isequal(max_feature_height, -_Int(1)) ? img_height : max_feature_height
    max_feature_width = isequal(max_feature_width, -_Int(1)) ? img_height : max_feature_width

    # Create features for all sizes and locations
    features = create_features(img_height, img_width, min_feature_width, max_feature_width, min_feature_height, max_feature_height)
    num_features = length(features)
    num_classifiers = isequal(num_classifiers, -_Int(1)) ? num_features : num_classifiers

    # create an empty array with dimensions (num_imgs, numFeautures)
    votes = Matrix{Int8}(undef, num_features, num_imgs)

    notify_user("Loading images ($(num_pos) positive and $(num_neg) negative images) and calculating their scores...")
    p = Progress(length(image_files), 1) # minimum update interval: 1 second
    num_processed = 0
    batch_size = 10
    # get votes for images
    map(partition(image_files, batch_size)) do batch
        ii_imgs = load_image.(batch; scale=scale, scale_to=scale_to)
        @threads for t in 1:length(batch)
            votes[:, num_processed + t] .= get_vote.(features, Ref(ii_imgs[t]))
            # map!(f -> get_vote(f, ii_imgs[t]), view(votes, :, num_processed + t), features)
            next!(p) # increment progress bar
        end
        num_processed += length(batch)
    end
    print("\n") # for a new line after the progress bar

    return votes, features
end

function learn(
    positive_path::AbstractString,
    negative_path::AbstractString,
    features::Array{HaarLikeObject, 1},
    votes::Matrix{Int8},
    num_classifiers::Integer=-one(Int32)
    )
    # get number of positive and negative images (and create a global variable of the total number of images——global for the @everywhere scope)
    num_pos = length(filtered_ls(positive_path))
    num_neg = length(filtered_ls(negative_path))
    num_imgs = num_pos + num_neg

    # Initialise weights $w_{1,i} = \frac{1}{2m}, \frac{1}{2l}$, for $y_i=0,1$ for negative and positive examples respectively
    pos_weights = ones(num_pos) / (2 * num_pos)
    neg_weights = ones(num_neg) / (2 * num_neg)

    # Concatenate positive and negative weights into one `weights` array
    weights = vcat(pos_weights, neg_weights)
    labels = vcat(ones(Int8, num_pos), ones(Int8, num_neg) * -one(Int8))

    # get number of features
    num_features = length(features)
    feature_indices = Array(1:num_features)
    num_classifiers = isequal(num_classifiers, -1) ? num_features : num_classifiers

    # select classifiers
    notify_user("Selecting classifiers...")
    classifiers = HaarLikeObject[]
    p = Progress(num_classifiers, 1) # minimum update interval: 1 second
    classification_errors = Vector{Float64}(undef, length(feature_indices))

    for t in 1:num_classifiers
        # normalize the weights $w_{t,i}\gets \frac{w_{t,i}}{\sum_{j=1}^n w_{t,j}}$
        weights .*= inv(sum(weights))

        # For each feature j, train a classifier $h_j$ which is restricted to using a single feature.  The error is evaluated with respect to $w_j,\varepsilon_j = \sum_i w_i\left|h_j\left(x_i\right)-y_i\right|$
        @threads for j in 1:length(feature_indices)
            classification_errors[j] = sum(1:num_imgs) do img_idx
                labels[img_idx] !== votes[feature_indices[j], img_idx] ? weights[img_idx] : zero(Float64)
            end
        end

        # choose the classifier $h_t$ with the lowest error $\varepsilon_t$
        best_error, min_error_idx = findmin(classification_errors)
        best_feature_idx = feature_indices[min_error_idx]
        best_feature = features[best_feature_idx]

        # set feature weight
        best_feature = features[best_feature_idx]
        feature_weight = β(best_error)
        best_feature.weight = feature_weight

        # append selected features
        classifiers = push!(classifiers, best_feature)

        # update image weights $w_{t+1,i}=w_{t,i}\beta_{t}^{1-e_i}$
        sqrt_best_error = @fastmath(sqrt(best_error / (one(best_error) - best_error)))
        inv_sqrt_best_error = @fastmath(sqrt((one(best_error) - best_error)/best_error))
        @inbounds for i in 1:num_imgs
            if labels[i] !== votes[best_feature_idx, i]
                weights[i] *= inv_sqrt_best_error
            else
                weights[i] *= sqrt_best_error
            end
        end

        # remove feature (a feature can't be selected twice)
        filter!(e -> e ∉ best_feature_idx, feature_indices) # note: without unicode operators, `e ∉ [a, b]` is `!(e in [a, b])`
        resize!(classification_errors, length(feature_indices))
        next!(p) # increment progress bar
    end

    print("\n") # for a new line after the progress bar

    return classifiers

end

function learn(
    positive_path::AbstractString,
    negative_path::AbstractString,
    num_classifiers::Int64=-1,
    min_feature_width::Int64=1,
    max_feature_width::Int64=-1,
    min_feature_height::Int64=1,
    max_feature_height::Int64=-1;
    scale::Bool = false,
    scale_to::Tuple = (200, 200)
)

    votes, features = get_feature_votes(
        positive_path,
        negative_path,
        num_classifiers,
        min_feature_width,
        max_feature_width,
        min_feature_height,
        max_feature_height,
        scale = scale,
        scale_to = scale_to
    )

    return learn(positive_path, negative_path, features, votes, num_classifiers)
end

"""
    create_features(
        img_height::Integer,
        img_width::Integer,
        min_feature_width::Integer,
        max_feature_width::Integer,
        min_feature_height::Integer,
        max_feature_height::Integer
    ) -> Array{HaarLikeObject, 1}

Iteratively creates the Haar-like feautures

# Arguments

- `img_height::Integer`: The height of the image
- `img_width::Integer`: The width of the image
- `min_feature_width::Integer`: The minimum width of the feature (used for computation efficiency purposes)
- `max_feature_width::Integer`: The maximum width of the feature
- `min_feature_height::Integer`: The minimum height of the feature
- `max_feature_height::Integer`: The maximum height of the feature

# Returns

- `features::AbstractArray`: an array of Haar-like features found for an image
"""
function create_features(
    img_height::Int64,
    img_width::Int64,
    min_feature_width::Int64,
    max_feature_width::Int64,
    min_feature_height::Int64,
    max_feature_height::Int64
)
    notify_user("Creating Haar-like features...")
    features = HaarLikeObject[]

    if img_width < max_feature_width || img_height < max_feature_height
        error("""
        Cannot possibly find classifiers whose size is greater than the image itself [(width,height) = ($img_width,$img_height)].
        """)
    end

    for feature in values(feature_types) # (feature_types are just tuples)
        feature_start_width = max(min_feature_width, first(feature))
        for feature_width in feature_start_width:first(feature):(max_feature_width)
            feature_start_height = max(min_feature_height, last(feature))
            for feature_height in feature_start_height:last(feature):(max_feature_height)
                for x in 1:(img_width - feature_width)
                    for y in 1:(img_height - feature_height)
                        push!(features, HaarLikeObject(feature, (x, y), feature_width, feature_height, 0, 1))
                        push!(features, HaarLikeObject(feature, (x, y), feature_width, feature_height, 0, -1))
                    end # end for y
                end # end for x
            end # end for feature height
        end # end for feature width
    end # end for feature in feature types

    println("...finished processing; ", length(features), " features created.\n")

    return features
end