text:
"""
`FinancialPortfolios`

A minimalist Julia package for working with simple portfolios of financial assets.
"""
module FinancialPortfolios

export FinancialPortfolio
export portfolioreturn, update!
export positions

"""
`FinancialPortfolio`

A (possibly named) collection of portfolio weights.
Weights can be a dictionary-like object matching asset identifiers with weights or just a vector of weights.

_**Example**_
```
w = [0.5, 0.25, 0.25]
FP_vec = FinancialPortfolio(w)

nm = ["a", "b", "c"]
plain_dict = Dict(Iterators.zip(nm,w))
FP_plaindict = FinancialPortfolio(plain_dict)

using Dictionaries
fancy_dict = dictionary(plain_dict)
FP_fancydict = FinancialPortfolio(fancy_dict)
```
"""
struct FinancialPortfolio{T}
    positions::T
    FinancialPortfolio{T}(x) where {T} = new(x)
end
# T should be <:AbstractDict or Dictionaries.AbstractDictionary or AbstractVector
# NamedTuple works except is immutable so breaks update!
FinancialPortfolio(x::T) where {T} = FinancialPortfolio{T}(copyandnormalize(x))

"""
`weights(fp::FinancialPortfolio)`

Recover the weights of the assets in the portfolio (as a fraction of total value).
"""
positions(fp::FinancialPortfolio) = fp.positions

"""
`portfolioreturn(fp::FinancialPortfolio,ret)`

Compute the weighted portfolio return. If the portfolio weights have names or identifiers, the returns should as well.
The portfolio is allowed to contain only a subset of the supplied asset returns, but not the other way around.

_**Example**_
```
w = [0.5, 0.25, 0.25]
r = [0.1, 0.1, -0.1]
FP = FinancialPortfolio(w)
portfolioreturn(FP,r)
```
"""
function portfolioreturn(fp::FinancialPortfolio,ret)
    portreturn = 0.0
    for k in eachindex(fp.positions)
        portreturn += fp.positions[k] * ret[k]
    end
    return portreturn
end

"""
`update!(fp::FinancialPortfolio,ret)`

Updates the weights of the portfolio in place. Returns the period portfolio return.

The portfolio is allowed to contain only a subset of the supplied asset returns, but not the other way around.

_**Example**_
```
w = [0.5, 0.25, 0.25]
r = [0.1, 0.1, -0.1]
FP = FinancialPortfolio(w)
update!(FP,r)
FP
```
"""
function update!(fp::FinancialPortfolio,ret)
    portreturn = portfolioreturn(fp,ret)

    for k in eachindex(fp.positions)
        fp.positions[k] = fp.positions[k] * (1+ret[k])
    end
    normalize!(fp)
    return portreturn
end

##################### UNEXPORTED AND/OR EXTENDED BASE METHODS
"""
`FinancialPortfolios.normalize!(fp::FinancialPortfolio)`

Normalize portfolio weights to sum to 1. Used during [`update!`](@ref).
"""
function normalize!(fp::FinancialPortfolio)
    fac = sum(values(fp.positions))
    for k in eachindex(fp.positions)
        fp.positions[k] = fp.positions[k] / fac
    end
    return nothing
end

"""
`FinancialPortfolios.copyandnormalize(x)`

Normalize portfolio weights to sum to 1. Used during construction.
"""
function copyandnormalize(x0)
    x = copy(x0)
    fac = sum(values(x))
    for k in eachindex(x)
        x[k] = x[k] / fac
    end
    return x
end

"""
`keys(fp::FinancialPortfolio)`

Recover the names or identifiers of the assets in the portfolio.
"""
Base.keys(fp::FinancialPortfolio) = collect(keys(fp.positions))

Base.length(fp::FinancialPortfolio) = length(fp.positions)

Base.copy(fp::FinancialPortfolio) = FinancialPortfolio(copy(fp.positions))

function Base.similar(fp::FinancialPortfolio)
    ew = 1 / length(fp.positions)
    # initialize an equal weight portfolio with same keys
    fp2 = copy(fp)
    for k in eachindex(fp2.positions)
        fp2.positions[k] = ew
    end
    return fp2
end

function Base.:(==)(fp1::FinancialPortfolio{T1},fp2::FinancialPortfolio{T2}) where {T1,T2}
    # short circuit if these are the exact same object
    fp1.positions === fp2.positions && return true

    # otherwise check for equality of types, keys, values
    typematch = T1 == T2
    !typematch && return false

    namematch = sort(keys(fp1)) == sort(keys(fp2))
    !namematch && return false

    for k in eachindex(fp1.positions)
        vm_i = fp1.positions[k] == fp2.positions[k]
        if !vm_i
            return false
        end
    end

    return true
end

end # module


text:
@testset "Periods" begin
    s = 1.0seconds
    m = 1.0minutes
    h = 1.0hours
    d = 1.0days
    y = 1.0years
    c = 1.0centuries
    @test s == Period{seconds}(1.0)
    @test m == Period{minutes}(1.0)
    @test h == Period{hours}(1.0)
    @test d == Period{days}(1.0)
    @test y == Period{years}(1.0)
    @test c == Period{centuries}(1.0)

    @test seconds(s) == 1.0seconds
    @test seconds(m) == 60.0seconds
    @test seconds(h) == 3600.0seconds
    @test seconds(d) == 86400.0seconds
    @test seconds(y) == 3.15576e7seconds
    @test seconds(c) == 3.15576e9seconds

    @test minutes(s) == (1.0 / 60.0)minutes
    @test minutes(m) == 1.0minutes
    @test minutes(h) == 60.0minutes
    @test minutes(d) == 1440.0minutes
    @test minutes(y) == 525960.0minutes
    @test minutes(c) == 5.2596e7minutes

    @test hours(s) == (1.0 / 3600.0)hours
    @test hours(m) == (1.0 / 60.0)hours
    @test hours(h) == 1.0hours
    @test hours(d) == 24.0hours
    @test hours(y) == 8766.0hours
    @test hours(c) == 876600.0hours

    @test days(s) == (1.0 / 86400.0)days
    @test days(m) == (1.0 / 1440.0)days
    @test days(h) == (1.0 / 24.0)days
    @test days(d) == 1.0days
    @test days(y) == 365.25days
    @test days(c) == 36525.0days

    @test years(s) == (1.0 / 3.15576e7)years
    @test years(m) == (1.0 / 525960.0)years
    @test years(h) == (1.0 / 8766.0)years
    @test years(d) == (1.0 / 365.25)years
    @test years(y) == 1.0years
    @test years(c) == 100.0years

    @test centuries(s) == (1.0 / 3.15576e9)centuries
    @test centuries(m) == (1.0 / 5.2596e7)centuries
    @test centuries(h) == (1.0 / 876600.0)centuries
    @test centuries(d) == (1.0 / 36525.0)centuries
    @test centuries(y) == (1.0 / 100.0)centuries
    @test centuries(c) == 1.0centuries

    @test zero(Period{years}) == 0.0years
    @test zero(1years) == 0years
    @test zero(1.0years) == 0.0years

    @test unit(1years) == years
    @test collect(1seconds:3seconds) == [1seconds, 2seconds, 3seconds]
    @test collect(1.0seconds:3.0seconds) == [1.0seconds, 2.0seconds, 3.0seconds]
end

ext: jl


text:
#!julia
include("../src/Tsp.jl")
using Tsp.Model
using Tsp.Solvers
using Tsp.Debug
using Base.Dates

function main(START)
    cities, seconds = load_cities(STDIN)
    timelimit = seconds > 0 ? timeLimit = Second(seconds) : seconds

    # @log @printf("timelimit: %s\n", timeLimit)
    # println("loaded $(length(cities)) cities")
    # tic()
    grid, geometry = make_grid(cities)
    @info begin
        @printf("Grid: %s, %s (%g x %g u; %g u^2) \n", geometry.sw, geometry.ne, geometry.width, geometry.height, geometry.area)
        @printf("Density: %g\n", length(cities) / geometry.area)
        print_chunk_density(grid)
    end
    # toc()
    # x = Tsp.Debug.@config(Debug.DEBUG)
    # @log @printf("loaded grid\n")
    local route::Route
    # method = Solvers.genetic
    # solve = Solvers.solver(method)
    # if method == Solvers.tabu
    #     rounds, restart_after, quit_after, min_tabu_duration, max_tabu_duration, n_neighbours, jump_radius = tabu_params(length(cities), length(grid.chunks))

    #     route = solve(grid,
    #                         rounds=rounds,
    #                         restart_after=restart_after,
    #                         quit_after=quit_after,
    #                         min_tabu_duration=min_tabu_duration,
    #                         max_tabu_duration=max_tabu_duration,
    #                         n_neighbours=n_neighbours,
    #                         jump_radius=jump_radius)

    # elseif method == Solvers.annealing
    #     # tic()
    #     s0=nearest_neighbour(grid.cities)
    #     @info open("initial", "w") do f
    #         println(f, "x,y,id")
    #         for c in s0.cities
    #             println(f, "$(c.x),$(c.y),$(c.id)")
    #         end
    #     end
    #     route = solve(
    #         grid=grid,
    #         s0=s0,
    #         # s0=Route([grid.cities[1], shuffle(grid.cities[2:end-1])..., grid.cities[1]]),
    #         T0=100.0,
    #         limit=10000000,
    #         starttime=START,
    #         timelimit=timelimit)
    #     println("initial solution: ", s0.score)
    # end
    # println("final solution: ", route.score)

    # @info open("solution", "w") do f
    #     println(f, "x,y,id")
    #     for c in route.cities
    #         println(f, "$(c.x),$(c.y),$(c.id)")
    #     end
    # end
    route = Tsp.Solvers.Genetic.find_optimal(
        grid, START, timelimit,
        pop_size=100,
        max_generations=500,
        init_mutation_p=.05,
        stall_limit=100,
        use_nn=true
    )
    println(STDOUT, route.score)
    for city in route.cities
            println(STDERR, city.id)
    end
    # @log println(Dates.CompoundPeriod(now()-START))

end

main(now())

function tabu_params(n_cities, n_chunks)
    quit_after = 3
    rounds = 100
    if n_cities < 100
        n_neighbours = Int(.5 * (n_cities^2))
        restart_after = 75
    elseif n_cities < 1000
        n_neighbours = n_cities * 10
        restart_after = 50
    elseif n_cities < 10000
        rounds = 500
        n_neighbours = max(2500, round(Int, .5 * n_cities))
        restart_after = 45
        quit_after = 3
    else
        rounds = 200
        n_neighbours = 3000
        restart_after = 20
        quit_after = 3
    end

    min_tabu_duration = div(restart_after, 5)
    rounds = restart_after * quit_after * 2
    max_tabu_duration = min(ceil(Int, min_tabu_duration * 1.5), round(Int(rounds/quit_after)))
    jump_radius = round(Int, n_cities/(n_chunks * 2)) # half of avg chunk density
    return (rounds, restart_after, quit_after, min_tabu_duration, max_tabu_duration, n_neighbours, jump_radius)
end

function annealing_params(n_cities, density)
    T0=100.0
    return T0
end


text:
## evaluate binary operations of symbolic objects
## XXX -- this may prove too narryw with the use of ::Sym
+{T<:SymbolicObject}(x::T, y::SymbolicObject) = pycall(sympy[:Add], Sym, x, y)::T
*{T<:SymbolicObject, S<:SymbolicObject}(x::T, y::S) = pycall(sympy[:Mul], Sym, x, y)::T
-{T<:SymbolicObject, S<:SymbolicObject}(x::T, y::S) = x + (-y)
-(x::SymbolicObject)                    =  (-1) * x
/(x::Sym, y::Sym) = x * pycall(sympy[:Pow], Sym, y, -1)::Sym
^(x::SymbolicObject, y::Rational) = x^convert(Sym,y)
^(x::SymbolicObject, y::Integer) = x^convert(Sym,y) # no Union{Integer, Rational}, as that has ambiguity
^(x::Sym, y::Sym) = pycall(sympy[:Pow], Sym, x, y)::Sym
//(x::SymbolicObject, y::Int) = x / Sym(y)
//(x::SymbolicObject, y::Rational) = x / Sym(y)
//(x::SymbolicObject, y::SymbolicObject) = x / y

\(x::SymbolicObject, y::SymbolicObject) = (y'/x')' # ?

inv(x::SymbolicObject) = x\one(x)


text:
################################################################################
## Loading packages that are used.                                            ##
################################################################################
using CSV, DataFrames, DataFramesMeta
using Gadfly
using ZipFile
using FTPClient
#################################################################################

#################################################################################
## Loading data from CDPR FTP server.                                          ##
#################################################################################
"""
    download_info(year1[, year2])

Saves the zip files from CDPR server into your computer memory
"""
function download_files(year1::Int, year2=Union{Nothing,Int})
    ftp_init()
    options = RequestOptions(hostname="transfer.cdpr.ca.gov/pub/outgoing/pur_archives/")
    global files = DataFrame(year = [], files = [])

    if year2 != nothing
        first_year = minimum([year1, year2])
        second_year = maximum([year1, year2])
        years = collect(UnitRange(first_year, second_year))
    else
        years = [year1]
    end

    for i in 1:length(years)
        name = "pur" * string(years[i]) *".zip"
        resp = (year = years[i], files = ftp_get(options, name))
        push!(files, resp)
    end
    ftp_cleanup()
    return files
end
#################################################################################

#################################################################################
## Unzipping and loading files into memory                                      ##
#################################################################################
"""
    get_files(year, kind)

Unizips and loads the desires files from the zip file from the year and save it in
a Dataframe.

kind = :chemicals loads the file "chemical.txt", which contains a list with
all approved chemicals for each year and their chem_codes.
kind = :sites loads the file "site.txt", which contains a list with all sites
(commodities) that can be reported in CA for that year and their site_codes.
kind = :counties loads the the file "county.txt", which contains a list with
all counties in CA and their county_cd (codes).
kind = :products loads the files "product.txt", which contains a list with all approved
chemical products in CA for any year. The resulting dataframe contains
"""
function get_files(year::Int, kind::Symbol)
    for i in 1:size(files)[1]
        if files.year[i] == 2016
            reading = ZipFile.Reader(files.files[i].body)
            for i in 1:length(reading.files)
                if kind == :chemicals && reading.files[i].name == "pur2016/chemical.txt"
                    chemicals = CSV.File(IOBuffer(read(reading.files[i])))
                    return chemicals
                    break
                end
                if kind == :sites && reading.files[i].name == "pur2016/site.txt"
                    sites = CSV.File(IOBuffer(read(reading.files[i])))
                    return sites
                    break
                end
                if kind == :counties && reading.files[i].name == "pur2016/county.txt"
                    counties = CSV.File(IOBuffer(read(reading.files[i])))
                    return counties
                    break
                end
                if kind == :products && reading.files[i].name == "pur2016/product.txt"
                    products = CSV.File(IOBuffer(read(reading.files[i])))
                    return products
                    break
                end
            end
        end

        if files.year[i] == year
            reading = ZipFile.Reader(files.files[i].body)
            for i in 1:length(reading.files)
                if kind == :chemicals && reading.files[i].name == "chemical.txt"
                    chemicals = CSV.File(IOBuffer(read(reading.files[i])))
                    return chemicals
                    break
                end
                if kind == :sites && reading.files[i].name == "site.txt"
                    sites = CSV.File(IOBuffer(read(reading.files[i])))
                    return sites
                    break
                end
                if kind == :counties && reading.files[i].name == "county.txt"
                    counties = CSV.File(IOBuffer(read(reading.files[i])))
                    return counties
                    break
                end
                if kind == :products && reading.files[i].name == "product.txt"
                    products = CSV.File(IOBuffer(read(reading.files[i])))
                    return products
                    break
                end
            end
        end
    end
end

# Function to extract the differences between the approved chemicals between years
function compare_chem_years(year1::Int, year2::Int)
    join(get_files(maximum([year1, year2]), :chemicals),
         get_files(minimum([year1, year2]), :chemicals),
         on =:chem_code, kind = :anti)
end

# Function to extract the differences between the approved sites between years.
function compare_site_years(year1::Int, year2::Int)
    join(get_files(maximum([year1, year2]), :sites),
         get_files(minimum([year1, year2]), :sites),
         on =:site_code, kind = :anti)
end

function compare_product_years(year1::Int, year2::Int)
    join(get_files(maximum([year1, year2]), :products),
         get_files(minimum([year1, year2]), :products),
         on =:site_code, kind = :anti)
end

# find sites using any word (site) as a string.
function find_product_code(product::AbstractString, year::Int)
    products = get_files(year, :products)
    filtered = DataFrame(prod=[], site_name=[])
    for i in 1:size(sites,1)
        matching = match(Regex(uppercase(site)), sites[i,2])
        if matching != nothing
            push!(filtered, (site_code = sites[i,1], site_name=sites[i,2]))
        end
    end
    return filtered
end

# find sites using any word (site) as a string.
function find_site_code(site::AbstractString, year::Int)
    sites = get_files(year, :sites)
    filtered = DataFrame(site_code=[], site_name=[])
    for i in 1:size(sites,1)
        matching = match(Regex(uppercase(site)), sites[i,2])
        if matching != nothing
            push!(filtered, (site_code = sites[i,1], site_name=sites[i,2]))
        end
    end
    return filtered
end

# Extract the county code
function find_county_code(county_name::AbstractString, year::Int)
    counties = get_files(year, :counties)
    county_code = []
    for i in 1:length(counties[:couty_name])
        if counties[i,2] == uppercase(county_name)
            county_code = counties[i,1]
            break
        end
    end
    return county_code
end

################################################################################
## Dealing with the PUR.txt file                                              ##
################################################################################
# Function to load pur file from disk, using the county name as a string and the year
function load_pur(county_name::AbstractString, year::Int, dropmissing::Bool = true)
    county_code = find_county_code(county_name, year)

    string_year = string(year)

    pur = ()
    for i in 1:length(files.files)
        if files.year[i] == 2016
            reading = ZipFile.Reader(files.files[i].body)
            for i in 1:length(reading.files)
                if reading.files[i].name == "pur2016/udc" * string_year[3:4] * "_"* string(county_code) * ".txt"
                    pur = CSV.File(IOBuffer(read(reading.files[i])), escapechar='#',
                    types=Dict(18 => Union{Int64, Char, Missing},20 => Union{Int64, Char, Missing}))
                end
            end
        end
        if files.year[i] == year
            reading = ZipFile.Reader(files.files[i].body)
            for i in 1:length(reading.files)
                if reading.files[i].name == "udc" * string_year[3:4] * "_"* string(county_code) * ".txt"
                    pur = CSV.File(IOBuffer(read(reading.files[i])), escapechar='#',
                    types=Dict(18 => Union{Int64, Char, Missing},20 => Union{Int64, Char, Missing}))
                end
            end
        end
    end

    pur_df = DataFrame(pur)

    pur_df = @select(pur_df, :use_no, :prodno, :chem_code, :prodchem_pct, :lbs_chm_used,
    :lbs_prd_used, :amt_prd_used, :unit_of_meas, :acre_planted, :unit_planted,
    :acre_treated, :unit_treated, :applic_cnt, :applic_dt, :grower_id, :site_code)

    # add column for year
    pur_df[:year] = repeat([year], outer = size(pur_df, 1))

    # add chem name
    pur_df[:chem_name] = repeat(["0"],outer = size(pur_df, 1))

    if dropmissing == true
        pur_df = dropmissing!(pur_df)
    end

    chem_codes = get_files(year, :chemicals)
    chemcodes = Dict()
    for i in 1:size(chem_codes,1)
        push!(chemcodes, (chem_codes[i,1] => chem_codes[i,3]))
    end

    for i in 1:size(pur_df,1)
        chemical = pur_df[i,:chem_code]
        if typeof(chemical) != Missing
            pur_df[i,:chem_name] = chemcodes[chemical]
        end
    end
    return pur_df
end

# Select a site code using the find_site_code function
function select_site_code(df::DataFrame, site_code::Int)
    @where(df, :site_code .== site_code)
end

function select_insecticide(df::DataFrame, insecticide_name::AbstractString)
    insecticide_name = string(uppercase(insecticide_name))
    @where(df, :chem_name .== insecticide_name)
end

# Function to append the data from different years
function append_years(county_name::AbstractString, first_year::Int, last_year::Union{Nothing,Int},
                      dropmissing::Bool = true)
    if last_year==nothing
        total = load_pur(county_name, first_year, false)
        return total
    end

    if dropmissing == false
        total = load_pur(county_name, first_year, false)
        for i in first_year+1:last_year
            pur = load_pur(county_name, i, false)
            total= append!(total, pur)
        end
        return total
    end

    total = load_pur(county_name, first_year)
    for i in first_year+1:last_year
        pur = load_pur(county_name, i)
        total= append!(total, pur)
    end
    return total
end

# Function to obtain the information in insecticide use from a given commodity and county
function insecticide_use(county_name::AbstractString,
                         insecticide_name::AbstractString,
                         site_code::Int,
                         first_year::Int, last_year::Union{Nothing,Int},
                         dropmissing::Bool=true)
    total = append_years(county_name, first_year, last_year, dropmissing)
    total = select_site_code(total, site_code)
    total = select_insecticide(total, insecticide_name)
    total = by(total, [:year]) do df
    DataFrame(lbs_chm_used = sum(df[:lbs_chm_used]),
              N = size(df,1))
    end
    return total
end

function plot_insecticide_use(ounty_name::AbstractString,
                         insecticide_name::AbstractString,
                         site_code::Int,
                         first_year::Int, last_year::Union{Nothing,Int},
                         dropmissing::Bool=true)
    result = insecticide_use(county_name, insecticide_name, site_code, first_year, last_year,
                             dropmissing)
    plot(result, x = :year, y = :lbs_chm_used, Geom.line)
end


text:
#=

  Mr Smith problem in Julia ConstraintSolver.jl

  From an IF Prolog example (http://www.ifcomputer.de/)
  """
  The Smith family and their three children want to pay a visit but they
  do not all have the time to do so. Following are few hints who will go
  and who will not:
      o If Mr Smith comes, his wife will come too.
      o At least one of their two sons Matt and John will come.
      o Either Mrs Smith or Tim will come, but not both.
      o Either Tim and John will come, or neither will come.
      o If Matt comes, then John and his father will
        also come.
  """

  The answer should be:
    Mr_Smith_comes      =  0
    Mrs_Smith_comes     =  0
    Matt_comes          =  0
    John_comes          =  1
    Tim_comes           =  1

  Model created by Hakan Kjellerstrand, hakank@gmail.com
  See also my Julia page: http://www.hakank.org/julia/

=#

using ConstraintSolver, JuMP
using Cbc, GLPK, Ipopt
const CS = ConstraintSolver
include("constraints_utils.jl")

function mr_smith(print_solutions=true,all_solutions=true)

    cbc_optimizer = optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 0)
    glpk_optimizer = optimizer_with_attributes(GLPK.Optimizer)
    ipopt_optimizer = optimizer_with_attributes(Ipopt.Optimizer)

    model = Model(optimizer_with_attributes(CS.Optimizer,   "all_solutions"=> all_solutions,
                                                            # "all_optimal_solutions"=>all_solutions,
                                                            "logging"=>[],

                                                            "traverse_strategy"=>:BFS,
                                                            # "traverse_strategy"=>:DFS,
                                                            # "traverse_strategy"=>:DBFS,

                                                            # "branch_split"=>:Smallest,
                                                            # "branch_split"=>:Biggest,
                                                            "branch_split"=>:InHalf,

                                                            # https://wikunia.github.io/ConstraintSolver.jl/stable/options/#branch_strategy-(:Auto)
                                                            "branch_strategy" => :IMPS, # default
                                                            # "branch_strategy" => :ABS, # Activity Based Search
                                                            # "activity.decay" => 0.999, # default 0.999
                                                            # "activity.max_probes" => 10, # default, 10
                                                            # "activity.max_confidence_deviation" => 20, # default 20

                                                            # "simplify"=>false,
                                                            # "simplify"=>true, # default

                                                            "time_limit"=>6,

                                                            # "backtrack" => false, # default true
                                                            # "backtrack_sorting" => false, # default true

                                                            # "lp_optimizer" => cbc_optimizer,
                                                            # "lp_optimizer" => glpk_optimizer,
                                                            # "lp_optimizer" => ipopt_optimizer,
                                        ))

    n = 5
    @variable(model, x[1:n], Bin)
    mr_smith,mrs_smith,matt,john,tim = x
    people = ["mr_smith","mrs_smith","matt","john","tim"]

    # If Mr Smith comes, his wife will come too.
    @constraint(model,mr_smith => {mrs_smith == 1})

    # At least one of their two sons Matt and John will come.
    @constraint(model,matt + john >= 1)

    # Either Mrs Smith or Tim will come, but not both.
    # Mrs_Smith + Tim #= 1,
    @constraint(model,mrs_smith + tim == 1)

    # Either Tim and John will come, or neither will come.
    @constraint(model,tim == john)

    # If Matt comes, then John and his father will also come.
    @variable(model, john_and_mr_smith, Bin)
    @constraint(model, john_and_mr_smith := {john + mr_smith == 2})
    @constraint(model,matt => {john_and_mr_smith == 1})

    # Solve the problem
    optimize!(model)

    status = JuMP.termination_status(model)
    # println("status:$status")
    if status == MOI.OPTIMAL
        num_sols = MOI.get(model, MOI.ResultCount())
        println("num_sols:$num_sols\n")
        if print_solutions
            for sol in 1:num_sols
                # println("solution #$sol")
                x_val = convert.(Integer,JuMP.value.(x; result=sol))
                println("x:$x_val")
                println("These will come:", [people[i] for i in 1:n if x_val[i] == 1])
            end
        end
    else
        println("status:$status")
    end

    return status
end

@time mr_smith()

max_stars_repo_name: Wikunia/hakank

max_issues_repo_name: Wikunia/hakank

max_forks_repo_name: Wikunia/hakank


text:
module GPARs

using AbstractGPs
using LinearAlgebra
using Random

include("gpar.jl")

export GPAR

end


text:
using LinearAlgebra
using SparseArrays
using Test
import Random

using JuMP
using GLPK
using LightGraphs
using GraphPlot
import SCS
using Colors: @colorant_str

# TODO explain problem

function compute_cut_value(g::LightGraphs.AbstractGraph, w::AbstractMatrix, vertex_subset)
    return sum(w[i,j] for i in vertices(g) for j in neighbors(g, i) if in(i, vertex_subset)  != in(j, vertex_subset))
end

g = SimpleGraph(6)
add_edge!(g, 1, 2)
add_edge!(g, 1, 5)
add_edge!(g, 2, 3)
add_edge!(g, 2, 5)
add_edge!(g, 3, 4)
add_edge!(g, 4, 5)
add_edge!(g, 4, 6);

GraphPlot.gplot(g, nodelabel=1:6)

w = spzeros(6, 6)

for e in edges(g)
   (i, j) = Tuple(e)
   w[i,j] = 1
end

@test compute_cut_value(g, w, (2, 5, 6)) ≈ 5
@show compute_cut_value(g, w, (2, 5, 6));

n = nv(g)

linear_max_cut = Model(GLPK.Optimizer)

@variable(linear_max_cut, z[1:n,1:n], Bin)
@variable(linear_max_cut, x[1:n], Bin)
@constraint(linear_max_cut, [i = 1:n, j = 1:n; (i, j) in edges(g)], z[i,j] <= x[i] + x[j])
@constraint(linear_max_cut, [i = 1:n, j = 1:n; (i, j) in edges(g)], z[i,j] <= 2 - (x[i] + x[j]))

@objective(linear_max_cut, Max, dot(w, z))

optimize!(linear_max_cut)

@test objective_value(linear_max_cut) ≈ 6.0
@show objective_value(linear_max_cut);

x_linear = value.(x)
@show x_linear;
@show [(i,j) for i in 1:n-1 for j in i+1:n if JuMP.value.(z)[i,j] > 0.5];

nodecolor = [colorant"lightseagreen", colorant"orange"]
all_node_colors = [nodecolor[round(Int, xi + 1)]  for xi in x_linear]
GraphPlot.gplot(g, nodelabel=1:6, nodefillc=all_node_colors)

sdp_max_cut = Model(optimizer_with_attributes(SCS.Optimizer, "verbose" => 0))

@variable(sdp_max_cut, Y[1:n,1:n] in PSDCone())
@constraint(sdp_max_cut, [i = 1:n], Y[i,i] == 1)
@objective(sdp_max_cut, Max, 1/4 * sum(w[i,j] * (1 - Y[i,j]) for i in 1:n for j in 1:n))

optimize!(sdp_max_cut)

@test objective_value(linear_max_cut) >= 6.0
@show objective_value(linear_max_cut);

F = svd(value.(Y))
U = F.U * Diagonal(sqrt.(F.S))

Random.seed!(33)

x = randn(size(U, 2))
xhat = sign.(U * x) .> 0

@show collect(zip(round.(Int, x_linear), xhat))
@test LinearAlgebra.norm1(round.(Int, x_linear) .- xhat) == 0 || LinearAlgebra.norm1(round.(Int, x_linear) .- xhat) == n

nodecolor = [colorant"lightseagreen", colorant"orange"]
all_node_colors = [nodecolor[xi + 1] for xi in xhat]
GraphPlot.gplot(g, nodelabel=1:6, nodefillc=all_node_colors)


text:
function test(name)
    print("Since you passed the setup, we can say 'Hello, $(name)'")
end

function greet()
    print("I am a 2nd function")
end


text:
include((@__DIR__)*"/../preamble.jl")

include("reflecting_model/model_def.jl")
include("default_params.jl")

tvec = CSV.read((@__DIR__)*"/hitting_times/data/point_mass/order_1_model_dg1.csv",DataFrame).t

pth = mkpath((@__DIR__)*"/hitting_times/data/errors")
for ic_string in ["point_mass","exp"]
    ks_errors = DataFrame(DG_limit=[])
    l1_errors = DataFrame(DG_limit=[])
    ks_error_row = zeros(1)
    l1_error_row = zeros(1)

    sim_cdf = CSV.read("hitting_times_model/hitting_times/data/"*ic_string*"/sims_cdf_evaluated.csv",DataFrame)
    for order in 1:2:21
        for (c,model_string) in enumerate(["dg1"])
            approx_data = CSV.read((@__DIR__)*"/hitting_times/data/"*
                ic_string*"/order_"*string(order)*"_model_"*model_string*".csv",DataFrame)
            ks_error = max(
                maximum(abs.(approx_data.phase_1-sim_cdf.phase_1)),
                maximum(abs.(approx_data.phase_2-sim_cdf.phase_2)),
            )
            l1_cdf = sum(abs.(approx_data.phase_1[1:end-1]+approx_data.phase_1[2:end]-
                    sim_cdf.phase_1[1:end-1]-sim_cdf.phase_1[2:end]).*diff(approx_data.t))/2 +
                sum(abs.(approx_data.phase_2[1:end-1]+approx_data.phase_2[2:end]-
                    sim_cdf.phase_2[1:end-1]-sim_cdf.phase_2[2:end]).*diff(approx_data.t))/2

            ks_error_row[c] = ks_error
            l1_error_row[c] = l1_cdf
        end
        push!(ks_errors,ks_error_row)
        push!(l1_errors,l1_error_row)
    end
    CSV.write(pth*"/ks_"*ic_string*".csv",ks_errors)
    CSV.write(pth*"/l1_"*ic_string*".csv",l1_errors)
end


text:
d = fresh(pokerlayout())
println("A new poker deck:")
println(pretty(d))

shuffle!(d)
println()
println("The deck shuffled:")
println(pretty(d))

n = 4
println()
println("Deal ", n, " hands:")
for i in 1:n
    h = deal!(d)
    println(pretty(h))
end

println()
println("And now the deck contains:")
println(pretty(d))

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
@shorthands meanplot
@shorthands autocorplot
@shorthands mixeddensity
@shorthands traceplot

struct _TracePlot; c; val; end
struct _MeanPlot; c; val;  end
struct _DensityPlot; c; val;  end
struct _HistogramPlot; c; val;  end
struct _AutocorPlot; lags; val;  end

# define alias functions for old syntax
const translationdict = Dict(
                        :traceplot => _TracePlot,
                        :meanplot => _MeanPlot,
                        :density => _DensityPlot,
                        :histogram => _HistogramPlot,
                        :autocorplot => _AutocorPlot
                      )

const supportedplots = push!(collect(keys(translationdict)), :mixeddensity)

@recipe function f(c::AbstractChains, i::Int; colordim = :chain, barbounds = (0, Inf), maxlag = nothing)
    st = get(plotattributes, :seriestype, :traceplot)

    if colordim == :parameter
        title := "Chain $(c.chains[i])"
        labels := c.names
        val = c.value[:, :, i]
    elseif colordim == :chain
        title := c.names[i]
        labels := map(k -> "Chain $(c.chains[k])", 1:size(c)[3])
        val = c.value[:, i, :]
    else
        throw(ArgumentError("`colordim` must be one of `:chain` or `:parameter`"))
    end

    if st == :mixeddensity
        discrete = MCMCChain.indiscretesupport(c, barbounds)
        st = if colordim == :chain
            discrete[i] ? :histogram : :density
        else
            # NOTE: It might make sense to overlay histograms and density plots here.
            :density
        end
        seriestype := st
    end

    if st == :autocorplot
        lags = 0:(maxlag === nothing ? round(Int, 10 * log10(length(c.range))) : maxlag)
        ac = MCMCChain.autocor(c, lags=collect(lags))
        val = colordim == :parameter ? ac.value[:, :, i]' : ac.value[i, :, :]
        _AutocorPlot(lags, val)
    elseif st ∈ supportedplots
        translationdict[st](c, val)
    else
        c.range, val
    end
end

@recipe function f(p::_DensityPlot)
    xaxis --> "Sample value"
    yaxis --> "Density"
    [collect(skipmissing(p.val[:,k])) for k in 1:size(p.val, 2)]
end

@recipe function f(p::_HistogramPlot)
    xaxis --> "Sample value"
    yaxis --> "Frequency"
    fillalpha --> 0.7
    [collect(skipmissing(p.val[:,k])) for k in 1:size(p.val, 2)]
end

@recipe function f(p::_MeanPlot)
    seriestype := :path
    xaxis --> "Iteration"
    yaxis --> "Mean"
    p.c.range, MCMCChain.cummean(p.val)
end

@recipe function f(p::_AutocorPlot)
    seriestype := :path
    xaxis --> "Lag"
    yaxis --> "Autocorrelation"
    p.lags, p.val
end

@recipe function f(p::_TracePlot)
    seriestype := :path
    xaxis --> "Iteration"
    yaxis --> "Sample value"
    p.c.range, p.val
end

@recipe function f(c::MCMCChain.AbstractChains;
                   width = 500,
                   height = 250,
                   colordim = :chain
                  )

    ptypes = get(plotattributes, :seriestype, (:traceplot, :mixeddensity))

    # sanity check
    ptypes = ptypes isa AbstractVector || ptypes isa Tuple ? ptypes : (ptypes,)
    @assert all(map(ptype -> ptype ∈ supportedplots, ptypes))

    nrows, nvars, nchains = size(c)
    ntypes = length(ptypes)
    N = colordim == :chain ? nvars : nchains
    layout := (N, ntypes)
    size --> (ntypes*width, N*height)
    indices = reshape(1:N*ntypes, ntypes, N)'

    legend --> false

    for (j, ptype) in enumerate(ptypes)
        for i in 1:N
            @series begin
                subplot := indices[i, j]
                colordim := colordim
                seriestype := ptype
                c, i
            end
        end
    end
end

function plot(c::AbstractChains, psyms::Vector{Symbol}; args...)
    @assert all(map(psym -> haskey(translationdict, psym), psyms))
    @warn "This syntax is deprecated, please use plot(c; seriestype = $psyms) instead."
    return plot(c; seriestype = map(psym -> translationdict[psym], psyms))
end

function plot(c::AbstractChains, psym::Symbol; args...)
    @assert haskey(translationdict, psym)
    @warn "This syntax is deprecated, please use plot(c, seriestype = $(translationdict[psym])) instead"
    return plot(c; seriestype = psym, args...)
end


text:
import Base: iterate
import Base: length

struct BitPackedIterator
	arr::Vector{UInt8}
	bitwidth::UInt8
end

function iterate(bpiter::BitPackedIterator, pos=0)
	if pos >= sizeof(UInt)*length(bpiter.arr)
		return nothing
	end

	OnesMask = UInt8(2^bpiter.bitwidth - 1)

	bytes_to_skip, bits_to_shift = divrem(pos, sizeof(UInt))
	value = (bpiter.arr[1 + bytes_to_skip] >> bits_to_shift) & OnesMask

	# say pos=6 and bitwidth is 3 then after shifting
	# 6 spots, there are only two bits left and so we need one more bit
	# from the next byte
	bits_left_in_byte_post_shift = sizeof(UInt) - bits_to_shift

	if bits_left_in_byte_post_shift < bpiter.bitwidth
		# in order to look into the next byte you need to have the next byte
		if (1 + bytes_to_skip + 1) <= length(bpiter.arr)
			n_bits_needed_from_next_byte = bpiter.bitwidth - bits_left_in_byte_post_shift

			mask = 2^n_bits_needed_from_next_byte-1
			bits_from_next_byte = bpiter.arr[1 + bytes_to_skip + 1] & mask

			value = value | (bits_from_next_byte << n_bits_needed_from_next_byte)
		end
	end

	new_pos = pos + bpiter.bitwidth
	return value, new_pos
end

IteratorSize(BitPackedIterator) = HasLength()

length(bpiter::BitPackedIterator) = ceil(sizeof(UInt8)*length(bpiter.arr)/ bpiter.bitwidth)


text:
using OrdinaryDiffEq
using DiffEqCallbacks
using QuasiMonteCarlo
using StatsBase
using Plots

function sir_ode!(du,u,p,t)
    (S,I,R) = u
    (β,c,γ) = p
    N = S+I+R
    @inbounds begin
        du[1] = -β*c*I/N*S
        du[2] = β*c*I/N*S - γ*I
        du[3] = γ*I
    end
    nothing
end;

tmax = 10000.0
tspan = (0.0,tmax)

cb_ss = TerminateSteadyState();

N = 1000.0;
u0 = [990.0,10.0,0.0];
p = [0.05,10.0,0.25]; # β,c,γ

n_samples = 1000 # Number of samples
# Parameters are β, c, γ, I₀
lb = [0.01, 5.0, 0.1, 1.0]
ub = [0.1, 20.0, 1.0, 50.0];

pu0 = QuasiMonteCarlo.sample(n_samples,lb,ub,LatinHypercubeSample());

prob_ode = ODEProblem(sir_ode!,u0,tspan,p);

f1 = function(pu0)
  p = pu0[1:3]
  I0 = pu0[4]
  u0 = [N-I0,I0,0.0]
  prob = remake(prob_ode;p=p,u=u0)
  sol = solve(prob, ROS34PW3(),callback=cb_ss)
  [maximum(sol[2,:]), sol.t[argmax(sol[2,:])], sol[end][3]]
end;

m_serial = [f1(pu0[:,i]) for i in 1:n_samples]
m_serial = hcat(m_serial...); # convert into matrix

pf1 = function (pu0)
  p = pu0[1:3,:]
  I0 = pu0[4,:]
  prob_func(prob,i,repeat) = remake(prob;p=p[:,i],u=[N-I0[i],I0[i],0.0])
  ensemble_prob = EnsembleProblem(prob_ode,prob_func=prob_func)
  sol = solve(ensemble_prob,ROS34PW3(),EnsembleThreads();trajectories=size(p,2))
  out = zeros(3,size(p,2))
  for i in 1:size(p,2)
    out[1,i] = maximum(sol[i][2,:])
    out[2,i] = sol[i].t[argmax(sol[i][2,:])]
    out[3,i] = sol[i][end][3]
  end
  out
end;

m_parallel = pf1(pu0);

l = @layout [a b; c d]
pl1 = scatter(pu0[1,:],m_parallel[1,:],title="Peak infected",xlabel="β",ylabel="Number")
pl2 = scatter(pu0[2,:],m_parallel[1,:],title="Peak infected",xlabel="c",ylabel="Number")
pl3 = scatter(pu0[3,:],m_parallel[1,:],title="Peak infected",xlabel="γ",ylabel="Number")
pl4 = scatter(pu0[4,:],m_parallel[1,:],title="Peak infected",xlabel="I₀",ylabel="Number")
plot(pl1,pl2,pl3,pl4,layout=l,legend=false)

l = @layout [a b; c d]
pl1 = scatter(pu0[1,:],m_parallel[2,:],title="Peak time",xlabel="β",ylabel="Time")
pl2 = scatter(pu0[2,:],m_parallel[2,:],title="Peak time",xlabel="c",ylabel="Time")
pl3 = scatter(pu0[3,:],m_parallel[2,:],title="Peak time",xlabel="γ",ylabel="Time")
pl4 = scatter(pu0[4,:],m_parallel[2,:],title="Peak time",xlabel="I₀",ylabel="Time")
plot(pl1,pl2,pl3,pl4,layout=l,legend=false)

l = @layout [a b; c d]
pl1 = scatter(pu0[1,:],m_parallel[3,:],title="Final size",xlabel="β",ylabel="Number")
pl2 = scatter(pu0[2,:],m_parallel[3,:],title="Final size",xlabel="c",ylabel="Number")
pl3 = scatter(pu0[3,:],m_parallel[3,:],title="Final size",xlabel="γ",ylabel="Number")
pl4 = scatter(pu0[4,:],m_parallel[3,:],title="Final size",xlabel="I₀",ylabel="Number")
plot(pl1,pl2,pl3,pl4,layout=l,legend=false)

corkendall(pu0',m_parallel')


text:
# cudnnActivationMode_t
const CUDNN_ACTIVATION_SIGMOID = Cint(0)
const CUDNN_ACTIVATION_RELU = Cint(1)
const CUDNN_ACTIVATION_TANH = Cint(2)
const CUDNN_ACTIVATION_CLIPPED_RELU = Cint(3)
const CUDNN_ACTIVATION_ELU = Cint(4)

"""
    ActivationDesc

* coef: floating point number to specify the clipping threashold when the activation
mode is set to CUDNN_ACTIVATION_CLIPPED_RELU or to specify the alpha coefficient
when the activation mode is set to CUDNN_ACTIVATION_ELU
"""
mutable struct ActivationDesc
    ptr::Cptr

    function ActivationDesc(mode::Cint, coef::Float64)
        ref = Ref{Cptr}()
        @cudnn :cudnnCreateActivationDescriptor (Ptr{Cptr},) ref
        desc = new(ref[])
        finalizer(desc) do x
            @cudnn :cudnnDestroyActivationDescriptor (Cptr,) x.ptr
        end

        @cudnn(:cudnnSetActivationDescriptor,
            (Cptr,Cint,Cint,Cdouble),
            desc, mode, CUDNN_NOT_PROPAGATE_NAN, coef)
        desc
    end
end

const ACTIVATION_DESCS = Dict{Tuple,ActivationDesc}()

Base.unsafe_convert(::Type{Cptr}, desc::ActivationDesc) = desc.ptr

function activation(x::CuArray{T}, mode::Cint, coef::Float64) where T
    actdesc = get!(ACTIVATION_DESCS, (mode,coef)) do
        ActivationDesc(mode, coef)
    end
    xdesc = TensorDesc(x, 4)
    y = similar(x)
    @cudnn(:cudnnActivationForward,
        (Cptr,Cptr,Cptr,Cptr,Cptr,Cptr,Cptr,Cptr),
        gethandle(), actdesc, T[1], xdesc, x, T[0], xdesc, y)
    y
end

sigmoid(x::CuArray) = activation(x, CUDNN_ACTIVATION_SIGMOID, 0.0)
relu(x::CuArray) = activation(x, CUDNN_ACTIVATION_RELU, 0.0)
tanh(x::CuArray) = activation(x, CUDNN_ACTIVATION_TANH, 0.0)

function ∇activation!(y::CuArray{T}, dy, x, dx, mode, coef) where T
    actdesc = ActivationDesc(mode, coef)
    xdesc = TensorDesc(x, 4)
    @cudnn(:cudnnActivationBackward,
        (Cptr,Cptr,Cptr,Cptr,Cptr,Cptr,Cptr,
        Cptr,Cptr,Cptr,Cptr,Cptr),
        gethandle(), actdesc, T[1], xdesc, y, xdesc, dy, xdesc, x, T[1], xdesc, dx)
end

∇sigmoid!(y, dy, x, dx) = ∇activation!(y, dy, x, dx, CUDNN_ACTIVATION_SIGMOID, 0.0)
∇relu!(y, dy, x, dx) = ∇activation!(y, dy, x, dx, CUDNN_ACTIVATION_RELU, 0.0)
∇tanh!(y, dy, x, dx) = ∇activation!(y, dy, x, dx, CUDNN_ACTIVATION_TANH, 0.0)


text:
module Parameters

export bgc_params_default, phyt_params_default
export default_PARF, default_temperature
export update_bgc_params, update_phyt_params

using CUDA
using PlanktonIndividuals.Grids

using PlanktonIndividuals: AbstractMode, CarbonMode, QuotaMode, MacroMolecularMode

include("param_default.jl")
include("param_update.jl")

end


text:
function solve(P::Matrix{T}, q::Vector{T}, A::Matrix{T}, b::Vector{T}, x::Vector{T}; r_min::T=zero(T), r_max::T=T(Inf), max_iter=Inf, kwargs...) where T
    if r_min < 1e-9; r_min = -one(T); end
    if abs(r_min - r_max) < 1e-9
        return solve_boundary(P, q, A, b, (r_min + r_max)/2, x; max_iter=max_iter, kwargs...)
    end

    boundary_data, interior_data = create_data(P, q, A, b, r_min, r_max, x; kwargs...)

    active_boundary = which_boundary(x, r_min, r_max)
    if active_boundary == :upper
        data = boundary_data
        data.r = r_max
    elseif active_boundary == :lower
        data = boundary_data
        data.r = r_min
    else
        data = interior_data
    end

    if data.verbosity > 0
        print_header(data)
        print_info(data)
    end

    while !data.done && data.iteration <= max_iter
        iterate!(data)
        if (isa(data, Data) && data.done)
            active_boundary = which_boundary(data.x, r_min, r_max)
            if (active_boundary == :upper && data.μ < -1e-8) ||
                (active_boundary == :lower && data.μ > 1e-8)
                # Don't terminate if the Lagrange multiplier of the spherical constraint is negative
                data.done = false
            end
        end

        if data.verbosity > 0
            mod(data.iteration, 10*data.printing_interval) == 0 && print_header(data)
            (mod(data.iteration, data.printing_interval) == 0 || data.done) && print_info(data)
        end

        active_boundary = which_boundary(data.x, r_min, r_max)
        if isa(data, GeneralQP.Data) && active_boundary != :interior
            # The QP solver has hit the constraint ‖x‖ = r_min or r_max
            # Thus we switch to the etrs_boundary algorithm
            update_data!(boundary_data, interior_data)
            data = boundary_data
        elseif isa(data, Data)
            if (active_boundary == :upper && data.μ < 0 && data.μ < minimum(data.λ)) ||
                (active_boundary == :lower && data.μ > 0 && -data.μ < minimum(data.λ))
                # The etrs_boundary algorithm is in a "stationary" point
                # with the most negative langrange multiplier being the
                # one for the constraint ‖x‖ = r_max or r_min
                # Thus we switch to the GeneralQP solver
                update_data!(interior_data, boundary_data)
                data = interior_data
            end
        end
    end

    if isa(data, Data)
        return data.x, [data.λ; data.μ]
    else
        return data.x, [data.λ; zero(T)]
    end
end

function create_data(P, q, A, b, r_min, r_max, x_init; kwargs...)
    interior_data = GeneralQP.Data(P, q, A, b, x_init; kwargs..., r_max=r_max, r_min=r_min)
    # interior_data.verbosity = 0
    boundary_data = Data(P, q, A, b, r_max, x_init,
        interior_data.F.QR, interior_data.working_set, interior_data.ignored_set,
        interior_data.A_shuffled, interior_data.b_shuffled; kwargs...)
    return boundary_data, interior_data
end

function update_data!(qp::GeneralQP.Data, etrs::Data)
    qp.iteration = etrs.iteration
    if qp.verbosity > 0
        @info "Changing to QP solver."
    end
    # We only have to recompute qp.F.LDL
    # Adding a (temporarily) constraint parallel to x ensures that the
    # projected hessian is p.d. if a shift from r_max to the interior happens
    add_constraint!(etrs.F, etrs.x/norm(etrs.x))
    qp.x = etrs.x
    qp.F.m = qp.F.QR.n - qp.F.QR.m
    qp.F.artificial_constraints = 0  # No need for them now :)
    GeneralQP.update_views!(qp.F)
    try
        qp.F.U.data .= cholesky!(Symmetric(etrs.F.ZPZ[qp.F.m:-1:1, qp.F.m:-1:1])).U
        qp.F.d .= 1
        # remove the temporal constraint parallel to x
        remove_constraint!(qp.F, qp.F.QR.m)
    catch error
        if isa(error, PosDefException)
            # This can only happen during a shift from r_min to the interior
            if qp.verbosity > 0
                @info "Introducing artifical constraints in QP's working set"
            end
            qp.F = NullspaceHessianLDL(qp.F.P, Matrix(view(qp.A, qp.working_set, :)'))
            etrs.F.QR = qp.F.QR # Link the qr of the qp solver and the constant norm solver
            if qp.F.m == 0 # To many artificial constraints...
                remove_constraint!(qp.F, 0)
            end
        else
            throw(error)
        end
    end
    GeneralQP.update_views!(qp)
    if qp.verbosity > 0
        print_header(qp)
        print_info(qp)
    end
end

function update_data!(etrs::Data, qp::GeneralQP.Data)
    if etrs.verbosity > 0
        @info "Changing to constant-norm QP solver."
    end
    etrs.iteration = qp.iteration
    etrs.x = qp.x
    etrs.r = norm(etrs.x)
    # We only have to recompute etrs.F.ZPZ
    etrs.F.m = etrs.F.QR.n - etrs.F.QR.m
    GeneralQP.update_views!(etrs.F)
    etrs.F.ZPZ .= etrs.F.Z'*etrs.F.P*etrs.F.Z
    etrs.F.ZPZ .= (etrs.F.ZPZ + etrs.F.ZPZ')/2
    if etrs.F.m <= 1
        etrs.y = etrs.F.Z'*etrs.x
        etrs.x0 = etrs.x - etrs.F.Z*etrs.y
        etrs.Zq = etrs.F.Z'*(etrs.q + etrs.F.P*etrs.x0) # Reduced q
        if norm(etrs.y) <= 1e-10
            @warn "LICQ failed. Terminating"
            etrs.done = true
        else
            idx = check_kkt!(etrs)
            !etrs.done && idx <= length(etrs.working_set) && remove_constraint!(etrs, idx)
        end
    end
    GeneralQP.update_views!(etrs)
    if etrs.verbosity > 0
        print_header(etrs)
        print_info(etrs)
    end
end

function iterate!(data::GeneralQP.Data{T}) where T
    GeneralQP.iterate!(data)
end

function print_info(data::GeneralQP.Data{T}) where T
    GeneralQP.print_info(data)
end

function print_header(data::GeneralQP.Data{T}) where T
    GeneralQP.print_header(data)
end

function which_boundary(x, r_min, r_max)
    tol = 1e-8
    if norm(x) - r_max >= -tol
        @assert norm(x) - r_max <= tol "Iterate infeasible for upper norm constraint"
        return :upper
    elseif norm(x) - r_min <= tol
        @assert norm(x) - r_min >= -tol "Iterate infeasible for lower norm constraint"
        return :lower
    else
        return :interior
    end
end


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

f(x) = cos(x)

k = 20
x = iterate(f, 1.0, k) # converges around ≈ 0.73

fp = fixedpoint(x, ϵ = 1e-3) # n = 17


text:
# -*- coding: utf-8 -*-
include("gerda_data_processing/init.jl")
include("gerda_data_processing/fct.jl")
include("fct-sim.jl")

plots_base_path = "pulses/sim/raw_0vbb/plots/"
base_path_raw   = "pulses/sim/raw_0vbb/"
base_path       = "pulses/sim/0vbb/"

parameters = JSON.parsefile("res/parameters.json")
bool_AE = false
std_intersect = 2;

# Search for available channels
channels = []
for c in 0:1:36
    files = glob(base_path_raw*string(c)*"-"*channel_to_name[c]*"/*w_filter-cl*.h5");
    if length(files) > 0
        push!(channels, c)
    end
end

for ch in channels
#     IJulia.clear_output(true)
    Base.run(`clear`)
    str_ch = lpad(ch, 2, "0");
    @info("Processing Ch" * str_ch * " | " * channel_to_name[ch])
    @info(">---------------------------------------------------------------<")
    @info("Read in baselines")
    files  = glob("data_pulses/baselines_filtered/"*str_ch*"-"*channel_to_name[ch]*"/*.h5");

    bl = Table(energy=[], multiplicity=[], timestamp=[], channel=[], waveform=[])

    @showprogress for file in files
        append!(bl, HDF5.h5open(file, "r") do h5f
            LegendHDF5IO.readdata(h5f, "data")
        end)
    end

    # get number of pulses

    files = glob(base_path_raw*string(ch)*"-"*channel_to_name[ch]*"/*w_filter-cl*.h5");

    number_of_pulses = 0;
    key_str = []
    for file in files
        push!(key_str, split(basename(file), "_evts_")[1])
    end
    key_str = unique(key_str)
    for key in key_str
        key_files = glob(base_path_raw*string(ch)*"-"*channel_to_name[ch]*"/"*key*"*.h5")
        last_file = split(basename(key_files[end]), ".h5")[1]
        skipped = 0
        for single_file in key_files
            if length(split(basename(single_file), "_skipped-")) > 1
                skipped += parse(Int64, split(split(basename(single_file), "_skipped-")[end], ".h5")[1])
#                 println(basename(single_file))
#                 return parse(Int64, split(split(split(last_file, "_evts_")[2], "_")[1], "-")[end])
            end
        end
        number_of_pulses += parse(Int64, split(split(split(last_file, "_evts_")[2], "_")[1], "-")[end]) #- skipped
    end
    #
    ## work-around for skipped events.. just increasing the number of indices
    number_of_pulses += 10000
    @info("Find number of simulation pulses: " * string(number_of_pulses))
    @info("Find number of noise baselines: " * string(size(bl,1)))
    # get first_range
    if number_of_pulses <= size(bl,1)
        first_range = [1]
    else
        mult_factor = ceil(number_of_pulses/size(bl,1))
        first = 1
        last = 201
        step = Int(floor((last-first)/(mult_factor)))
        first_range = first:step:last
    end

    @info("This results in "*string( size(bl, 1) *length(first_range) )*" baselines for "*string(number_of_pulses)*" pulses")
    @info(">---------------------------------------------------------------<")
    @info("Indexing the baselines")
    indices = Table(bl=[], first=[])
    for i in 0:1:number_of_pulses
        bl_id = i%size(bl,1) + 1
        first = first_range[Int(ceil((i + 1) / size(bl,1)))]
        append!(indices, Table(bl=[bl_id], first=[first]))
    end

    pulse_id = 1
    @info("Add baselines to pulses")

    prog = Progress(length(files), dt=0.5,
                    barglyphs=BarGlyphs('|','█', ['▁' ,'▂' ,'▃' ,'▄' ,'▅' ,'▆', '▇'],' ','|',),
                    barlen=10)
    for file in files
        create_file = false
        ch_str = lpad(ch, 2, "0");
        filename  = base_path*ch_str*"-"*channel_to_name[ch]*"/"
        if !isdir(filename)
            mkpath(filename)
        end
        filename *= channel_to_name[ch] * "_w_noise" * split(basename(file), channel_to_name[ch])[2]
        if isfile(filename)
            if stat(filename).size <= stat(file).size/2
                create_file = true
            end
        else
            create_file = true
        end
        if create_file
            data_w_bl = Table(energy=[], multiplicity=[], timestamp=[], run=[], channel=[], waveform=[])

            data = HDF5.h5open(file, "r") do h5f
                LegendHDF5IO.readdata(h5f, "data")
            end

            for event in data
                time = event.waveform.time
                pulse = event.waveform.value ./ get_avg_maximum(event.waveform.value, 10)
                pulse .*= sum(event.energy).val

                bl_id = indices[pulse_id].bl
                first = indices[pulse_id].first

                last = first + length(pulse) -1
                pulse .+= bl[bl_id].waveform.value[first:1:last]

                append!(data_w_bl, Table(energy=[event.energy],
                                        multiplicity=[event.multiplicity],
                                        timestamp=[event.timestamp],
                                        run=[event.run],
                                        channel=[event.channel],
                                        waveform=[RDWaveform(time, pulse)]));
                pulse_id += 1
            end

            HDF5.h5open(filename, "w") do h5f
                LegendHDF5IO.writedata(h5f, "data", Table(energy        = VectorOfArrays(Array{typeof(data_w_bl[1].energy),1}(data_w_bl.energy)),
                                                        multiplicity = Int32.(data_w_bl.multiplicity),
                                                        timestamp    = Int32.(data_w_bl.timestamp),
                                                        run          = Int32.(data_w_bl.run),
                                                        channel      = Int32.(data_w_bl.channel),
                                                        waveform     = ArrayOfRDWaveforms( Array{typeof(data_w_bl[1].waveform), 1}(data_w_bl.waveform) )))
            end
        else
            pulse_id += size(HDF5.h5open(file, "r") do h5f
                LegendHDF5IO.readdata(h5f, "data")
            end, 1)

        end
#         IJulia.clear_output(true)
        Base.run(`clear`)
        @info(string(pulse_id)*" pulses done")
        next!(prog)
    end
end


text:
## Definitions of algorithms ##

# WrayHERK is scheme C in Liska and Colonius (JCP 2016)
# BraseyHairerHERK is scheme B in Liska and Colonius (JCP 2016)
# LiskaIFHERK is scheme A in Liska and Colonius (JCP 2016)

for (Alg,Order) in [(:WrayHERK,3),(:BraseyHairerHERK,3),(:LiskaIFHERK,2),(:IFHEEuler,1)]
    @eval struct $Alg{solverType} <: ConstrainedOrdinaryDiffEqAlgorithm
      maxiter :: Int
      tol :: Float64
    end

    @eval $Alg(;saddlesolver=Direct,maxiter=4,tol=eps(Float64)) = $Alg{saddlesolver}(maxiter,tol)

    @eval export $Alg

    @eval alg_order(alg::$Alg) = $Order
end

# LiskaIFHERK

@cache struct LiskaIFHERKCache{sc,ni,solverType,uType,rateType,expType1,expType2,saddleType,pType,TabType} <: ConstrainedODEMutableCache{sc,solverType}
  u::uType
  uprev::uType # qi
  k1::rateType # w1
  k2::rateType # w2
  k3::rateType # w3
  utmp::uType  # cache
  udiff::uType
  dutmp::rateType # cache for rates
  fsalfirst::rateType
  Hhalfdt::expType1
  Hzero::expType2
  S::saddleType
  ptmp::pType
  k::rateType
  tab::TabType
end

struct LiskaIFHERKConstantCache{sc,ni,solverType,T,T2} <: ConstrainedODEConstantCache{sc,solverType}
  ã11::T
  ã21::T
  ã22::T
  ã31::T
  ã32::T
  ã33::T
  c̃1::T2
  c̃2::T2
  c̃3::T2

  function LiskaIFHERKConstantCache{sc,ni,solverType}(T, T2) where {sc,ni,solverType}
    ã11 = T(1//2)
    ã21 = T(√3/3)
    ã22 = T((3-√3)/3)
    ã31 = T((3+√3)/6)
    ã32 = T(-√3/3)
    ã33 = T((3+√3)/6)
    c̃1 = T2(1//2)
    c̃2 = T2(1.0)
    c̃3 = T2(1.0)
    new{sc,ni,solverType,T,T2}(ã11,ã21,ã22,ã31,ã32,ã33,c̃1,c̃2,c̃3)
  end
end

LiskaIFHERKCache{sc,ni,solverType}(u,uprev,k1,k2,k3,utmp,udiff,dutmp,fsalfirst,
                                Hhalfdt,Hzero,S,ptmp,k,tab) where {sc,ni,solverType} =
        LiskaIFHERKCache{sc,ni,solverType,typeof(u),typeof(k1),typeof(Hhalfdt),typeof(Hzero),
                        typeof(S),typeof(ptmp),typeof(tab)}(u,uprev,k1,k2,k3,utmp,udiff,dutmp,fsalfirst,
                                                          Hhalfdt,Hzero,S,ptmp,k,tab)

function alg_cache(alg::LiskaIFHERK{solverType},u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,
                   tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Val{true}) where {solverType}

  u isa ArrayPartition || error("u must be of type ArrayPartition")

  y, z = state(u), constraint(u)

  utmp, udiff = (zero(u) for i in 1:2)
  k1, k2, k3, dutmp, fsalfirst, k = (zero(rate_prototype) for i in 1:6)

  sc = isstatic(f)
  ni = needs_iteration(f,u,p,rate_prototype)

  tab = LiskaIFHERKConstantCache{sc,ni,solverType}(constvalue(uBottomEltypeNoUnits),
                                                constvalue(tTypeNoUnits))

  L = _fetch_ode_L(f)
  Hhalfdt = exp(L,-dt/2,y)
  Hzero = exp(L,zero(dt),y)

  S = []
  push!(S,SaddleSystem(Hhalfdt,f,p,p,dutmp,solverType))
  push!(S,SaddleSystem(Hzero,f,p,p,dutmp,solverType))

  LiskaIFHERKCache{sc,ni,solverType}(u,uprev,k1,k2,k3,utmp,udiff,dutmp,fsalfirst,
                                  Hhalfdt,Hzero,S,deepcopy(p),k,tab)
end

function alg_cache(alg::LiskaIFHERK{solverType},u,rate_prototype,
                                  uEltypeNoUnits,uBottomEltypeNoUnits,
                                  tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,
                                  p,calck,::Val{false}) where {solverType}
  LiskaIFHERKConstantCache{isstatic(f),
                           needs_iteration(f,u,p,rate_prototype),
                           solverType}(constvalue(uBottomEltypeNoUnits),
                                          constvalue(tTypeNoUnits))
end

# IFHEEuler

@cache struct IFHEEulerCache{sc,ni,solverType,uType,rateType,expType,saddleType,pType} <: ConstrainedODEMutableCache{sc,solverType}
  u::uType
  uprev::uType # qi
  k1::rateType # w1
  utmp::uType  # cache
  udiff::uType
  dutmp::rateType # cache for rates
  fsalfirst::rateType
  Hdt::expType
  S::saddleType
  ptmp::pType
  k::rateType
end

struct IFHEEulerConstantCache{sc,ni,solverType} <: ConstrainedODEConstantCache{sc,solverType}

end

IFHEEulerCache{sc,ni,solverType}(u,uprev,k1,utmp,udiff,dutmp,fsalfirst,
                                Hdt,S,ptmp,k) where {sc,ni,solverType} =
        IFHEEulerCache{sc,ni,solverType,typeof(u),typeof(k1),typeof(Hdt),
                        typeof(S),typeof(ptmp)}(u,uprev,k1,utmp,udiff,dutmp,fsalfirst,
                                                              Hdt,S,ptmp,k)

function alg_cache(alg::IFHEEuler{solverType},u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,
                   tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Val{true}) where {solverType}

  u isa ArrayPartition || error("u must be of type ArrayPartition")

  y, z = state(u), constraint(u)

  utmp, udiff = (zero(u) for i in 1:2)
  k1, dutmp, fsalfirst, k = (zero(rate_prototype) for i in 1:4)

  sc = isstatic(f)
  ni = needs_iteration(f,u,p,rate_prototype)

  Hdt = exp(_fetch_ode_L(f),-dt,y)

  S = []
  push!(S,SaddleSystem(Hdt,f,p,p,dutmp,solverType))

  IFHEEulerCache{sc,ni,solverType}(u,uprev,k1,utmp,udiff,dutmp,fsalfirst,
                                  Hdt,S,deepcopy(p),k)
end

function alg_cache(alg::IFHEEuler{solverType},u,rate_prototype,
                                  uEltypeNoUnits,uBottomEltypeNoUnits,
                                  tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,
                                  p,calck,::Val{false}) where {solverType}
  IFHEEulerConstantCache{isstatic(f),needs_iteration(f,u,p,rate_prototype),solverType}()
end

#######

function initialize!(integrator,cache::LiskaIFHERKCache)
    @unpack k,fsalfirst = cache

    integrator.fsalfirst = fsalfirst
    integrator.fsallast = k
    integrator.kshortsize = 2
    resize!(integrator.k, integrator.kshortsize)
    integrator.k[1] = integrator.fsalfirst
    integrator.k[2] = integrator.fsallast
    integrator.f.odef(integrator.fsalfirst, integrator.uprev, integrator.p, integrator.t) # Pre-start fsal
    integrator.f.param_update_func(integrator.p,integrator.uprev,integrator.p,integrator.t)
    integrator.destats.nf += 1

end

@muladd function perform_step!(integrator,cache::LiskaIFHERKCache{sc,ni,solverType},repeat_step=false) where {sc,ni,solverType}
    @unpack t,dt,uprev,u,f,p,alg,opts = integrator
    @unpack internalnorm = opts
    @unpack maxiter, tol = alg
    @unpack k1,k2,k3,utmp,udiff,dutmp,fsalfirst,Hhalfdt,Hzero,S,ptmp,k = cache
    @unpack ã11,ã21,ã22,ã31,ã32,ã33,c̃1,c̃2,c̃3 = cache.tab
    @unpack param_update_func = f

    init_err = float(1)
    init_iter = ni ? 1 : maxiter

    # aliases to the state and constraint parts
    ytmp, ztmp, xtmp = state(utmp), constraint(utmp), aux_state(utmp)
    yprev = state(uprev)
    y, z, x = state(u), constraint(u), aux_state(u)
    pold_ptr = p
    pnew_ptr = ptmp

    ttmp = t
    u .= uprev

    ## Stage 1
    _ode_r1!(k1,f,u,pold_ptr,ttmp)
    integrator.destats.nf += 1
    @.. k1 *= dt*ã11
    @.. utmp = uprev + k1
    ttmp = t + dt*c̃1

    u .= utmp

    # if applicable, update p, construct new saddle system here, using Hhalfdt
    # and solve system. Solve iteratively if saddle operators depend on
    # constrained part of the state.
    err, numiter = init_err, init_iter
    while err > tol && numiter <= maxiter
      udiff .= u
      param_update_func(pnew_ptr,u,pold_ptr,ttmp)
      S[1] = SaddleSystem(S[1],Hhalfdt,f,pnew_ptr,pold_ptr,cache)
      _constraint_r2!(utmp,f,u,pnew_ptr,ttmp) # only updates the z part
      mainvector(u) .= S[1]\mainvector(utmp)
      @.. udiff -= u
      numiter += 1
      err = internalnorm(udiff,ttmp)
    end
    zero_vec!(xtmp)
    B1_times_z!(utmp,S[1])
    pold_ptr = ptmp
    pnew_ptr = p

    ldiv!(yprev,Hhalfdt,yprev)
    ldiv!(state(k1),Hhalfdt,state(k1))

    @.. k1 = (k1-utmp)/(dt*ã11)  # r1(y,t) - B1T*z

    ## Stage 2
    _ode_r1!(k2,f,u,pold_ptr,ttmp)
    integrator.destats.nf += 1
    @.. k2 *= dt*ã22
    @.. utmp = uprev + k2 + dt*ã21*k1
    ttmp = t + dt*c̃2

    u .= utmp

    # if applicable, update p, construct new saddle system here, using Hhalfdt
    err, numiter = init_err, init_iter
    while err > tol && numiter <= maxiter
      udiff .= u
      param_update_func(pnew_ptr,u,pold_ptr,ttmp)
      S[1] = SaddleSystem(S[1],Hhalfdt,f,pnew_ptr,pold_ptr,cache)
      _constraint_r2!(utmp,f,u,pnew_ptr,ttmp)
      mainvector(u) .= S[1]\mainvector(utmp)
      @.. udiff -= u
      numiter += 1
      err = internalnorm(udiff,ttmp)
    end
    zero_vec!(xtmp)
    B1_times_z!(utmp,S[1])
    pold_ptr = p
    pnew_ptr = ptmp

    ldiv!(yprev,Hhalfdt,yprev)
    ldiv!(state(k1),Hhalfdt,state(k1))
    ldiv!(state(k2),Hhalfdt,state(k2))

    @.. k2 = (k2-utmp)/(dt*ã22)

    ## Stage 3
    _ode_r1!(k3,f,u,pold_ptr,ttmp)
    integrator.destats.nf += 1
    @.. k3 *= dt*ã33
    @.. utmp = uprev + k3 + dt*ã32*k2 + dt*ã31*k1
    ttmp = t + dt

    u .= utmp

    # if applicable, update p, construct new saddle system here, using Hzero (identity)
    err, numiter = init_err, init_iter
    while err > tol && numiter <= maxiter
      udiff .= u
      param_update_func(pnew_ptr,u,pold_ptr,ttmp)
      S[2] = SaddleSystem(S[2],Hzero,f,pnew_ptr,pold_ptr,cache)
      _constraint_r2!(utmp,f,u,pnew_ptr,t+dt)
      mainvector(u) .= S[2]\mainvector(utmp)
      @.. udiff -= u
      numiter += 1
      err = internalnorm(udiff,ttmp)
      #println("error = ",err)
    end
    @.. z /= (dt*ã33)

    # Final steps
    param_update_func(p,u,ptmp,t+dt)
    f.odef(integrator.fsallast, u, p, t+dt)
    integrator.destats.nf += 1

    return nothing
end

function initialize!(integrator,cache::LiskaIFHERKConstantCache)
  integrator.kshortsize = 2
  integrator.k = typeof(integrator.k)(undef, integrator.kshortsize)
  integrator.fsalfirst = integrator.f.odef(integrator.uprev, integrator.p, integrator.t) # Pre-start fsal
  integrator.p = integrator.f.param_update_func(integrator.uprev,integrator.p,integrator.t)
  integrator.destats.nf += 1

  # Avoid undefined entries if k is an array of arrays
  integrator.fsallast = zero(integrator.fsalfirst)
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
end

@muladd function perform_step!(integrator,cache::LiskaIFHERKConstantCache{sc,ni,solverType},repeat_step=false) where {sc,ni,solverType}
    @unpack t,dt,uprev,f,p,opts,alg = integrator
    @unpack internalnorm = opts
    @unpack maxiter, tol = alg
    @unpack ã11,ã21,ã22,ã31,ã32,ã33,c̃1,c̃2,c̃3 = cache
    @unpack param_update_func = f

    init_err = float(1)
    init_iter = ni ? 1 : maxiter

    # set up some cache variables
    yprev = state(uprev)
    udiff = deepcopy(uprev)
    ducache = deepcopy(uprev)
    ptmp = deepcopy(p)
    L = _fetch_ode_L(f)
    Hhalfdt = exp(L,-dt/2,state(uprev))
    Hzero = exp(L,zero(dt),state(uprev))
    pold_ptr = p

    ## Stage 1
    ttmp = t
    k1 = _ode_r1(f,uprev,pold_ptr,ttmp)
    integrator.destats.nf += 1
    @.. k1 *= dt*ã11
    utmp = @.. uprev + k1
    ttmp = t + dt*c̃1

    # if applicable, update p, construct new saddle system here, using Hhalfdt
    # and solve system. Solve iteratively if saddle operators depend on
    # constrained part of the state.
    err, numiter = init_err, init_iter
    u = deepcopy(utmp)
    while err > tol && numiter <= maxiter
      udiff .= u
      ptmp = param_update_func(u,pold_ptr,ttmp)
      pnew_ptr = ptmp
      S = SaddleSystem(Hhalfdt,f,pnew_ptr,pold_ptr,ducache,solverType)
      constraint(utmp) .= constraint(_constraint_r2(f,u,pnew_ptr,ttmp))
      mainvector(u) .= S\mainvector(utmp)
      B1_times_z!(ducache,S)

      @.. udiff -= u
      numiter += 1
      err = internalnorm(udiff,ttmp)
    end
    zero_vec!(aux_state(utmp))
    state(utmp) .= state(ducache)
    pold_ptr = ptmp

    ldiv!(yprev,Hhalfdt,yprev)
    ldiv!(state(k1),Hhalfdt,state(k1))

    @.. k1 = (k1-utmp)/(dt*ã11)  # r1(y,t) - B1T*z

    ## Stage 2
    k2 = _ode_r1(f,u,pold_ptr,ttmp)
    integrator.destats.nf += 1
    @.. k2 *= dt*ã22
    @.. utmp = uprev + k2 + dt*ã21*k1
    ttmp = t + dt*c̃2

    # if applicable, update p, construct new saddle system here, using Hhalfdt
    err, numiter = init_err, init_iter
    u .= utmp
    while err > tol && numiter <= maxiter
      udiff .= u
      ptmp = param_update_func(u,pold_ptr,ttmp)
      S = SaddleSystem(Hhalfdt,f,ptmp,pold_ptr,ducache,solverType)
      constraint(utmp) .= constraint(_constraint_r2(f,u,ptmp,ttmp))
      mainvector(u) .= S\mainvector(utmp)
      B1_times_z!(ducache,S)
      @.. udiff -= u
      numiter += 1
      err = internalnorm(udiff,ttmp)
    end
    zero_vec!(aux_state(utmp))
    state(utmp) .= state(ducache)
    pold_ptr = ptmp

    ldiv!(yprev,Hhalfdt,yprev)
    ldiv!(state(k1),Hhalfdt,state(k1))
    ldiv!(state(k2),Hhalfdt,state(k2))

    @.. k2 = (k2-utmp)/(dt*ã22)

    ## Stage 3
    k3 = _ode_r1(f,u,pold_ptr,ttmp)
    integrator.destats.nf += 1
    @.. k3 *= dt*ã33
    @.. utmp = uprev + k3 + dt*ã32*k2 + dt*ã31*k1
    ttmp = t + dt

    # if applicable, update p, construct new saddle system here, using Hzero (identity)
    err, numiter = init_err, init_iter
    u .= utmp
    while err > tol && numiter <= maxiter
      udiff .= u
      ptmp = param_update_func(u,pold_ptr,ttmp)
      S = SaddleSystem(Hzero,f,ptmp,pold_ptr,ducache,solverType)
      constraint(utmp) .= constraint(_constraint_r2(f,u,ptmp,ttmp))
      mainvector(u) .= S\mainvector(utmp)
      @.. udiff -= u
      numiter += 1
      err = internalnorm(udiff,ttmp)
      #println("error = ",err)
    end
    z = constraint(u)
    @.. z /= (dt*ã33)

    # Final steps
    integrator.p = param_update_func(u,ptmp,t)
    k = f.odef(u, integrator.p, t+dt)
    integrator.destats.nf += 1

    integrator.fsallast = k
    integrator.k[1] = integrator.fsalfirst
    integrator.k[2] = integrator.fsallast
    integrator.u = u

    return nothing
end

####

function initialize!(integrator,cache::IFHEEulerCache)
    @unpack k,fsalfirst = cache

    integrator.fsalfirst = fsalfirst
    integrator.fsallast = k
    integrator.kshortsize = 2
    resize!(integrator.k, integrator.kshortsize)
    integrator.k[1] = integrator.fsalfirst
    integrator.k[2] = integrator.fsallast
    integrator.f.odef(integrator.fsalfirst, integrator.uprev, integrator.p, integrator.t) # Pre-start fsal
    integrator.f.param_update_func(integrator.p,integrator.uprev,integrator.p,integrator.t)
    integrator.destats.nf += 1

end

@muladd function perform_step!(integrator,cache::IFHEEulerCache{sc,ni,solverType},repeat_step=false) where {sc,ni,solverType}
    @unpack t,dt,uprev,u,f,p,opts,alg = integrator
    @unpack internalnorm = opts
    @unpack k1,utmp,udiff,dutmp,fsalfirst,Hdt,S,ptmp,k = cache
    @unpack maxiter, tol = alg
    @unpack param_update_func = f

    init_err = float(1)
    #init_iter = ni ? 1 : maxiter
    init_iter = maxiter  # First-order method does not need iteration

    # aliases to the state and constraint parts
    ytmp, ztmp = state(utmp), constraint(utmp)
    z = constraint(u)
    pold_ptr = p
    pnew_ptr = ptmp

    ttmp = t
    u .= uprev

    _ode_r1!(k1,f,u,pold_ptr,ttmp)
    integrator.destats.nf += 1
    @.. k1 *= dt
    @.. utmp = uprev + k1
    ttmp = t + dt

    # if applicable, update p, construct new saddle system here, using Hdt
    err, numiter = init_err, init_iter
    u .= utmp
    while err > tol && numiter <= maxiter
      udiff .= u
      param_update_func(pnew_ptr,u,pold_ptr,ttmp)
      S[1] = SaddleSystem(S[1],Hdt,f,pnew_ptr,pold_ptr,cache)
      _constraint_r2!(utmp,f,u,pnew_ptr,t+dt)
      mainvector(u) .= S[1]\mainvector(utmp)
      @.. udiff -= u
      numiter += 1
      err = internalnorm(udiff,ttmp)
      #println("error = ",err)
    end
    @.. z /= dt

    # Final steps
    param_update_func(p,u,pold_ptr,t)
    f.odef(integrator.fsallast, u, p, t+dt)
    integrator.destats.nf += 1

    return nothing
end

function initialize!(integrator,cache::IFHEEulerConstantCache)
  integrator.kshortsize = 2
  integrator.k = typeof(integrator.k)(undef, integrator.kshortsize)
  integrator.fsalfirst = integrator.f.odef(integrator.uprev, integrator.p, integrator.t) # Pre-start fsal
  integrator.p = integrator.f.param_update_func(integrator.uprev,integrator.p,integrator.t)
  integrator.destats.nf += 1

  # Avoid undefined entries if k is an array of arrays
  integrator.fsallast = zero(integrator.fsalfirst)
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
end

@muladd function perform_step!(integrator,cache::IFHEEulerConstantCache{sc,ni,solverType},repeat_step=false) where {sc,ni,solverType}
  @unpack t,dt,uprev,f,p,opts,alg = integrator
  @unpack internalnorm = opts
  @unpack maxiter, tol = alg
  @unpack param_update_func = f

  init_err = float(1)
  #init_iter = ni ? 1 : maxiter
  init_iter = maxiter  # First-order method does not need iteration

  # set up some cache variables
  udiff = deepcopy(uprev)
  ducache = deepcopy(uprev)
  ptmp = deepcopy(p)
  L = _fetch_ode_L(f)
  Hdt = exp(L,-dt,state(uprev))
  pold_ptr = p
  pnew_ptr = ptmp

  k1 = _ode_r1(f,uprev,pold_ptr,t)
  integrator.destats.nf += 1
  @.. k1 *= dt
  utmp = @.. uprev + k1

  # if applicable, update p, construct new saddle system here, using Hdt
  err, numiter = init_err, init_iter
  u = deepcopy(utmp)
  while err > tol && numiter <= maxiter
    udiff .= u
    pnew_ptr = param_update_func(u,pold_ptr,t+dt)
    S = SaddleSystem(Hdt,f,pnew_ptr,pold_ptr,ducache,solverType)
    constraint(utmp) .= constraint(_constraint_r2(f,u,pnew_ptr,t+dt))
    mainvector(u) .= S\mainvector(utmp)
    @.. udiff -= u
    numiter += 1
    err = internalnorm(udiff,t+dt)
    #println("error = ",err)
  end
  z = constraint(u)
  @.. z /= dt

  # Final steps
  integrator.p = param_update_func(u,pold_ptr,t)
  k = f.odef(u, integrator.p, t+dt)
  integrator.destats.nf += 1

  integrator.fsallast = k
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
  integrator.u = u
end


text:
using JuMP
using Clp
COSTS = [80 215; 100 108; 102 68]
LIMITS = [1000, 1500, 1200]
DEMANDS = [2300 1400]
m = Model(Clp.Optimizer)
@variable(m, x[1:3, 1:2] ≥ 0)
@objective(m, Min, sum(COSTS .* x))
@constraint(m, sum(x, dims=2) .≤ LIMITS)
@constraint(m, sum(x, dims=1) .== DEMANDS)
print(m)
optimize!(m)
status = termination_status(m)
println("Solution status: ", status)
println("Objective value: ", objective_value(m))


text:
let
	# Test Optim.cg for all differentiable functions in Optim.UnconstrainedProblems.examples
	for (name, prob) in Optim.UnconstrainedProblems.examples
		if prob.isdifferentiable
			df = DifferentiableFunction(prob.f, prob.g!)
			res = Optim.optimize(df, prob.initial_x, ConjugateGradient())
				@assert norm(Optim.minimizer(res) - prob.solutions) < 1e-2
		end
	end

	let
	objective(X, B) = sum((X.-B).^2)/2

	function objective_gradient!(X, G, B)
	    for i = 1:length(G)
	        G[i] = X[i]-B[i]
	    end
	end

	srand(1)
	B = rand(2,2)
	df = Optim.DifferentiableFunction(X -> objective(X, B), (X, G) -> objective_gradient!(X, G, B))
	results = Optim.optimize(df, rand(2,2), ConjugateGradient())
	@assert Optim.converged(results)
	@assert Optim.minimum(results) < 1e-8
	end
end


text:
# Helpers
id_or_chain(obj::AbstractString, kind::AbstractString) = kind * "id=" * obj * "&"
id_or_chain(obj::AbstractVector{<:AbstractString}, kind::AbstractString) = reduce(*, kind .* "id=" .* obj .* "&")
function period(startdate::TimeType, enddate::TimeType, dataset::AbstractString)
    starts = string.(startdate:ifelse(dataset ∈ ["GSOM", "GSOY"], Year(10), Year(1)):enddate)
    return "&startdate=" .* starts .* "&enddate=" .* push!(starts[2:end], string(enddate))
end
for_assignment(::Type{Union{Missing, T}}, value::Missing) where T <: Union{AbstractString, Number, TimeType} = missing
for_assignment(::Type{Union{Missing, T}}, value::AbstractString) where T <: AbstractString = convert(T, value)
for_assignment(::Type{Union{Missing, T}}, value::Number) where T <: Number = convert(T, value)
for_assignment(::Type{Union{Missing, T}}, value::AbstractString) where T <: TimeType = T(convert(String, value))
for_assignment(::Type{T}, value::Missing) where T <: Union{AbstractString, Number, TimeType} = missing
for_assignment(::Type{T}, value::AbstractString) where T <: AbstractString = convert(T, value)
for_assignment(::Type{T}, value::Number) where T <: Number = convert(T, value)
for_assignment(::Type{T}, value::AbstractString) where T <: TimeType = T(convert(String, value))

isvalid_cdotoken(obj::AbstractString) = occursin(r"[A-Za-z]{32}", obj)
isvalid_datsets(obj::AbstractVector{<:AbstractString}) = all(isvalid_datsets, obj)
isvalid_datsets(obj::AbstractString) = obj ∈ DATASETS
isvalid_datacategories(obj::AbstractVector{<:AbstractString}) = all(isvalid_datacategories, obj)
isvalid_datacategories(obj::AbstractString) = obj ∈ DATACATEGORIES
isvalid_datatypes(obj::AbstractVector{<:AbstractString}) = all(elem -> occursin(r"^[A-Z0-9-]{3,}$", elem), obj)
isvalid_datatypes(obj::AbstractString) = occursin(r"^[A-Z0-9-]{3,}$", obj)
isvalid_locationcategories(obj::AbstractVector{<:AbstractString}) = all(isvalid_locationcategories, obj)
isvalid_locationcategories(obj::AbstractString) = obj ∈ LOCATIONCATEGORIES
isvalid_locations(obj::AbstractVector{<:AbstractString}) = all(isvalid_locations, obj)
isvalid_locations(obj::AbstractString) = occursin(r"^(CITY|CLIM_DIV|CLIM_REG|CNTRY|CNTY|FIPS|HYD_ACC|HYD_CAT|HYD_REG|HYD_SUB|ST|US_TERR|ZIP):[A-Z0-9]{1,}$", obj)
isvalid_stations(obj::AbstractVector{<:AbstractString}) = all(isvalid_stations, obj)
isvalid_stations(obj::AbstractString) = occursin(r"^(COOP|GHCND|NEXRAD|WBAN):[a-zA-Z0-9-_]{1,}$", obj)
isvalid_extent(obj::AbstractVector{<:AbstractFloat}) = length(obj) ∈ [0, 4]

# Struct
abstract type Endpoint end
struct CDO_Data <: Endpoint
    token::String
    url::Vector{String}
    function CDO_Data(CDO_token::AbstractString,
                      dataset::AbstractString,
                      startdate::TimeType,
                      enddate::TimeType,
                      datatypes::Union{AbstractString, AbstractVector{<:AbstractString}},
                      locations::Union{AbstractString, AbstractVector{<:AbstractString}},
                      stations::Union{AbstractString, AbstractVector{<:AbstractString}},
                      metric::Bool)
        url = "https://www.ncdc.noaa.gov/cdo-web/api/v2/data?datasetid=" * dataset * "&" *
              id_or_chain(datatypes, "datatype") *
              id_or_chain(locations, "location") *
              id_or_chain(stations, "station") *
              "units=" * ifelse(metric, "metric", "standard") .*
              period(startdate, enddate, dataset) .*
              "&limit=1000&offset=1"
        return new(CDO_token, url)
    end
end
struct CDO_DataCategory <: Endpoint
    token::String
    url::String
    CDO_DataCategory(CDO_token::AbstractString, datacategory::AbstractString) = new(CDO_token, "https://www.ncdc.noaa.gov/cdo-web/api/v2/datacategories/" * datacategory)
end
struct CDO_DataCategories <: Endpoint
    token::String
    url::String
    function CDO_DataCategories(CDO_token::AbstractString,
                                datasets::Union{AbstractString, AbstractVector{<:AbstractString}},
                                locations::Union{AbstractString, AbstractVector{<:AbstractString}},
                                stations::Union{AbstractString, AbstractVector{<:AbstractString}},
                                startdate::Date,
                                enddate::Date)
        return new(CDO_token,
                   "https://www.ncdc.noaa.gov/cdo-web/api/v2/datacategories?" *
                   id_or_chain(datasets, "dataset") *
                   id_or_chain(locations, "location") *
                   id_or_chain(stations, "station") *
                   "startdate=" * string(startdate) *
                   "&enddate=" * string(enddate) *
                   "&limit=1000&offset=1")
    end
end
struct CDO_Dataset <: Endpoint
    token::String
    url::String
    CDO_Dataset(CDO_token::AbstractString, dataset::AbstractString) = new(CDO_token, "https://www.ncdc.noaa.gov/cdo-web/api/v2/datasets/" * dataset)
end
struct CDO_Datasets <: Endpoint
    token::String
    url::String
    function CDO_Datasets(CDO_token::AbstractString,
                          datatypes::Union{AbstractString, AbstractVector{<:AbstractString}},
                          locations::Union{AbstractString, AbstractVector{<:AbstractString}},
                          stations::Union{AbstractString, AbstractVector{<:AbstractString}},
                          startdate::Date,
                          enddate::Date)
        return new(CDO_token,
                   "https://www.ncdc.noaa.gov/cdo-web/api/v2/datasets?" *
                   id_or_chain(datatypes, "datatype") *
                   id_or_chain(locations, "location") *
                   id_or_chain(stations, "station") *
                   "startdate=" * string(startdate) *
                   "&enddate=" * string(enddate) *
                   "&limit=1000&offset=1")
    end
end
struct CDO_DataType <: Endpoint
    token::String
    url::String
    CDO_DataType(CDO_token::AbstractString, datatype::AbstractString) = new(CDO_token, "https://www.ncdc.noaa.gov/cdo-web/api/v2/datatypes/" * datatype)
end
struct CDO_DataTypes <: Endpoint
    token::String
    url::String
    function CDO_DataTypes(CDO_token::AbstractString,
                           datasets::Union{AbstractString, AbstractVector{<:AbstractString}},
                           datacategories::Union{AbstractString, AbstractVector{<:AbstractString}},
                           locations::Union{AbstractString, AbstractVector{<:AbstractString}},
                           stations::Union{AbstractString, AbstractVector{<:AbstractString}},
                           startdate::Date,
                           enddate::Date)
        return new(CDO_token,
                   "https://www.ncdc.noaa.gov/cdo-web/api/v2/datatypes?" *
                   id_or_chain(datasets, "dataset") *
                   id_or_chain(datacategories, "datacategory") *
                   id_or_chain(locations, "location") *
                   id_or_chain(stations, "station") *
                   "startdate=" * string(startdate) *
                   "&enddate=" * string(enddate) *
                   "&limit=1000&offset=1")
        end
end
struct CDO_LocationCategory <: Endpoint
    token::String
    url::String
    CDO_LocationCategory(token::AbstractString, locationcategory::AbstractString) = new(token, "https://www.ncdc.noaa.gov/cdo-web/api/v2/locationcategories/" * locationcategory)
end
struct CDO_LocationCategories <: Endpoint
    token::String
    url::String
    function CDO_LocationCategories(CDO_token::AbstractString,
                                    datasets::Union{AbstractString, AbstractVector{<:AbstractString}},
                                    startdate::Date,
                                    enddate::Date)
        return new(CDO_token,
                   "https://www.ncdc.noaa.gov/cdo-web/api/v2/locationcategories?" *
                   id_or_chain(datasets, "dataset") *
                   "startdate=" * string(startdate) *
                   "&enddate=" * string(enddate) *
                   "&limit=1000&offset=1")
    end
end
struct CDO_Location <: Endpoint
    token::String
    url::String
    CDO_Location(CDO_token::AbstractString, location::AbstractString) = new(CDO_token, "https://www.ncdc.noaa.gov/cdo-web/api/v2/locations/" * location)
end
struct CDO_Locations <: Endpoint
    token::String
    url::String
    function CDO_Locations(CDO_token::AbstractString,
                           datasets::Union{AbstractString, AbstractVector{<:AbstractString}},
                           locationcategories::Union{AbstractString, AbstractVector{<:AbstractString}},
                           datacategories::Union{AbstractString, AbstractVector{<:AbstractString}},
                           startdate::Date,
                           enddate::Date)
        return new(CDO_token,
                   "https://www.ncdc.noaa.gov/cdo-web/api/v2/locations?" *
                   id_or_chain(datasets, "dataset") *
                   id_or_chain(locationcategories, "locationcategory") *
                   id_or_chain(datacategories, "datacategory") *
                   "startdate=" * string(startdate) *
                   "&enddate=" * string(enddate) *
                   "&limit=1000&offset=1")
        end
end
struct CDO_Station <: Endpoint
    token::String
    url::String
    CDO_Station(CDO_token::AbstractString, station::AbstractString) = new(CDO_token, "https://www.ncdc.noaa.gov/cdo-web/api/v2/stations/" * station)
end
struct CDO_Stations <: Endpoint
    token::String
    url::String
    function CDO_Stations(CDO_token::AbstractString,
                          datasets::Union{AbstractString, AbstractVector{<:AbstractString}},
                          locations::Union{AbstractString, AbstractVector{<:AbstractString}},
                          datacategories::Union{AbstractString, AbstractVector{<:AbstractString}},
                          datatypes::Union{AbstractString, AbstractVector{<:AbstractString}},
                          extent::Union{AbstractVector{<:AbstractFloat}},
                          startdate::Date,
                          enddate::Date)
        return new(CDO_token,
                   "https://www.ncdc.noaa.gov/cdo-web/api/v2/stations?" *
                   id_or_chain(datasets, "dataset") *
                   id_or_chain(locations, "location") *
                   id_or_chain(datacategories, "datacategory") *
                   id_or_chain(datatypes, "datatype") *
                   "extend=" * reduce(*, string.(extent) .* ",")[1:end - 1] *
                   "&startdate=" * string(startdate) *
                   "&enddate=" * string(enddate) *
                   "&limit=1000&offset=1")
    end
end

# Constants
const CDO_Single = Union{CDO_DataCategory, CDO_Dataset, CDO_DataType, CDO_LocationCategory, CDO_Location, CDO_Station}
const CDO_Meta = Union{CDO_Data, CDO_DataCategories, CDO_Datasets, CDO_DataTypes, CDO_LocationCategories, CDO_Locations, CDO_Stations}
const DATASETS = ["GHCND", "GSOM", "GSOY", "NEXRAD2", "NEXRAD3", "NORMAL_ANN",
                  "NORMAL_DLY", "NORMAL_HLY", "NORMAL_MLY", "PRECIP_15", "PRECIP_HLY"]
const LOCATIONCATEGORIES = ["CITY", "CLIM_DIV", "CLIM_REG", "CNTRY", "CNTY", "HYD_ACC",
                            "HYD_CAT", "HYD_REG", "HYD_SUB", "ST", "US_TERR", "ZIP"]
const DATACATEGORIES = ["ANNAGR", "ANNDD", "ANNPRCP", "ANNTEMP", "AUAGR", "AUDD",
                        "AUPRCP", "AUTEMP", "COMP", "COMPAGR", "DD", "DUALPOLMOMENT",
                        "ECHOTOP", "EVAP", "HYDROMETEOR", "LAND", "MISC", "OTHER", "OVERLAY",
                        "PRCP", "PRES", "REFLECTIVITY", "SKY", "SPAGR", "SPDD", "SPPRCP", "SPTEMP",
                        "SUAGR", "SUDD", "SUN", "SUPRCP", "SUTEMP", "TEMP", "VELOCITY", "VERTINTLIQUID",
                        "WATER", "WIAGR", "WIDD", "WIND", "WIPRCP", "WITEMP", "WXTYPE"]
const CDO_NonValidToken = ArgumentError("The CDO_token is not valid.")
const CDO_NonValidDatasets = ArgumentError("The provided datasets argument is not valid. Valid datasets are: " * reduce(*, DATASETS .* ", ")[1:end - 2] * ".")
const CDO_NonValidDataCategories = ArgumentError("The data category argument is not valid. Valid data categories are: " * reduce(*, DATACATEGORIES .* ", ")[1:end - 2] * ".")
const CDO_NonValidDataTypes = ArgumentError("The provided datatypes argument is not valid. A valid data type has form r\"^[A-Z0-9-]{3,}\$\"")
const CDO_NonValidLocationCategories = ArgumentError("The provided location categories argument is not valid. Valid location categories are: CITY, CLIM_DIV, CLIM_REG, CNTRY, CNTY, HYD_ACC, HYD_CAT, HYD_REG, HYD_SUB, ST, US_TERR, ZIP.")
const CDO_NonValidLocations = ArgumentError("The provided locations argument is not valid. A valid location has form r\"^(CITY|CLIM_DIV|CLIM_REG|CNTRY|CNTY|FIPS|HYD_ACC|HYD_CAT|HYD_REG|HYD_SUB|ST|US_TERR|ZIP):[A-Z0-9]{1,}\$\"")
const CDO_NonValidStations = ArgumentError("The provided locations argument is not valid. A valid location has form r\"(COOP|GHCND|NEXRAD|WBAN):[a-zA-Z0-9-]{1,}\$\"")
const CDO_NonValidExtent = ArgumentError("Extent must be of length 0 or 4.")

const Names = [[:date, :datatype, :station, :attributes, :value],
               [:name, :id],
               [:mindate, :maxdate, :name, :datacoverage, :id],
               [:uid, :mindate, :maxdate, :name, :datacoverage, :id],
               [:mindate, :maxdate, :datacoverage, :id],
               [:mindate, :maxdate, :name, :datacoverage, :id],
               [:elevation, :mindate, :maxdate, :latitude, :name, :datacoverage, :id, :elevationUnit, :longitude]
               ]
const Types = [[DateTime, String, String, Union{String, Missing}, Float64],
               [String, String],
               [Date, Date, String, Union{Float64, Missing}, String],
               [String, Date, Date, String, Float64, String],
               [Date, Date, Union{Float64, Missing}, String],
               [Date, Date, String, Float64, String],
               [Union{Float64, Missing}, Date, Date, Float64, String, Float64, String, Union{String, Missing}, Float64]
               ]

"""
    names(::Endpoint)::Vector{Symbol}

Return the names of the values from a json values stream in order for the endpoint.
"""
names(::Endpoint) = error("names is not defined for $(typeof(obj)).")
names(::CDO_Data) = Names[1]
names(::CDO_DataCategory) = Names[2]
names(::CDO_DataCategories) = Names[2]
names(::CDO_Dataset) = Names[3]
names(::CDO_Datasets) = Names[4]
names(::CDO_DataType) = Names[5]
names(::CDO_DataTypes) = Names[6]
names(::CDO_LocationCategory) = Names[2]
names(::CDO_LocationCategories) = Names[2]
names(::CDO_Location) = Names[6]
names(::CDO_Locations) = Names[6]
names(::CDO_Station) = Names[7]
names(::CDO_Stations) = Names[7]

"""
    types(::Endpoint)::Vector{Symbol}

Return the types of the values from a json values stream in order for the endpoint.
"""
types(::Endpoint) = error("types is not defined for $(typeof(obj)).")
types(::CDO_Data) = Types[1]
types(::CDO_DataCategory) = Types[2]
types(::CDO_DataCategories) = Types[2]
types(::CDO_Dataset) = Types[3]
types(::CDO_Datasets) = Types[4]
types(::CDO_DataType) = Types[5]
types(::CDO_DataTypes) = Types[3]
types(::CDO_LocationCategory) = Types[2]
types(::CDO_LocationCategories) = Types[2]
types(::CDO_Location) = Types[6]
types(::CDO_Locations) = Types[6]
types(::CDO_Station) = Types[7]
types(::CDO_Stations) = Types[7]

"""
    exception(::Endpoint)::DataFrame
    exception(::Endpoint)::ArgumentError(msg::AbstractString)

In case of no results for a data request, return either an empty DataFrame of
the expected type indicating no data for the query parameters or an ArgumentError
if it is the result can be deduced to be an error in the argument.
"""
exception(obj::Endpoint) = DataFrame(types(obj), names(obj), 0)
exception(obj::CDO_DataCategory) = ArgumentError(getfield(obj, :url)[57:end] * " is not a valid data category. For a complete list of valid data categories run `cdo_datacategories(CDO_token::AbstractString)`.")
exception(obj::CDO_Dataset) = ArgumentError(getfield(obj, :url)[51:end] * " is not a valid dataset. For a complete list of valid datasets run `cdo_datasets(CDO_token::AbstractString)`.")
exception(obj::CDO_DataType) = ArgumentError(getfield(obj, :url)[52:end] * " is not a valid data type. For a complete list of valid data types run `cdo_datatypes(CDO_token::AbstractString)`.")
exception(obj::CDO_LocationCategory) = ArgumentError(getfield(obj, :url)[61:end] * " is not a valid location category. For a complete list of valid location categories run `cdo_locationcategories(CDO_token::AbstractString)`.")
exception(obj::CDO_Location) = ArgumentError(getfield(obj, :url)[52:end] * " is not a valid location. For a complete list of valid locations run `cdo_locations(CDO_token::AbstractString)`.")
exception(obj::CDO_Station) = ArgumentError(getfield(obj, :url)[51:end] * " is not a valid weather station. For a complete list of valid stations run `cdo_stations(CDO_token::AbstractString)`.")

"""
    skeleton(endpoint::CDO_Single, jsontext::AbstractString)

Return a dataframe with the parse values of a JSON according to the endpoint.
In case of no observations an empty dataframe is returned.
If an empty JSON indicates an ArgumentError this one is throw.
"""
function skeleton(endpoint::CDO_Single, jsontext::AbstractString)
    jsontext == "{}" && return exception(endpoint)
    Ts, Ns = (types(endpoint), names(endpoint))
    Keys = string.(Ns)
    output = DataFrame(Ts, Ns, 1)
    vals = value(jsontext)
    for (col, T, key) ∈ zip(Ns, Ts, Keys)
        output[col] = for_assignment(T, get(vals, key, missing))
    end
    return output
end

"""
    parse(obj::Endpoint)

Return the results for the endpoint.
"""
parse(obj::Endpoint) = error("parse is not defined for $(typeof(obj)).")
function parse(obj::CDO_Single)
    response = request("GET", getfield(obj, :url), [("token", getfield(obj, :token))])
    jsontext = String(getfield(response, :body))
    return skeleton(obj, jsontext)
end
function parse(obj::CDO_Meta)
    header = [("token", getfield(obj, :token))]
    Ts, Ns = (types(obj), names(obj))
    Keys = string.(Ns)
    url = getfield(obj, :url)
    response = request("GET", url, header)
    jsontext = String(getfield(response, :body))
    jsontext == "{}" && return DataFrame(Ts, Ns, 0)
    json = value(jsontext)
    Count = Int64(json["metadata"]["resultset"]["count"])
    urls = url[1:end - 1] .* string.(1:1000:Count)[2:end]
    idx = 1
    output = DataFrame(Ts, Ns, Count)
    for elem ∈ json["results"]
        for (col, T, key) ∈ zip(Ns, Ts, Keys)
            output[idx, col] = for_assignment(T, get(elem, key, missing))
        end
        idx += 1
    end
    for suburl ∈ urls
        response = request("GET", suburl, header)
        jsontext = String(getfield(response, :body))
        json = value(jsontext)["results"]
        for elem ∈ json
            for (col, T, key) ∈ zip(Ns, Ts, Keys)
                output[idx, col] = for_assignment(T, get(elem, key, missing))
            end
            idx += 1
        end
    end
    return output
end
function parse(obj::CDO_Data)
    header = [("token", getfield(obj, :token))]
    urls = getfield(obj, :url)
    Ns, Ts = (names(obj), types(obj))
    Keys = string.(Ns)
    output = DataFrame(Ts, Ns, 0)
    for url ∈ urls
        response = request("GET", url, header)
        jsontext = String(getfield(response, :body))
        jsontext == "{}" && return DataFrame(Ts, Ns, 0)
        json = value(jsontext)
        Count = Int64(json["metadata"]["resultset"]["count"])
        urls = url[1:end - 1] .* string.(1:1000:Count)[2:end] .* "&includemetadata=false"
        tmp = DataFrame(Ts, Ns, Count)
        idx = 1
        for elem ∈ json["results"]
            for (col, T, key) ∈ zip(Ns, Ts, Keys)
                tmp[idx, col] = for_assignment(T, get(elem, key, missing))
            end
            idx += 1
        end
        for suburl ∈ urls
            response = request("GET", suburl, header)
            jsontext = String(getfield(response, :body))
            json = value(jsontext)["results"]
            for elem ∈ json
                for (col, T, key) ∈ zip(Ns, Ts, Keys)
                    tmp[idx, col] = for_assignment(T, get(elem, key, missing))
                end
                idx += 1
            end
        end
        append!(output, tmp)
    end
    return output
end