text:
"""
    simple_extr_val_sel(data::ClustData,
                        extreme_value_descr_ar::Array{SimpleExtremeValueDescr,1};
                        rep_mod_method::String="feasibility")

Selects simple extreme values and returns modified data, extreme values, and the corresponding indices.

Inputs options for `rep_mod_method`:
- `rep_mod_method`::String : `feasibility`,`append`
"""
function simple_extr_val_sel(data::ClustData,
                             extr_value_descr_ar::Array{SimpleExtremeValueDescr,1};
                             rep_mod_method::String="feasibility"
                             )
  idcs = simple_extr_val_ident(data,extr_value_descr_ar)
  extr_vals = extreme_val_output(data,idcs;rep_mod_method=rep_mod_method)
  # for append method: modify data to be clustered to only contain the values that are not extreme values
  if rep_mod_method=="feasibility"
    data_mod = data
  elseif rep_mod_method=="append"
    data_mod = input_data_modification(data,idcs)
  else
    error("rep_mod_method - "*rep_mod_method*" - does not exist")
  end
  return data_mod,extr_vals,idcs
end

"""
    simple_extr_val_sel(data::ClustData,
                        extreme_value_descr::SimpleExtremeValueDescr;
                        rep_mod_method::String="feasibility")
Wrapper function for only one simple extreme value.
Selects simple extreme values and returns modified data, extreme values, and the corresponding indices.
"""
function simple_extr_val_sel(data::ClustData,
                             extr_value_descr::SimpleExtremeValueDescr;
                             rep_mod_method::String="feasibility"
                             )
  return simple_extr_val_sel(data,[extr_value_descr];rep_mod_method=rep_mod_method)
end

"""
    simple_extr_val_ident(data::ClustData,
                               extreme_value_descr_ar::Array{SimpleExtremeValueDescr,1})

Identifies multiple simple extreme values from the data and returns array of column indices of extreme value within data
- data_type: any attribute from the attributes contained within *data*
- extremum: "min" or "max"
- peak_def: "absolute" or "integral"
"""
function simple_extr_val_ident(data::ClustData,
                               extreme_value_descr_ar::Array{SimpleExtremeValueDescr,1})
  idcs = Array{Int,1}()
  # for each desired extreme value description, finds index of that extreme value within data
  for i=1:length(extreme_value_descr_ar)
    append!(idcs,simple_extr_val_ident(data,extreme_value_descr_ar[i]))
  end
  return idcs
end

"""
    simple_extr_val_ident(data::ClustData,
                               extreme_value_descr::SimpleExtremeValueDescr)

Wrapper function for only one simple extreme value:
identifies a single simple extreme value from the data and returns column index of extreme value
- `data_type`: any attribute from the attributes contained within *data*
- `extremum`: "min" or "max"
- `peak_def`: "absolute" or "integral"
- `consecutive_periods`: number of consecutive_periods combined to analyze
"""
function simple_extr_val_ident(data::ClustData,
                               extreme_value_descr::SimpleExtremeValueDescr)
  return simple_extr_val_ident(data, extreme_value_descr.data_type; extremum=extreme_value_descr.extremum, peak_def=extreme_value_descr.peak_def, consecutive_periods=extreme_value_descr.consecutive_periods)
end

"""
    simple_extr_val_ident(clust_data::ClustData,
                               data_type::String;
                               extremum::String="max",
                               peak_def::String="absolute",
                               consecutive_periods::Int=1)

Identifies a single simple extreme period from the data and returns column index of extreme period
- `data_type`: any attribute from the attributes contained within *data*
- `extremum`: "min" or "max"
- `peak_def`: "absolute" or "integral"
- `consecutive_periods`: The number of consecutive periods that are summed to identify a maximum or minimum. A rolling approach is used: E.g. for a value of `consecutive_periods`=2: 1) 1st & 2nd periods summed, 2) 2nd & 3rd period summed, 3) 3rd & 4th ... The min/max of the 1), 2), 3)... is determined and the two periods indices, where the  min/max were identified, are returned
"""
function simple_extr_val_ident(clust_data::ClustData,
                               data_type::String;
                               extremum::String="max",
                               peak_def::String="absolute",
                               consecutive_periods::Int=1)
  data=clust_data.data[data_type]
  delta_period=consecutive_periods-1
  # set data to be compared
  if peak_def=="absolute" && consecutive_periods==1
      data_eval = data
  elseif peak_def=="integral"
    # The number of consecutive_periods is substracted by one as k:k+period
    data_eval=zeros(1,(size(data,2)-delta_period))
    for k in 1:(size(data,2)-delta_period)
      data_eval[1,k] = sum(data[:,k:(k+delta_period)])
    end
  else
    error("peak_def - "*peak_def*" and consecutive_periods $consecutive_periods - not defined")
  end
  # find minimum or maximum index. Second argument returns cartesian indices, second argument of that is the column (period) index
  if extremum=="max"
    idx_k = findmax(data_eval)[2][2]
  elseif extremum=="min"
    idx_k = findmin(data_eval)[2][2]
  else
    error("extremum - "*extremum*" - not defined")
  end
  idx=collect(idx_k:(idx_k+delta_period))
  return idx
end

"""
    input_data_modification(data::ClustData,
                                 extr_val_idcs::Array{Int,1})

Returns ClustData structs with extreme vals and with remaining input data [data-extreme_vals].
Gives extreme vals the weight that they had in data.
This function is needed for the append method for representation modification
! the k-ids have to be monoton increasing - don't modify clustered data !
"""
function input_data_modification(data::ClustData,
                                 extr_val_idcs::Array{Int,1})
  unique_extr_val_idcs = unique(extr_val_idcs)
  K_dn = data.K- length(unique_extr_val_idcs)
  data_dn=Dict{String,Array}()
  index=setdiff(1:data.K,extr_val_idcs)
  for dt in keys(data.data)
    data_dn[dt] = data.data[dt][:,index] #take all columns but the ones that are extreme vals. If index occurs multiple times, setdiff only treats it as one.
  end
  weights_dn = data.weights[index]
  #take all columns but the ones that are extreme vals
  deltas_dn= data.delta_t[:,index]
  #deepcopy to change k_ids
  k_ids_dn=deepcopy(data.k_ids)
  #check for uniqueness and right sorting (however just those one representing)
  k_ids_check=k_ids_dn[findall(k_ids_dn.!=0)]
  allunique(k_ids_check) || error("the provided clust_data.k_ids are not unique - The clust_data is probably the result of a clustering already.")
  sort(k_ids_check)==k_ids_check || error("the provided clust_data.k_ids are not monoton increasing - The clust_data is probably the result of a clustering already.")
  #get all k-ids that are represented within this clust-data
  k_ids_dn_data=k_ids_dn[findall(data.k_ids.!=0)]
  for k in sort(extr_val_idcs)
    #reduce the following k_ids by one for all of the following k-ids (the deleted column will reduce the following column-indices by one for each deleted column)
    k_ids_dn_data[k:end].-=1
    #set this k_id to zero, as it corresponding column is being removed from the data
    k_ids_dn_data[k]=0
  end
  #just modify the k_ids that are also represented within this clust-data (don't reduce 0 to -1...)
  k_ids_dn[findall(data.k_ids.!=0)]=k_ids_dn_data
  #return the new Clust Data
  return ClustData(data.region,data.years,K_dn,data.T,data_dn,weights_dn,k_ids_dn;delta_t=deltas_dn,mean=data.mean,sdv=data.sdv)
end

"""
    input_data_modification(data::ClustData,extr_val_idcs::Int)

wrapper function for a single extreme val.
returns ClustData structs with extreme vals and with remaining input data [data-extreme_vals].
Gives extreme vals the weight that they had in data.
"""
function input_data_modification(data::ClustData,extr_val_idcs::Int)
  return input_data_modification(data,[extr_val_idcs])
end

"""
    extreme_val_output(data::ClustData,
                            extr_val_idcs::Array{Int,1};
                            rep_mod_method="feasibility")

Takes indices as input and returns ClustData struct that contains the extreme vals from within data.
"""
function extreme_val_output(data::ClustData,
                            extr_val_idcs::Array{Int,1};
                            rep_mod_method="feasibility")
  unique_extr_val_idcs = unique(extr_val_idcs)
  K_ed = length(unique_extr_val_idcs)
  data_ed=Dict{String,Array}()
  for dt in keys(data.data)
    data_ed[dt] = data.data[dt][:,unique_extr_val_idcs]
  end
  weights_ed=[]
  #initiate new k-ids-ed that don't represent any original time-period
  k_ids_ed=zeros(Int,size(data.k_ids))
  if rep_mod_method == "feasibility"
    weights_ed = zeros(length(unique_extr_val_idcs))
    #no representation is done of the original time-period, it's just for feasibility
  elseif rep_mod_method == "append"
    weights_ed = data.weights[unique_extr_val_idcs]
    # if original time series period isn't represented by any extreme period it has value 0
    # get all the indices that acutally represent original time-series
    index_k_ids_data=findall(data.k_ids.!=0)
    # each original time series period which is represented recieves the number of it's extreme period in this extreme value output
    k_ids_ed_data=zeros(size(index_k_ids_data))
    k_ids_ed_data[unique_extr_val_idcs]=collect(1:K_ed)
    # assign it to the full original time-series
    k_ids_ed[index_k_ids_data]=k_ids_ed_data
  else
    error("rep_mod_method - "*rep_mod_method*" - does not exist")
  end
  delta_t_ed=data.delta_t[:,unique_extr_val_idcs]
  extr_vals = ClustData(data.region,data.years,K_ed,data.T,data_ed,weights_ed,k_ids_ed;delta_t=delta_t_ed,mean=data.mean,sdv=data.sdv)
  return extr_vals
end

"""
   extreme_val_output(data::ClustData,extr_val_idcs::Array{Int,1};rep_mod_method="feasibility")

wrapper function for a single extreme val.
Takes indices as input and returns ClustData struct that contains the extreme vals from within data.
"""
function extreme_val_output(data::ClustData,
                            extr_val_idcs::Int;
                            rep_mod_method="feasibility")
  return extreme_val_output(data,[extr_val_idcs];rep_mod_method=rep_mod_method)
end

"""
    representation_modification(extr_vals::ClustData,clust_data::ClustData)

Merges the clustered data and extreme vals into one ClustData struct. Weights are chosen according to the rep_mod_method
"""
function representation_modification(extr_vals::ClustData,
                                     clust_data::ClustData)
  K_mod = clust_data.K + extr_vals.K
  data_mod=Dict{String,Array}()
  for dt in keys(clust_data.data)
    data_mod[dt] = [clust_data.data[dt] extr_vals.data[dt]]
  end
  weights_mod = [clust_data.weights; extr_vals.weights]
  # Add extra columns to delta_t
  delta_t_mod = [clust_data.delta_t extr_vals.delta_t]
  # originial time series periods are represented by periods in clust_data
  k_ids_mod=deepcopy(clust_data.k_ids)
  # if this particular original time series period is though represented in the extreme values, the new period number of the extreme value (clust_data.K+old number) is assigned to this original time series period - in case of feasibility they are all zero and nothing is changed
  k_ids_mod[findall(extr_vals.k_ids.!=0)]=extr_vals.k_ids[findall(extr_vals.k_ids.!=0)].+clust_data.K
  return ClustData(clust_data.region,clust_data.years,K_mod,clust_data.T,data_mod,weights_mod,k_ids_mod;delta_t=delta_t_mod,mean=clust_data.mean,sdv=clust_data.sdv)
end

"""
    representation_modification(extr_vals::ClustData,clust_data::ClustData)

Merges the clustered data and extreme vals into one ClustData struct. Weights are chosen according to the rep_mod_method
"""
function representation_modification(extr_vals_array::Array{ClustData,1},
                                     clust_data::ClustData,
                                     )
    for extr_vals in extr_vals_array
      clust_data=representation_modification(extr_vals,clust_data)
    end
    return clust_data
end

"""
    representation_modification(full_data::ClustData,clust_data::ClustData,extr_val_idcs::Array{Int,1};rep_mod_method::String="feasibility")

Merges the clustered data and extreme vals into one ClustData struct. Weights are chosen according to the rep_mod_method
"""
function representation_modification(full_data::ClustData,
                                     clust_data::ClustData,
                                     extr_val_idcs::Array{Int,1};
                                     rep_mod_method::String="feasibility")
  extr_vals = extreme_val_output(full_data,extr_val_idcs;rep_mod_method=rep_mod_method)
  return representation_modification(extr_vals,clust_data;rep_mod_method=rep_mod_method)
end

"""
    representation_modification(full_data::ClustData,clust_data::ClustData,extr_val_idcs::Int;rep_mod_method::String="feasibility")

wrapper function for a single extreme val.
Merges the clustered data and extreme vals into one ClustData struct. Weights are chosen according to the rep_mod_method
"""
function representation_modification(full_data::ClustData,
                                     clust_data::ClustData,
                                     extr_val_idcs::Int;
                                     rep_mod_method::String="feasibility")
  return representation_modification(full_data,clust_data,[extr_val_idcs];rep_mod_method=rep_mod_method)
end


text:
using JuMP, Gurobi, Graphs, JLD, DataFrames
include("DataImport.jl")
include("TypeInfo.jl")
include("contin_data.jl")
include("load24_data.jl")
include("BFsub.jl")
include("DelayMatrix.jl")
#### Optimization Parameters Control Area ####
gap = 1e-2
delayMode = 1
##############################################
########## Simulation Case Setting ##########
N_line = 14
N_time = 24
N_storage = 3
N_bus = 15
N_scen_normal = 1
N_scen_contin = 1
N_scen = N_scen_normal + N_scen_contin
N_event = 1
##############################################
filename_Node = "Node.csv"
filename_Generator = "Generator.csv"
filename_Line = "Line.csv"
filename_Storage = "Storage1MWhx3.csv"
filename_Station = "Station.csv"
filename_SMP = "DataMiner-Export_2018-01-10-121427.csv"
filename_Load24 = "load24data-PJM.csv"

load24 = load24_data(filename_Load24)
# N_days = size(load24,1)

contin = contin_data(N_line,N_time,N_scen_normal,N_scen_contin,N_event)
scen_data = ScenarioData(contin,load24[1:N_scen_normal,:])

status,obj_value,simulation_time,ESlocation,SoC,u,z,v,a,pg,qg,IC,OC,EC,OCinter,OCgen,OCloadshed,OCstorage,x,sigma_d,sigma_l = BFsub(scen_data,gap,delayMode,filename_Node, filename_Generator, filename_Line, filename_Storage, filename_Station, filename_SMP);

@save "result.jld"

ext: jl


text:
"""
    X = project_sequences(msa::AbstractMultipleSequenceAlignment; fracvar::Real = 0.9)

Perform a classical multidimensional scaling analysis to project the sequences in `msa` to a space
in which pairwise distances approximately reproduce `100 - percentsimilarity(seq1, seq2)`.
The dimensionality is chosen to reconstruction `fracvar` of the variance.
"""
function project_sequences(msa::AbstractMultipleSequenceAlignment; fracvar::Real = 0.9)
    sim = percentsimilarity(msa)
    D = 100 .- Matrix(sim)
    f = fit(MDS, D; distances=true)
    # Capture sufficient variance
    cl = cumsum(f.λ)
    cl ./= cl[end]
    nd = findfirst(>=(fracvar), cl)
    X = transform(f)
    return X[1:nd, :]
end

const reduced_code = ReducedAlphabet("(AILMV)(NQST)(RHK)(DE)(FWY)CGP")

function _entropy(v)
    count = countitems(v)
    e = 0.0
    for (_, n) in count
        p = n / length(v)
        e += - p * log(p)
    end
    return e
end

"""
    columnwise_entropy(msa::AbstractMultipleSequenceAlignment, aacode = ReducedAlphabet("(AILMV)(NQST)(RHK)(DE)(FWY)CGP")

Compute the entropy of each column in an MSA. Low entropy indicates high conservation.

Unmatched entries (`'-'` residues) contribute to the entropy calculation as if they were an ordinary residue.
"""
function columnwise_entropy(msa::AbstractMultipleSequenceAlignment, aacode=reduced_code)
    resnum = map(r -> aacode[r], getresidues(msa))
    return map(_entropy, eachcol(resnum))
end

"""
    residue_centroid(r::PDBResidue)

Compute the mean position of all atoms in `r` excluding hydrogen.

Residue centers may yield a more reliable measure of "comparable residues" than the α-carbons (see `CAmatrix`) because
they incorporate the orientation of the side chain relative to the overall fold.

See also [`residue_centroid_matrix`](@ref).
"""
residue_centroid(r::PDBResidue) = mean(atom -> atom.coordinates, Iterators.filter(atom -> atom.element != "H", r.atoms))

"""
    residue_centroid_matrix(seq)

Return a matrix of all residue centroids as columns. See also [`residue_centroid`](@ref).
"""
residue_centroid_matrix(seq) = reduce(hcat, map(residue_centroid, seq))

"""
    align(fixedpos::AbstractMatrix{Float64}, moving::AbstractVector{PDBResidue}, sm::SequenceMapping)
    align(fixed::AbstractVector{PDBResidue}, moving::AbstractVector{PDBResidue}, sm::SequenceMapping)

Return a rotated and shifted version of `moving` so that the α-carbons of residues `moving[sm]` have least
mean square error deviation from positions `fixedpos` or those of residues `fixed`.
"""
function align(fixedpos::AbstractMatrix{Float64}, moving::AbstractVector{PDBResidue}, sm::SequenceMapping; seqname=nothing)
    length(sm) == size(fixedpos, 2) || throw(DimensionMismatch("reference has $(size(fixedpos, 2)) positions, but `sm` has $(length(sm))"))
    movres = moving[sm]
    keep = map(!isnothing, movres)
    if !all(keep)
        @warn "missing $(length(keep)-sum(keep)) anchors for sequence $seqname"
    end
    movpos = residue_centroid_matrix(movres[keep])
    fixedpos = fixedpos[:,keep]
    refmean = mean(fixedpos; dims=2)
    movmean = mean(movpos; dims=2)
    R = kabsch((fixedpos .- refmean)', (movpos .- movmean)')
    return change_coordinates(moving, (coordinatesmatrix(moving) .- movmean') * R .+ refmean')
end
align(fixed::AbstractVector{PDBResidue}, moving::AbstractVector{PDBResidue}, sm::SequenceMapping; kwargs...) =
    align(residue_centroid_matrix(fixed), moving, sm; kwargs...)

function mapclosest(mapto::AbstractVector{PDBResidue}, mapfrom::AbstractVector{PDBResidue})
    refcenters = residue_centroid_matrix(mapto)
    seqcenters = residue_centroid_matrix(mapfrom)
    D = pairwise(Euclidean(), refcenters, seqcenters; dims=2)
    assignment, _ = hungarian(D)
    return collect(zip(assignment, [j == 0 ? convert(eltype(D), NaN) : D[i, j] for (i, j) in enumerate(assignment)]))
    # d, m = findmin(D; dims=1)
    # return [mi[1] => di for (mi, di) in zip(vec(m), vec(d))]
end

"""
    chargelocations(pdb::AbstractVector{PDBResidue}; include_his::Bool=false)

Return a list of potential charge locations in the protein structure. Each is a tuple `(position, residueindex, AAname)`.
The positions are those of `N` (in positively-charged residues like Arg & Lys) or `O` (in negatively-charged residues like
Asp and Glu). N- and C-termini are not included in the list. While each residue will carry a net total charge of ±1,
the location of each potential charge will be listed (1 for Lys, 2 each for Arg, Asp, and Glu).

By default, histidine is not considered charged, but you can include it by setting `include_His=true`.
"""
function chargelocations(pdb::AbstractVector{PDBResidue}; include_His::Bool=false)
    argpat(atom) = startswith(atom.atom, "NH")  # arginine
    lyspat(atom) = atom.atom == "NZ"            # lysine
    asppat(atom) = startswith(atom.atom, "OD")  # aspartate
    glupat(atom) = startswith(atom.atom, "OE")  # glutamate
    hispat(atom) = startswith(atom.atom, r"N[DE]")  # histidine

    locs = Tuple{Coordinates,Int,String}[]
    for (i, res) in pairs(pdb)
        name = res.id.name
        if name == "ARG"
            idx = findfirst(argpat, res.atoms)
            push!(locs, (res.atoms[idx].coordinates, i, name))
            idx = findnext(argpat, res.atoms, idx+1)
            push!(locs, (res.atoms[idx].coordinates, i, name))
        elseif name == "LYS"
            idx = findfirst(lyspat, res.atoms)
            push!(locs, (res.atoms[idx].coordinates, i, name))
        elseif name == "ASP"
            idx = findfirst(asppat, res.atoms)
            push!(locs, (res.atoms[idx].coordinates, i, name))
            idx = findnext(asppat, res.atoms, idx+1)
            push!(locs, (res.atoms[idx].coordinates, i, name))
        elseif name == "GLU"
            idx = findfirst(glupat, res.atoms)
            push!(locs, (res.atoms[idx].coordinates, i, name))
            idx = findnext(glupat, res.atoms, idx+1)
            push!(locs, (res.atoms[idx].coordinates, i, name))
        elseif include_His && name == "HIS"
            idx = findfirst(hispat, res.atoms)
            push!(locs, (res.atoms[idx].coordinates, i, name))
            idx = findnext(hispat, res.atoms, idx+1)
            push!(locs, (res.atoms[idx].coordinates, i, name))
        end
    end
    return locs
end

positive_locations(chargelocs::AbstractVector{Tuple{Coordinates,Int,String}}) = [loc[1] for loc in chargelocs if loc[3] ∈ ("ARG", "LYS", "HIS")]
negative_locations(chargelocs::AbstractVector{Tuple{Coordinates,Int,String}}) = [loc[1] for loc in chargelocs if loc[3] ∈ ("ASP", "GLU")]


text:
config = Dict{Any, Any}(

    "DRIVER" => Dict(
        "casename"           => "Sandbox",
        "caseroot"           => joinpath(@__DIR__, "Sandbox", "caseroot"),
        "caserun"            => joinpath(@__DIR__, "Sandbox", "caserun"),
        "archive_root"       => joinpath(@__DIR__, "Sandbox", "archive"),
    ),

    "MODEL_MISC" => Dict(
        "timetype"               => "DateTimeNoLeap",
        "init_file"              => "BLANK_PROFILE",#nothing,#joinpath(@__DIR__, "ocn_init.nc"),
        "rpointer_file"          => "rpointer.hoom",
        "monthly_record"         => ["{ESSENTIAL}",],
        "daily_record"         => ["{ESSENTIAL}", "Ks_H_U", "Ks_H_V"],
        "enable_archive"         => true,
    ),

    "MODEL_CORE" => Dict(
        #"domain_file"                  => joinpath(@__DIR__, "CESM_domains", "domain.ocn.gx1v6.090206.nc"),
        #"cdata_file"                   => joinpath(@__DIR__, "forcing_gx1v6.nc"),
        "domain_file"                  => joinpath(@__DIR__, "domain.ocn_aqua.fv4x5_gx3v7.091218.nc"),
        "cdata_file"                   => joinpath(@__DIR__, "forcing.nc"),

        "cdata_beg_time"               => DateTimeNoLeap(1, 1, 1, 0, 0, 0),
        "cdata_end_time"               => DateTimeNoLeap(2, 1, 1, 0, 0, 0),
        "cdata_align_time"             => DateTimeNoLeap(1, 1, 1, 0, 0, 0),

        "z_w"                          => collect(Float64, 0:-10:-350),

        "substeps"                     => 8,
        "MLD_scheme"                   => "static",
        "Qflx"                         => "off",
        "Qflx_finding"                 => "off",
        "convective_adjustment"        => "off",
        "advection_scheme"             => "ekman_CO2012",

        "weak_restoring"               => "off",
        "τwk_TEMP"                     => 86400.0 * 365,
        "τwk_SALT"                     => 86400.0 * 365,

        "topo_file" => "",

        "τ_frz"                        => 3600.0,
        "Ekman_layers"      => 5,
        "Returnflow_layers" => 25,
    ),

)

max_stars_repo_name: meteorologytoday/IOM

max_issues_repo_name: meteorologytoday/IOM

max_forks_repo_name: meteorologytoday/IOM


text:
using Cubature
using Distributions
using PyCall
@pyimport matplotlib as mpl
mpl.rc("text", usetex=true)
mpl.use("Agg")
@pyimport matplotlib.ticker as mti
using PyPlot

Cg(x::Vector{Float64}, θ::Union{Int, Float64}) = exp(-((-log(x[1]))^θ+(-log(x[2]))^θ)^(1/θ))

function Ccl(x::Vector{Float64}, θ::Union{Int, Float64})
  if θ > 0
    return (x[1]^(-θ)+x[2]^(-θ)-1)^(-1/θ)
  else
    return (maximum([x[1]^(-θ)+x[2]^(-θ)-1, 0]))^(-1/θ)
  end
end

Cf(x::Vector{Float64}, θ::Union{Int, Float64}) =
  -log(1+((exp(-θ*x[1])-1)*(exp(-θ*x[2])-1))/(exp(-θ)-1))/θ

Camh(x::Vector{Float64}, θ::Union{Int, Float64}) = x[1]*x[2]/(1-θ*(1-x[1])*(1-x[2]))

function getcov(copula::Function, d, θa::Vector{Float64})
  covs = Float64[]
  a = 0.000000001
  b = 0.999999999
  f(y) = pdf.(d, quantile.(d, y))
  for θ in θa
    intf(x) = (copula(x, θ)-x[1]*x[2])/(f(x[1])*f(x[2]))
    push!(covs, hcubature(intf, [a,a],[b,b])[1]/var(d))
  end
  covs
end

function getρ(copula::Function, θa::Vector{Float64})
  covs = Float64[]
  for θ in θa
    push!(covs, 12*hcubature(x-> copula(x, θ), [0,0],[1,1])[1]-3)
  end
  covs
end

θ = [i for i in -0.5:0.1:1.]
function plotcor(copula::Function, θ::Vector{Float64})
  mpl.rc("font", family="serif", size = 7)
  fig, ax = subplots(figsize = (2.5, 2.))
  ax[:plot](θ, getcov(copula, Normal(0,1), θ), label= "normal dist", color = "r", linewidth = 0.5)
  ax[:plot](θ, getcov(copula, LogNormal(1,1), θ), label= "lognormal(1,1)", color = "blue", linewidth = 0.5)
  ax[:plot](θ, getcov(copula, Arcsine(0,1), θ), label= "arcsine(0,1)", color = "green", linewidth = 0.5)
  ax[:plot](θ, getcov(copula, TDist(5), θ), label= "tdist(5)", color = "grey", linewidth = 0.5)
  ax[:plot](θ, getcov(copula, Weibull(1,1), θ), label= "weibull(1,1)", color = "yellow", linewidth = 0.5)
  ax[:plot](θ, getρ(copula, θ), label = "uniform dist", color = "black", linewidth = 0.5)
  ax[:legend](fontsize = 4.5, loc = 4, ncol = 1)
  subplots_adjust(bottom = 0.16, top=0.92, left = 0.18, right = 0.92)
  PyPlot.ylabel("Pearson correlation", labelpad = -1.5)
  PyPlot.xlabel("\$ \\theta \$", labelpad = 0.6)
  fig[:savefig]("pics/cor"*string(copula)*".eps")
end

plotcor(Cg, vcat([1.01], [i for i in 1.2:0.2:10.]))
plotcor(Cf, [i for i in -10.:0.2:10.])
plotcor(Ccl, vcat([i for i in -1.:0.2:-0.1], [i for i in 0.1:0.2:10.]))
plotcor(Camh, [i for i in -1.:0.05:1.])


text:
using OrdinaryDiffEq, Test

function time_derivative(du,u,p,t)
  du[1] = -t
end
function time_derivative_analytic(u0,p,t)
  u0 .- t.^2 ./ 2
end

ff_time_derivative = ODEFunction(time_derivative,
                                 analytic=time_derivative_analytic)
u0 = [1.0]
tspan = (0.0,1.0)
prob = ODEProblem(ff_time_derivative,u0,tspan)

sol = solve(prob,Rosenbrock32(),reltol=1e-9,abstol=1e-9)
@test sol.errors[:final] < 1e-5
@show Rosenbrock23()
sol = solve(prob,Rosenbrock23(),reltol=1e-9,abstol=1e-9)
@test sol.errors[:final] < 1e-10
@show Rodas4
sol = solve(prob,Rodas4(),reltol=1e-9,abstol=1e-9)
@test sol.errors[:final] < 1e-10
@show Rodas5
sol = solve(prob,Rodas5(),reltol=1e-9,abstol=1e-9)
@test sol.errors[:final] < 1e-10
@show Veldd4
sol = solve(prob,Veldd4(),reltol=1e-9,abstol=1e-9)
@test sol.errors[:final] < 1e-10
@show KenCarp3
sol = solve(prob,KenCarp3(),reltol=1e-12,abstol=1e-12)
@test length(sol) > 2
@test sol.errors[:final] < 1e-10
@show KenCarp4
sol = solve(prob,KenCarp4(),reltol=1e-12,abstol=1e-12)
@test length(sol) > 2
@test sol.errors[:final] < 1e-10
@show KenCarp47
sol = solve(prob,KenCarp47(),reltol=1e-12,abstol=1e-12)
@test length(sol) > 2
@test sol.errors[:final] < 1e-10
@show KenCarp5
sol = solve(prob,KenCarp5(),reltol=1e-12,abstol=1e-12)
@test length(sol) > 2
@test sol.errors[:final] < 1e-10
@show KenCarp58
sol = solve(prob,KenCarp5(),reltol=1e-12,abstol=1e-12)
@test length(sol) > 2
@test sol.errors[:final] < 1e-10
@show TRBDF2
sol = solve(prob,TRBDF2(),reltol=1e-9,abstol=1e-9)
@test sol.errors[:final] < 1e-10
sol = solve(prob,ImplicitEuler(),dt=1/10)
@test sol.errors[:final] < 1e-1
sol = solve(prob,Trapezoid(),dt=1/10)
@test sol.errors[:final] < 1e-12
sol = solve(prob,Euler(),dt=1/100)
@test sol.errors[:final] < 6e-3

const CACHE_TEST_ALGS = [Euler(),Midpoint(),RK4(),SSPRK22(),SSPRK33(),SSPRK53(),
  SSPRK63(),SSPRK73(),SSPRK83(),SSPRK432(),SSPRK932(),SSPRK54(),SSPRK104(),CarpenterKennedy2N54(),
  BS3(),BS5(),DP5(),DP8(),Feagin10(),Feagin12(),Feagin14(),TanYam7(),
  Tsit5(),TsitPap8(),Vern6(),Vern7(),Vern8(),Vern9(),OwrenZen3(),OwrenZen4(),OwrenZen5()]

for alg in CACHE_TEST_ALGS
  sol = solve(prob,alg,dt=1/10)
  if !(typeof(alg) <: Euler)
    @test sol.errors[:final] < 4e-14
  end
end


text:
@info "adaptive_reweighting_tests"
using Flux, OptimizationFlux
using DiffEqFlux
using ModelingToolkit
using DiffEqBase
using Test, NeuralPDE
using Optimization
using OptimizationOptimJL
using Quadrature,Cubature, Cuba
using QuasiMonteCarlo
using SciMLBase
import ModelingToolkit: Interval, infimum, supremum
using DomainSets
using Random
#using Plots
@info "Starting Soon!"

nonadaptive_loss = NeuralPDE.NonAdaptiveLoss(pde_loss_weights=1, bc_loss_weights=1)
gradnormadaptive_loss = NeuralPDE.GradientScaleAdaptiveLoss(100, pde_loss_weights=1e3, bc_loss_weights=1)
adaptive_loss = NeuralPDE.MiniMaxAdaptiveLoss(100; pde_loss_weights=1, bc_loss_weights=1)
adaptive_losses = [nonadaptive_loss, gradnormadaptive_loss,adaptive_loss]
maxiters=4000
seed=60

## 2D Poisson equation
function test_2d_poisson_equation_adaptive_loss(adaptive_loss; seed=60, maxiters=4000)
    Random.seed!(seed)
    hid = 40
    chain_ = FastChain(FastDense(2,hid,Flux.σ),FastDense(hid,hid,Flux.σ),FastDense(hid,1))
    strategy_ =  NeuralPDE.StochasticTraining(256)
    @info "adaptive reweighting test outdir:, maxiters: $(maxiters), 2D Poisson equation, adaptive_loss: $(nameof(typeof(adaptive_loss))) "
    @parameters x y
    @variables u(..)
    Dxx = Differential(x)^2
    Dyy = Differential(y)^2

    # 2D PDE
    eq  = Dxx(u(x,y)) + Dyy(u(x,y)) ~ -sin(pi*x)*sin(pi*y)

    # Initial and boundary conditions
    bcs = [u(0,y) ~ 0.0, u(1,y) ~ -sin(pi*1)*sin(pi*y),
           u(x,0) ~ 0.0, u(x,1) ~ -sin(pi*x)*sin(pi*1)]
    # Space and time domains
    domains = [x ∈ Interval(0.0,1.0),
               y ∈ Interval(0.0,1.0)]

    initθ = Float64.(DiffEqFlux.initial_params(chain_))
    iteration = [0]
    discretization = NeuralPDE.PhysicsInformedNN(chain_,
                                                 strategy_;
                                                 init_params = initθ,
                                                 adaptive_loss = adaptive_loss,
                                                 logger = nothing,
                                                 iteration=iteration)

    @named pde_system = PDESystem(eq,bcs,domains,[x,y],[u(x, y)])
    prob = NeuralPDE.discretize(pde_system,discretization)
    phi = discretization.phi
    sym_prob = NeuralPDE.symbolic_discretize(pde_system,discretization)

    xs,ys = [infimum(d.domain):0.01:supremum(d.domain) for d in domains]
    analytic_sol_func(x,y) = (sin(pi*x)*sin(pi*y))/(2pi^2)
    u_real = reshape([analytic_sol_func(x,y) for x in xs for y in ys], (length(xs),length(ys)))

    callback = function (p,l)
        iteration[1] += 1
        if iteration[1] % 100 == 0
            @info "Current loss is: $l, iteration is $(iteration[1])"
        end
        return false
    end
    res = Optimization.solve(prob, ADAM(0.03); maxiters=maxiters, callback =callback)

    u_predict = reshape([first(phi([x,y],res.minimizer)) for x in xs for y in ys],(length(xs),length(ys)))
    diff_u = abs.(u_predict .- u_real)
    total_diff = sum(diff_u)
    total_u = sum(abs.(u_real))
    total_diff_rel = total_diff / total_u

    #p1 = plot(xs, ys, u_real, linetype=:contourf,title = "analytic");
    #p2 = plot(xs, ys, u_predict, linetype=:contourf,title = "predict");
    #p3 = plot(xs, ys, diff_u,linetype=:contourf,title = "error");
    #(plot=plot(p1,p2,p3), error=total_diff, total_diff_rel=total_diff_rel)
    (error=total_diff, total_diff_rel=total_diff_rel)
end

@info "testing that the adaptive loss methods roughly succeed"
test_2d_poisson_equation_adaptive_loss_no_logs_run_seediters(adaptive_loss) = test_2d_poisson_equation_adaptive_loss(adaptive_loss; seed=seed, maxiters=maxiters)
error_results_no_logs = map(test_2d_poisson_equation_adaptive_loss_no_logs_run_seediters, adaptive_losses)

# accuracy tests
@show error_results_no_logs[1][:total_diff_rel]
@show error_results_no_logs[2][:total_diff_rel]
@show error_results_no_logs[3][:total_diff_rel]
# accuracy tests, these work for this specific seed but might not for others
# note that this doesn't test that the adaptive losses are outperforming the nonadaptive loss, which is not guaranteed, and seed/arch/hyperparam/pde etc dependent
@test error_results_no_logs[1][:total_diff_rel] < 0.4
@test error_results_no_logs[2][:total_diff_rel] < 0.4
@test error_results_no_logs[3][:total_diff_rel] < 0.4

#plots_diffs[1][:plot]
#plots_diffs[2][:plot]
#plots_diffs[3][:plot]


text:
export Berlage,
Ormsby,
Ricker
include("Berlage.jl")
include("Ormsby.jl")
include("Ricker.jl")

max_stars_repo_name: fercarozzi/myseismicjulia

max_issues_repo_name: fercarozzi/myseismicjulia

max_forks_repo_name: fercarozzi/myseismicjulia


text:
# This script will generate a `userimg.jl` file.
# Append the generated `userimg.jl` file to your base/userimg.jl
# and build julia in order to add this package to your system image.

using SnoopCompile

compiles_csv = "compiles.csv"
SnoopCompile.@snoop compiles_csv begin
    include(Pkg.dir("BusinessDays", "test", "runtests.jl"))
end

blacklist = ["Base.Test", "yyyy-mm-dd", "CustomCalendar", "TestCalendar", "mm"]

data = SnoopCompile.read(compiles_csv)
pc = SnoopCompile.format_userimg(reverse!(data[2]), blacklist=blacklist)
unshift!(pc, "using BusinessDays")
SnoopCompile.write("userimg.jl", pc)
rm(compiles_csv)


text:
# these graph layout functions are experimental and mostly for decorative purposes
# adapted from various places, modified for use with Luxor

#= Example of spring layout

using Luxor, LightGraphs, Combinatorics

g = smallgraph(:moebiuskantor)

# spring layout

x, y = layoutgraph(Array(adjacency_matrix(g)),
    densityconstant    = 20,
    maxiterations      = 50,
    initialtemperature = 10)

@draw begin
    pts = Point.(zip(x, y))
    circle.(pts, 15, :fill)
    for p in combinations(pts, 2)
        randomhue()
        arrow(p[1], p[2],  decorate = () -> ngon(O, 10, 3, 0, :fill))
    end
end

=#

"""
    layoutgraph(adj_matrix::Array{T,2}; kwargs...)       <--
    layoutgraph(adj_list::AbstractVector; kwargs...)     not yet implemented

Run one of two graph layout algorithms.

1: Layout the graph in adjacency matrix `adj_matrix`, using
a spring-based method. Returns a tuple of `([xcoords],
[ycoords])`, which are scaled according to the supplied
BoundingBox.

### Keyword arguments

Lay out the graph in `adjmatrix` using the spring/repulsion
model of Fruchterman and Reingold (1991). Returns a tuple of
`([xcoords], [ycoords])`.

## Arguments

`am` An adjacency matrix of some type. Non-zero of the
eltype of the matrix is used to determine if a link exists,
but currently no sense of magnitude

`densityconstant` Constant to adjust the density of
resulting layout. The default value is 2.0.

`maxiterations` how many iterations for applying the forces

`initialtemperature` the initial temperature controls
movement per iteration. The default value is 2.0. Each
iteration uses a lower temperature (TEMP =
initialtemperature / iter). The idea is that the
displacements of vertices are limited to some maximum
temperature value, and this decreases over time. As the
layout becomes better, the amount of adjustment becomes
smaller.

`locs_x` the starting x values, between -1 and 1. The
default values are randomly selected.

`locs_y` the starting y values, between -1 and 1. The
default values are randomly selected.

`boundingbox` The Luxor BoundingBox into which the graph's
coordinates will fit. Defaults to 500 × 500.

"""
layoutgraph(am::Array{T,2} where T; kwargs...) = layout_spring(am; kwargs...)

"""
    layoutgraph(adj_matrix::Array{T,2}; kwargs...)
    layoutgraph(adj_list::AbstractVector; kwargs...)      <--

2: Layout the graph in adjacency list `adj_list`, assuming
the graph is a tree layout with a single root node. Return a
list of Luxor points that can be used to lay out the nodes
of the tree.

### Keyword arguments

```nodesize=ones(size(adjlist, 1)))```

`nodesize` is a vector of values that define the area
occupied by each node. No idea how it works yet...

Return a list of Luxor points that can be used to lay out a
tree defined in adjacency list `adjlist`. """

layoutgraph(al::AbstractVector; kwargs...) = layout_tree(al;
kwargs...)

"""
    layout_spring(adjmatrix::Array{T,2} where T;
        densityconstant = 2.0,
        maxiterations = 100,
        initialtemperature = 2.0,
        boundingbox=BoundingBox(O - (250, 250), O + (250, 250)))

"""
function layout_spring(adjmatrix::Array{T,2} where T;
    densityconstant = 2.0,
    maxiterations = 100,
    initialtemperature = 2.0,
    locs_x = 2rand(size(adjmatrix, 1)) .- 1.0,
    locs_y = 2rand(size(adjmatrix, 1)) .- 1.0,
    boundingbox =  BoundingBox(O - (250, 250), O + (250, 250)))

    N = size(adjmatrix, 1)
    if N != size(adjmatrix, 2)
        error("Adjacency matrix must be square but is $(size(adjmatrix)).")
    end

    # randomize results arrays, normalized to between -1 and 1

    # the optimal distance bewteen normalized vertices
    K = densityconstant * sqrt(4.0 / N)

    # store forces and apply at end of iteration all at once
    force_x = zeros(N)
    force_y = zeros(N)

    @inbounds for iter in 1:maxiterations
        # Calculate forces
        for i in 1:N
            force_vec_x = 0.0
            force_vec_y = 0.0
            for j in 1:N
                i == j && continue
                d_x = locs_x[j] - locs_x[i]
                d_y = locs_y[j] - locs_y[i]
                d   = sqrt(d_x^2 + d_y^2)
                if adjmatrix[i,j] != zero(eltype(adjmatrix)) || adjmatrix[j,i] != zero(eltype(adjmatrix))
                    F_d = d^2 / K
                else
                    # Just repulsive
                    F_d = -K^2 / d
                end
                force_vec_x += F_d * d_x
                force_vec_y += F_d * d_y
            end
            force_x[i] = force_vec_x
            force_y[i] = force_vec_y
        end
        # Cool down
        TEMP = initialtemperature / iter
        # apply but limit to temperature
        for i = 1:N
            force_mag  = sqrt(force_x[i] ^ 2 + force_y[i] ^ 2)
            scale      = min(force_mag, TEMP)/force_mag
            locs_x[i] += force_x[i] * scale
            locs_y[i] += force_y[i] * scale
        end
    end

    # fit to bounding box argument
    min_x, max_x = boxbottomleft(boundingbox).x, boxbottomright(boundingbox).x
    min_y, max_y = boxtopleft(boundingbox).y, boxbottomleft(boundingbox).y

    return (
        rescale.(locs_x, extrema(locs_x)..., min_x, max_x),
        rescale.(locs_y, extrema(locs_y)..., min_y, max_y)
        )
end

"""
   layout_tree(adjlist)

`adj_list` must not be cyclic, otherwise welcome to stack overflow...
"""
function layout_tree(adj_list::AbstractVector)
    return _layout_tree(adj_list)
end


text:
@testset "Read a Table Schema descriptor" begin

    DESCRIPTOR_MIN = Dict("fields" => [
        Dict( "name" => "id" ),
        Dict( "name" => "height", "type" => "integer" )
    ])

    @testset "Minimal from dictionary" begin
        s = Schema(DESCRIPTOR_MIN)
        @test length(s.fields) == 2
        f1 = s.fields[1]
        @test f1.name == "id"
        @test f1.typed == "string"
        f2 = s.fields[2]
        @test f2.typed == "integer"
        @test !f2.constraints.required
    end

    @testset "Parsed from a JSON string" begin
        s = Schema(joinpath(dirname(@__FILE__), "../data/schema_valid_infer.json"))
        @test length(s.fields) == 2
        @test s.fields[1].name == "id"
        @test !s.fields[2].constraints.required
    end

    @testset "Full descriptor from JSON" begin
        s = Schema(joinpath(dirname(@__FILE__), "../data/schema_valid_full.json"))
        @test length(s.fields) == 15
        @test length(s.primary_key) == 4
    end

    @testset "Missing values and constraints" begin
        s = Schema(joinpath(dirname(@__FILE__), "../data/schema_valid_missing.json"))
        @test length(s.fields) == 5
        @test length(s.primary_key) == 1
        @test length(s.missing_values) == 3
        @test s.fields[1].constraints.required
        @test !(s.fields[1].constraints.unique)
    end

end

@testset "Field value casting" begin
    import TableSchema: cast_by_type, cast_value, cast_row, CastError

    DESCRIPTOR_MAX = Dict(
        "fields" => [
            Dict( "name" => "id", "type" => "string",
                "constraints" => Dict( "required" => true )
            ),
            Dict( "name" => "height", "type" => "number" ),
            Dict( "name" => "age", "type" => "integer" ),
            Dict( "name" => "name", "type" => "string" ),
            Dict( "name" => "occupation", "type" => "string" )
        ],
        "primaryKey" => [ "id" ],
        "foreignKeys" => [ Dict(
            "fields" => [ "name" ],
            "reference" => Dict(
                "resource" => "",
                "fields" => [ "id" ]
            )
        ) ],
        "missingValues" => ["", "-", "null"]
    )

    @testset "Cast functions" begin
        @test cast_by_type(Nothing, "any", "", Dict()) == Nothing

        @test cast_by_type(true, "boolean", "", Dict())
        @test cast_by_type("true", "boolean", "", Dict())
        @test isa(cast_by_type("", "boolean", "", Dict()), CastError)

        @test cast_by_type(123, "integer", "", Dict()) == 123
        @test cast_by_type("123", "integer", "", Dict()) == 123
        @test isa(cast_by_type("a", "integer", "", Dict()), CastError)

        @test cast_by_type(123, "number", "", Dict()) == 123.0
        @test cast_by_type("12.3", "number", "", Dict()) == 12.3
        @test isa(cast_by_type("a", "number", "", Dict()), CastError)

        @test cast_by_type("{\"a\":1}", "object", "", Dict()) == Dict("a"=>1)
        @test isa(cast_by_type("[\"a\",1]", "object", "", Dict()), CastError)

        @test cast_by_type(123, "string", "", Dict()) == "123"
        @test cast_by_type("123", "string", "", Dict()) == "123"
        @test isa(cast_by_type([1,2], "string", "", Dict()), CastError)
    end

    @testset "Cast row" begin
        s = Schema(DESCRIPTOR_MAX)
        source = ["string", "10.2", "1", "string", "string"]
        target = ["string", 10.2, 1, "string", "string"]
        @test cast_value(s.fields[1], "test", false) == "test"
        @test cast_row(s, source, false, false) == target
    end

    @testset "Cast row too short or long" begin
        s = Schema(DESCRIPTOR_MAX)
        source = ["string", "10.2", "1", "string"]
        @test_throws CastError cast_row(s, source, false, false)
        source = ["string", "10.2", "1", "string", "string", "string"]
        @test_throws CastError cast_row(s, source, false, false)
    end

    @testset "Cast row incorrect type" begin
        s = Schema(DESCRIPTOR_MAX)
        source = ["string", "notedecimal", "1", [1,2], "string"]
        @test_throws CastError cast_row(s, source, false, false)
        try
            cast_row(s, source, false, false)
        catch e
            if isa(e, CastError)
                @test length(e.errors) == 2
            else
                throw(e)
            end
        end
    end
end

max_stars_repo_name: loleg/tableschema-jl

max_issues_repo_name: loleg/tableschema-jl

max_forks_repo_name: loleg/tableschema-jl


text:
using Reproduce

tst = ItemCollection("/home/matthewmcleod/Documents/Masters/curiosity/MountainCar/data")

function load_results(ic, param)
    FileIO.load(joinpath(itm.folder_str, "results.jld2")))
end

max_stars_repo_name: MatthewMcLeod/curiosity

max_issues_repo_name: MatthewMcLeod/curiosity

max_forks_repo_name: MatthewMcLeod/curiosity


text:
module StreamRequest

import ..Layer, ..request
using ..IOExtras
using ..Messages
using ..Streams
import ..ConnectionPool
using ..MessageRequest
import ..@debug, ..DEBUG_LEVEL, ..printlncompact, ..sprintcompact

"""
    request(StreamLayer, ::IO, ::Request, body) -> HTTP.Response

Create a [`Stream`](@ref) to send a `Request` and `body` to an `IO`
stream and read the response.

Send the `Request` body in a background task and begins reading the response
immediately so that the transmission can be aborted if the `Response` status
indicates that the server does not wish to receive the message body.
[RFC7230 6.5](https://tools.ietf.org/html/rfc7230#section-6.5).
"""
abstract type StreamLayer{Next <: Layer} <: Layer{Next} end
export StreamLayer

function request(::Type{StreamLayer{Next}}, io::IO, request::Request, body;
                 response_stream=nothing,
                 iofunction=nothing,
                 verbose::Int=0,
                 kw...)::Response where Next

    verbose == 1 && printlncompact(request)

    response = request.response
    http = Stream(response, io)
    startwrite(http)

    if verbose == 2
        println(request)
        if iofunction === nothing && request.body === body_is_a_stream
            println("$(typeof(request)).body: $(sprintcompact(body))")
        end
    end

    if !isidempotent(request)
        # Wait for pipelined reads to complete
        # before sending non-idempotent request body.
        startread(io)
    end

    aborted = false
    write_error = nothing
    try

        @sync begin
            if iofunction === nothing
                @async try
                    writebody(http, request, body)
                catch e
                    write_error = e
                    @info("Error in @async writebody task. " *
                          "Server likely closed the connection unexpectedly. " *
                          "Only an issue if unable to read the response and this error gets re-thrown. ",
                          exception=(write_error, catch_backtrace()))
                    isopen(io) && close(io)
                end
                yield()
                startread(http)
                readbody(http, response, response_stream)
            else
                iofunction(http)
            end

            if isaborted(http)
                close(io)
                aborted = true
            end
        end

    catch e
        if write_error !== nothing
            throw(write_error)
        else
            rethrow(e)
        end
    end

    closewrite(http)
    closeread(http)

    verbose == 1 && printlncompact(response)
    verbose == 2 && println(response)

    return response
end

function writebody(http::Stream, req::Request, body)

    if req.body === body_is_a_stream
        writebodystream(http, req, body)
        closebody(http)
    else
        write(http, req.body)
    end

    req.txcount += 1

    if isidempotent(req)
        closewrite(http)
    else
        @debug 2 "🔒  $(req.method) non-idempotent, " *
                 "holding write lock: $(http.stream)"
        # "A user agent SHOULD NOT pipeline requests after a
        #  non-idempotent method, until the final response
        #  status code for that method has been received"
        # https://tools.ietf.org/html/rfc7230#section-6.3.2
    end
end

function writebodystream(http, req, body)
    for chunk in body
        writechunk(http, req, chunk)
    end
end

function writebodystream(http, req, body::IO)
    req.body = body_was_streamed
    write(http, body)
end

writechunk(http, req, body::IO) = writebodystream(http, req, body)
writechunk(http, req, body) = write(http, body)

function readbody(http::Stream, res::Response, response_stream)
    if response_stream == nothing
        res.body = read(http)
    else
        res.body = body_was_streamed
        write(response_stream, http)
        close(response_stream)
    end
end

end # module StreamRequest


text:
module IntervalBasedForest

using DecisionTree: apply_tree_proba, build_tree, Node, prune_tree
using Statistics, StableRNGs

"""
   `TimeSeriesForestClassifier(m, Xmatrix::Array, yplain::Array)`

TimeSeriesForestClassifier builds a forest of trees from the training set
`(Xmatrix, yplain)` using random intervals and summary features.
"""
function TimeSeriesForestClassifier(m, Xmatrix::Array, yplain::Array)
    n_trees = m.n_trees
    Interval_features, Intervals = InvFeatureGen(Xmatrix, n_trees, m.min_interval,
                                                 m.random_state)
    forest = Array{Node,1}()
    @inbounds for i in range(1, stop=n_trees)
        tree = build_tree(yplain, view(Interval_features, :, :, i),
                         m.n_subfeatures,
                         m.max_depth,
                         m.min_samples_leaf,
                         m.min_samples_split,
                         m.min_purity_increase)
        if m.post_prune
            tree = prune_tree(tree, m.merge_purity_threshold)
        end
        push!(forest, tree)
    end
    forest, Intervals
end

"""
   `predict_new(Xmatrix, forest, Intervals::Array, integers_seen)`

Find predictions for all cases in Xmatrix. Built on top of `apply_tree_proba`.
"""
function predict_new(Xmatrix::Array, forest, Intervals::Array, integers_seen)
    n_instance, s_length = size(Xmatrix)
    n_trees = length(forest)
    Interval_features = InvFeatureGen(Xmatrix, Intervals, n_trees)
    sum = zeros(n_instance, length(integers_seen))
    @inbounds for i=1:n_trees
        sum += apply_tree_proba(forest[i], view(Interval_features, :, :, i), integers_seen)
    end
    return sum/n_trees
end

"""
   `InvFeatureGen(Xmatrix::Array, n_trees::Int, min_interval::Int, random_state)`

InvFeatureGen selects random intervals, generating the mean, standard deviation
and slope of the random intervals resulting 3*√m features.
"""
function InvFeatureGen(Xmatrix::Array, n_trees::Int, min_interval::Int, random_state)
    n_instance, series_length = size(Xmatrix)
    n_intervals = floor(Int64, sqrt(series_length))
    Interval_features = Array{Float64,3}(undef, n_instance, 3*n_intervals, n_trees)
    Intervals = zeros(Int64, n_trees, n_intervals, 2)
    @inbounds for i in range(1, stop = n_trees)
       @inbounds for j in range(1, stop = n_intervals)
           if typeof(random_state) == StableRNGs.LehmerRNG
               Intervals[i,j,1] = rand(random_state, 1:(series_length - min_interval))
               len = rand(random_state, 1:(series_length - Intervals[i,j,1]))
           else
               Intervals[i,j,1] = rand(1:(series_length - min_interval))
               len = rand(1:(series_length - Intervals[i,j,1]))
           end
           if len < min_interval
               len = min_interval
           end
           Intervals[i,j,2] = Intervals[i,j,1] + len
           Y = @views Xmatrix[:, Intervals[i,j,1]:Intervals[i,j,2]]
           x = Array(1:size(Y)[2])
           Interval_features[:, 3*j-2, i] =  mean(Y, dims=2)                              #means
           Interval_features[:, 3*j-1, i] =  std(Y, dims=2)                               #stds
           Interval_features[:, 3*j, i]   =  (mean(transpose(x).*Y, dims=2) -             #slope
                                         mean(x)*mean(Y, dims=2)) / (mean(x.*x) - mean(x)^2)
       end
    end
    return Interval_features, Intervals
end

function InvFeatureGen(Xmatrix::Array, Intervals::Array, n_trees::Int)
    n_instance, series_length = size(Xmatrix)
    n_intervals = floor(Int64, sqrt(series_length))
    Interval_features = Array{Float64,3}(undef, n_instance, 3*n_intervals, n_trees)
    @inbounds for i in range(1, stop = n_trees)
       @inbounds for j in range(1, stop = n_intervals)
           Y = @views Xmatrix[:, Intervals[i,j,1]:Intervals[i,j,2]]
           x = Array(1:size(Y)[2])
           Interval_features[:, 3*j-2, i] =  mean(Y, dims=2)                              #means
           Interval_features[:, 3*j-1, i] =  std(Y, dims=2)                               #stds
           Interval_features[:, 3*j, i]   =  (mean(transpose(x).*Y, dims=2) -             #slope
                                         mean(x)*mean(Y, dims=2)) / (mean(x.*x) - mean(x)^2)
       end
    end
    return Interval_features
end

end # TimeSeriesForestClassifier


text:
# UnitCommitment.jl: Optimization Package for Security-Constrained Unit Commitment
# Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
# Released under the modified BSD license. See COPYING.md for more details.

function _add_status_vars!(
    model::JuMP.Model,
    g::Unit,
    formulation_status_vars::Gar1962.StatusVars,
)::Nothing
    is_on = _init(model, :is_on)
    switch_on = _init(model, :switch_on)
    switch_off = _init(model, :switch_off)
    for t in 1:model[:instance].time
        if g.must_run[t]
            is_on[g.name, t] = 1.0
            switch_on[g.name, t] = (t == 1 ? 1.0 - _is_initially_on(g) : 0.0)
            switch_off[g.name, t] = 0.0
        else
            is_on[g.name, t] = @variable(model, binary = true)
            switch_on[g.name, t] = @variable(model, binary = true)
            switch_off[g.name, t] = @variable(model, binary = true)
        end
    end
    return
end

function _add_status_eqs!(
    model::JuMP.Model,
    g::Unit,
    formulation_status_vars::Gar1962.StatusVars,
)::Nothing
    eq_binary_link = _init(model, :eq_binary_link)
    eq_switch_on_off = _init(model, :eq_switch_on_off)
    is_on = model[:is_on]
    switch_off = model[:switch_off]
    switch_on = model[:switch_on]
    for t in 1:model[:instance].time
        if !g.must_run[t]
            # Link binary variables
            if t == 1
                eq_binary_link[g.name, t] = @constraint(
                    model,
                    is_on[g.name, t] - _is_initially_on(g) ==
                    switch_on[g.name, t] - switch_off[g.name, t]
                )
            else
                eq_binary_link[g.name, t] = @constraint(
                    model,
                    is_on[g.name, t] - is_on[g.name, t-1] ==
                    switch_on[g.name, t] - switch_off[g.name, t]
                )
            end
            # Cannot switch on and off at the same time
            eq_switch_on_off[g.name, t] = @constraint(
                model,
                switch_on[g.name, t] + switch_off[g.name, t] <= 1
            )
        end
    end
    return
end


text:
import .StaticArrays: MArray

import Base: Matrix, Array

Array{T, N}(::UndefInitializer, d::NTuple{N,StaticInteger}) where {T,N} =  MArray{Tuple{Int.(d)...}, T, N, prod(d)}(undef)
Array{T, 1}(::UndefInitializer, ::Tuple{StaticInteger{D}}) where {T, D} = MArray{Tuple{Int(D)}, T, 1, Int(D)}(undef)
Array{T, N}(::UndefInitializer, d::Vararg{StaticInteger,N}) where {T, N} = Array{T, N}(undef, d)
Matrix(::UndefInitializer, m::StaticInteger, n::StaticInteger) = Matrix{Any}(undef, m, n)
Vector(::UndefInitializer, n::StaticInteger) = Vector{Any}(undef, n)


text:
action_request_test(line0, char0, line1=line0, char1=char0; diags=[]) = LanguageServer.textDocument_codeAction_request(LanguageServer.CodeActionParams(LanguageServer.TextDocumentIdentifier("testdoc"), LanguageServer.Range(LanguageServer.Position(line0, char0), LanguageServer.Position(line1, char1)), LanguageServer.CodeActionContext(diags, missing)), server, server.jr_endpoint)

@testset "reexport" begin
    settestdoc("using Base.Meta\n")
    @test any(c.command == "ReexportModule" for c in action_request_test(0, 15))
    c = filter(c -> c.command == "ReexportModule", action_request_test(0, 15))[1]

    LanguageServer.workspace_executeCommand_request(LanguageServer.ExecuteCommandParams(missing, c.command, c.arguments), server, server.jr_endpoint)
end

@testset "inline expand" begin
    settestdoc("f(x) = x")
    @test any(c.command == "ExpandFunction" for c in action_request_test(0, 5))
    c = filter(c -> c.command == "ExpandFunction", action_request_test(0, 5))[1]

    LanguageServer.workspace_executeCommand_request(LanguageServer.ExecuteCommandParams(missing, c.command, c.arguments), server, server.jr_endpoint)

    settestdoc("f(x) = begin x end")
    @test any(c.command == "ExpandFunction" for c in action_request_test(0, 5))
    c = filter(c -> c.command == "ExpandFunction", action_request_test(0, 5))[1]

    LanguageServer.workspace_executeCommand_request(LanguageServer.ExecuteCommandParams(missing, c.command, c.arguments), server, server.jr_endpoint)
end

@testset "fixmissingref" begin
    doc = settestdoc("argtail\n")
    e = LanguageServer.mark_errors(doc)[1]
    @test any(c.command == "FixMissingRef" for c in action_request_test(0, 5, diags=[e]))
    c = filter(c -> c.command == "FixMissingRef", action_request_test(0, 5, diags=[e]))[1]

    LanguageServer.workspace_executeCommand_request(LanguageServer.ExecuteCommandParams(missing, c.command, c.arguments), server, server.jr_endpoint)
end

@testset "explicit import" begin
    doc = settestdoc("using Base.Meta\nMeta.quot")
    @test LanguageServer.find_using_statement(doc.cst.args[2].args[1]) !== nothing

    @test any(c.command == "ExplicitPackageVarImport" for c in action_request_test(1, 1))
    c = filter(c -> c.command == "ExplicitPackageVarImport", action_request_test(1, 1))[1]

    LanguageServer.workspace_executeCommand_request(LanguageServer.ExecuteCommandParams(missing, c.command, c.arguments), server, server.jr_endpoint)
end


text:
# 1 -> 2
#      |
#      |
#      v
#      3
station1 = EuclideanStation(1, EuclideanCoordinate(1, 4), "test", false)
station2 = EuclideanStation(2, EuclideanCoordinate(2, 4), "test", false)
station3 = EuclideanStation(4, EuclideanCoordinate(2, 2), "test", false)
line1 = Line("S1")
edge = ProcessedEdge(station1, station2, line1, 0, 1, false)
edge2 = ProcessedEdge(station2, station3, line1, 6, 1, false)
transit_map = InputGraph([station1, station2, station3],
    [edge, edge2],
    [line1])

result = optimize(GLPKSolverMIP(msg_lev = 3), transit_map, false)

@test result.stations[1].coordinate.y > result.stations[2].coordinate.y
@test result.stations[1].coordinate.x < result.stations[2].coordinate.x
@test result.stations[2].coordinate.x < result.stations[3].coordinate.x
@test result.stations[2].coordinate.y > result.stations[3].coordinate.y

max_stars_repo_name: brainfo/MetroMap

max_issues_repo_name: brainfo/MetroMap

max_forks_repo_name: brainfo/MetroMap


text:
module BarotropicQGQL

export
  Problem,
  set_zeta!,
  updatevars!,

  energy,
  enstrophy,
  dissipation,
  work,
  drag

using
  FFTW,
  Reexport

@reexport using FourierFlows

using LinearAlgebra: mul!, ldiv!
using FourierFlows: parsevalsum, parsevalsum2
import FFTW: rfft

abstract type BarotropicQGQLVars <: AbstractVars end
abstract type BarotropicQGQLForcedVars <: BarotropicQGQLVars end

const physicalvars = [:zeta, :psi, :u, :v, :uzeta, :vzeta, :U, :Zeta, :Psi]
const transformvars = [ Symbol(var, :h) for var in physicalvars ]
const forcedvars = [:Fh]
const stochforcedvars = [:prevsol]

nothingfunction(args...) = nothing

"""
    Problem(; parameters...)

Construct a BarotropicQGQL turbulence problem.
"""
function Problem(;
    # Numerical parameters
            nx = 256,
            Lx = 2π,
            ny = nx,
            Ly = Lx,
            dt = 0.01,
    # Physical parameters
            f0 = 1.0,
          beta = 0.0,
           eta = nothing,
    # Drag and/or hyper-/hypo-viscosity
            nu = 0.0,
           nnu = 1,
            mu = 0.0,
   # Timestepper and eqn options
       stepper = "RK4",
         calcF = nothingfunction,
    stochastic = false,
             T = Float64)

  # the grid
  gr  = TwoDGrid(nx, Lx, ny, Ly; T=T)
  x, y = gridpoints(gr)

  # topographic PV
  if eta==nothing
    eta = 0*x
    etah = rfft(eta)
  end

  if typeof(eta)!=Array{T, 2} #this is true if eta was passes in Problem as a function
    pr = Params(gr, f0, beta, eta, mu, nu, nnu, calcF)
  else
    pr = Params{T}(f0, beta, eta, rfft(eta), mu, nu, nnu, calcF)
  end
  vs = calcF == nothingfunction ? Vars(gr) : (stochastic ? StochasticForcedVars(gr) : ForcedVars(gr))
  eq = BarotropicQGQL.Equation(pr, gr)
  FourierFlows.Problem(eq, stepper, dt, gr, vs, pr)
end

InitialValueProblem(; kwargs...) = Problem(; kwargs...)
ForcedProblem(; kwargs...) = Problem(; kwargs...)

# ----------
# Parameters
# ----------

"""
    Params(g::TwoDGrid, f0, beta, FU, eta, mu, nu, nnu)

Returns the params for an unforced two-dimensional barotropic QG problem.
"""
struct Params{T} <: AbstractParams
  f0::T                      # Constant planetary vorticity
  beta::T                    # Planetary vorticity y-gradient
  eta::Array{T, 2}           # Topographic PV
  etah::Array{Complex{T}, 2} # FFT of Topographic PV
  mu::T                      # Linear drag
  nu::T                      # Viscosity coefficient
  nnu::Int                   # Hyperviscous order (nnu=1 is plain old viscosity)
  calcF!::Function           # Function that calculates the forcing on QGPV q
end

"""
    Params(g::TwoDGrid, f0, beta, eta::Function, mu, nu, nnu)

Constructor for Params that accepts a generating function for the topographic PV.
"""
function Params(g::AbstractGrid{T}, f0, beta, eta::Function, mu, nu, nnu, calcF) where T
  x, y = gridpoints(g)
  etagrid = eta(x, y)
  etah = rfft(etagrid)
  Params{T}(f0, beta, etagrid, etah, mu, nu, nnu, calcF)
end

# ---------
# Equations
# ---------

"""
    Equation(p, g)

Returns the equation for two-dimensional barotropic QG problem with params p and grid g.
"""
function Equation(p::Params, g::AbstractGrid{T}) where T
  L = @. -p.mu - p.nu*g.Krsq^p.nnu + im*p.beta*g.kr*g.invKrsq
  L[1, 1] = 0
  FourierFlows.Equation(L, calcN!, g)
end

# ----
# Vars
# ----

"""
    Vars(g)

Returns the vars for unforced two-dimensional barotropic QG problem with grid g.
"""
mutable struct Vars{T} <: BarotropicQGQLVars
  u::Array{T,2}
  v::Array{T,2}
  U::Array{T,2}
  uzeta::Array{T,2}
  vzeta::Array{T,2}
  zeta::Array{T,2}
  Zeta::Array{T,2}
  psi::Array{T,2}
  Psi::Array{T,2}
  Nz::Array{Complex{T},2}
  NZ::Array{Complex{T},2}
  uh::Array{Complex{T},2}
  vh::Array{Complex{T},2}
  Uh::Array{Complex{T},2}
  zetah::Array{Complex{T},2}
  Zetah::Array{Complex{T},2}
  psih::Array{Complex{T},2}
  Psih::Array{Complex{T},2}
end

function Vars(g::AbstractGrid{T}) where T
  @createarrays T (g.nx, g.ny) u v U uzeta vzeta zeta Zeta psi Psi
  @createarrays Complex{T} (g.nkr, g.nl) N NZ uh vh Uh zetah Zetah psih Psih
  Vars(u, v, U, uzeta, vzeta, zeta, Zeta, psi, Psi, N, NZ, uh, vh, Uh, zetah, Zetah, psih, Psih)
end

"""
    ForcedVars(g)

Returns the vars for forced two-dimensional barotropic QG problem with grid g.
"""
mutable struct ForcedVars{T} <: BarotropicQGQLForcedVars
  u::Array{T,2}
  v::Array{T,2}
  U::Array{T,2}
  uzeta::Array{T,2}
  vzeta::Array{T,2}
  zeta::Array{T,2}
  Zeta::Array{T,2}
  psi::Array{T,2}
  Psi::Array{T,2}
  Nz::Array{Complex{T},2}
  NZ::Array{Complex{T},2}
  uh::Array{Complex{T},2}
  vh::Array{Complex{T},2}
  Uh::Array{Complex{T},2}
  zetah::Array{Complex{T},2}
  Zetah::Array{Complex{T},2}
  psih::Array{Complex{T},2}
  Psih::Array{Complex{T},2}
  Fh::Array{Complex{T},2}
end

function ForcedVars(g::AbstractGrid{T}) where T
  v = Vars(g)
  Fh = zeros(Complex{T}, (g.nkr, g.nl))
  ForcedVars(getfield.(Ref(v), fieldnames(typeof(v)))..., Fh)
end

mutable struct StochasticForcedVars{T} <: BarotropicQGQLForcedVars
  u::Array{T,2}
  v::Array{T,2}
  U::Array{T,2}
  uzeta::Array{T,2}
  vzeta::Array{T,2}
  zeta::Array{T,2}
  Zeta::Array{T,2}
  psi::Array{T,2}
  Psi::Array{T,2}
  Nz::Array{Complex{T},2}
  NZ::Array{Complex{T},2}
  uh::Array{Complex{T},2}
  vh::Array{Complex{T},2}
  Uh::Array{Complex{T},2}
  zetah::Array{Complex{T},2}
  Zetah::Array{Complex{T},2}
  psih::Array{Complex{T},2}
  Psih::Array{Complex{T},2}
  Fh::Array{Complex{T},2}
  prevsol::Array{Complex{T},2}
end

function StochasticForcedVars(g::AbstractGrid{T}) where T
  v = ForcedVars(g)
  prevsol = zeros(Complex{T}, (g.nkr, g.nl))
  StochasticForcedVars(getfield.(Ref(v), fieldnames(typeof(v)))..., prevsol)
end

# -------
# Solvers
# -------

function calcN_advection!(N, sol, t, cl, v, p, g)
  # Note that U = sol[1, 1]. For all other elements ζ = sol
  Kr = [ g.kr[i] for i=1:g.nkr, j=1:g.nl]
  @. v.zetah = sol
  @. v.zetah[Kr .== 0] = 0
  @. v.Zetah = sol
  @. v.Zetah[abs.(Kr) .> 0] = 0

  @. v.uh =  im * g.l  * g.invKrsq * v.zetah
  @. v.vh = -im * g.kr * g.invKrsq * v.zetah
  @. v.Uh =  im * g.l  * g.invKrsq * v.Zetah

  ldiv!(v.zeta, g.rfftplan, v.zetah)
  ldiv!(v.u, g.rfftplan, v.uh)
  ldiv!(v.v, g.rfftplan, v.vh)

  ldiv!(v.Zeta, g.rfftplan, v.Zetah)
  ldiv!(v.U, g.rfftplan, v.Uh)

  @. v.uzeta = v.u*v.zeta # u*ζ
  @. v.vzeta = v.v*v.zeta # v*ζ

  mul!(v.uh, g.rfftplan, v.uzeta) # \hat{u*q}
  @. v.NZ = -im*g.kr*v.uh # -∂[u*q]/∂x
  mul!(v.vh, g.rfftplan, v.vzeta) # \hat{v*q}
  @. v.NZ += - im*g.l*v.vh # -∂[v*q]/∂y
  @. v.NZ[abs.(Kr) .> 0] = 0

  @. v.U = v.U*v.zeta # U*ζ
  @. v.u = v.u*v.Zeta # u*Ζ
  @. v.v = v.v*v.Zeta # v*Ζ

  mul!(v.uh, g.rfftplan, v.U + v.u) # \hat{U*ζ + u*Ζ}
  @. v.Nz = -im*g.kr*v.uh # -∂[U*ζ + u*Ζ]/∂x
  mul!(v.vh, g.rfftplan, v.v) # \hat{v*Z}
  @. v.Nz += - im*g.l*v.vh # -∂[v*Z]/∂y
  @. v.Nz[abs.(Kr) .== 0] = 0

  @. N = v.NZ + v.Nz
end

function calcN!(N, sol, t, cl, v, p, g)
  calcN_advection!(N, sol, t, cl, v, p, g)
  addforcing!(N, sol, t, cl, v, p, g)
  nothing
end

addforcing!(N, sol, t, cl, v::Vars, p, g) = nothing

function addforcing!(N, sol, t, cl, v::ForcedVars, p, g)
  p.calcF!(v.Fh, sol, t, cl, v, p, g)
  @. N += v.Fh
  nothing
end

function addforcing!(N, sol, t, cl, v::StochasticForcedVars, p, g)
  if t == cl.t # not a substep
    @. v.prevsol = sol # sol at previous time-step is needed to compute budgets for stochastic forcing
    p.calcF!(v.Fh, sol, t, cl, v, p, g)
  end
  @. N += v.Fh
  nothing
end

# ----------------
# Helper functions
# ----------------

"""
    updatevars!(v, s, g)

Update the vars in v on the grid g with the solution in sol.
"""
function updatevars!(sol, v, p, g)
  Kr = [ g.kr[i] for i=1:g.nkr, j=1:g.nl]
  sol[1, 1] = 0
  @. v.zetah = sol
  @. v.zetah[Kr .== 0] = 0
  @. v.Zetah = sol
  @. v.Zetah[abs.(Kr) .> 0] = 0

  @. v.Psih = -v.Zetah * g.invKrsq
  @. v.psih = -v.zetah * g.invKrsq
  @. v.uh = -im * g.l  * v.psih
  @. v.vh =  im * g.kr * v.psih
  @. v.Uh =  im * g.l  * v.Zetah * g.invKrsq

  ldiv!(v.zeta, g.rfftplan, deepcopy(v.zetah))
  ldiv!(v.Zeta, g.rfftplan, deepcopy(v.Zetah))
  ldiv!(v.psi, g.rfftplan, v.psih)
  ldiv!(v.Psi, g.rfftplan, v.Psih)
  ldiv!(v.u, g.rfftplan, deepcopy(v.uh))
  ldiv!(v.v, g.rfftplan, deepcopy(v.vh))
  ldiv!(v.U, g.rfftplan, deepcopy(v.Uh))

  nothing
end

updatevars!(prob) = updatevars!(prob.sol, prob.vars, prob.params, prob.grid)

"""
    set_zeta!(prob, zeta)
    set_zeta!(sol, v, g, zeta)

Set the solution sol as the transform of zeta and update variables v
on the grid g.
"""
function set_zeta!(sol, v, p, g, zeta)
  mul!(v.zetah, g.rfftplan, zeta)
  v.zetah[1, 1] = 0.0
  @. sol = v.zetah

  updatevars!(sol, v, p, g)
  nothing
end

set_zeta!(prob, zeta) = set_zeta!(prob.sol, prob.vars, prob.params, prob.grid, zeta)

"""
Calculate the domain-averaged kinetic energy.
"""
function energy(prob)
  sol, g = prob.sol, prob.grid
  0.5*(parsevalsum2(g.kr.*g.invKrsq.*sol, g)
        + parsevalsum2(g.l.*g.invKrsq.*sol, g))/(g.Lx*g.Ly)
end

"""
Returns the domain-averaged enstrophy.
"""
function enstrophy(prob)
  sol, v, g = prob.sol, prob.vars, prob.grid
  @. v.uh = sol
  v.uh[1, 1] = 0
  0.5*parsevalsum2(v.uh, g)/(g.Lx*g.Ly)
end

"""
    dissipation(prob)
    dissipation(sol, v, p, g)

Returns the domain-averaged dissipation rate. nnu must be >= 1.
"""
@inline function dissipation(sol, v, p, g)
  @. v.uh = g.Krsq^(p.nnu-1) * abs2(sol)
  v.uh[1, 1] = 0
  p.nu/(g.Lx*g.Ly)*parsevalsum(v.uh, g)
end

@inline dissipation(prob) = dissipation(prob.sol, prob.vars, prob.params, prob.grid)

"""
    work(prob)
    work(sol, v, p, g)

Returns the domain-averaged rate of work of energy by the forcing Fh.
"""
@inline function work(sol, v::ForcedVars, g)
  @. v.uh = g.invKrsq * sol * conj(v.Fh)
  1/(g.Lx*g.Ly)*parsevalsum(v.uh, g)
end

@inline function work(sol, v::StochasticForcedVars, g)
  @. v.uh = g.invKrsq * (v.prevsol + sol)/2.0 * conj(v.Fh) # Stratonovich
  # @. v.uh = g.invKrsq * v.prevsol * conj(v.Fh)             # Ito
  1/(g.Lx*g.Ly)*parsevalsum(v.uh, g)
end

@inline work(prob) = work(prob.sol, prob.vars, prob.grid)

"""
    drag(prob)
    drag(sol, v, p, g)

Returns the extraction of domain-averaged energy by drag mu.
"""
@inline function drag(sol, v, p, g)
  @. v.uh = g.Krsq^(-1) * abs2(sol)
  v.uh[1, 1] = 0
  p.mu/(g.Lx*g.Ly)*parsevalsum(v.uh, g)
end

@inline drag(prob) = drag(prob.sol, prob.vars, prob.params, prob.grid)

end # module


text:
function call(::Types.Help, args::Vector, storage, conn::IO)
    write(conn, Types.help_msg)
end


text:
#
# Hyperparameter search
#

function loss_reduction(calibration, validation, initial_parameters, kwargs)
    params, losses, mean_vars = eki(calibration, initial_parameters; kwargs...)
    valid_loss_start = validation(initial_parameters)
    valid_loss_final = validation(params)

    @info "parameters: $params"
    @info "final loss / initial loss: $(losses[end] / losses[1])"

    train_loss_reduction = losses[end] / losses[1]
    valid_loss_reduction = valid_loss_final / valid_loss_start
    return train_loss_reduction, valid_loss_reduction
end

## Stds within bounds

function plot_stds_within_bounds(calibration, validation, initial_parameters, directory; xrange=-3:0.25:5)
    loss_reductions = Dict()
    val_loss_reductions = Dict()
    for stds_within_bounds in xrange
        train_loss_reduction, val_loss_reduction = loss_reduction(calibration, validation, initial_parameters, (stds_within_bounds = stds_within_bounds,))
        loss_reductions[stds_within_bounds] = train_loss_reduction
        val_loss_reductions[stds_within_bounds] = val_loss_reduction
    end
    p = Plots.plot(title="Loss reduction versus prior std's within bounds (n_σ)", ylabel="Loss reduction (Final / Initial)", xlabel="prior std's spanned by bound width (n_σ)", legend=false, lw=3)
    plot!(loss_reductions, label="training", color=:purple, lw=4)
    plot!(val_loss_reductions, label="validation", color=:blue, lw=4)
    Plots.savefig(p, joinpath(directory, "stds_within_bounds.pdf"))
    println("loss-minimizing stds within bounds: $(argmin(val_loss_reductions))")
end

## Prior Variance
function plot_prior_variance(calibration, validation, initial_parameters, directory; xrange=0.1:0.1:1.0)
    loss_reductions = Dict()
    val_loss_reductions = Dict()
    for variance = xrange
        train_loss_reduction, val_loss_reduction = loss_reduction(calibration, validation, initial_parameters, (stds_within_bounds = variance,))
        loss_reductions[variance] = train_loss_reduction
        val_loss_reductions[variance] = val_loss_reduction
    end
    p = Plots.plot(title="Loss reduction versus prior variance", ylabel="Loss reduction (Final / Initial)", xlabel="Prior variance")
    plot!(loss_reductions, label="training", lw=4, color=:purple)
    plot!(val_loss_reductions, label="validation", lw=4, color=:blue)
    Plots.savefig(p, joinpath(directory, "variance.pdf"))
    v = argmin(var_val_loss_reductions)
    println("loss-minimizing variance: $(v)")
    return v
end

## Number of ensemble members

function plot_num_ensemble_members(calibration, validation, initial_parameters, directory; xrange=1:5:30)
    loss_reductions = Dict()
    val_loss_reductions = Dict()
    for N_ens = xrange
        train_loss_reduction, val_loss_reduction = loss_reduction(calibration, validation, initial_parameters, (N_ens = N_ens,))
        loss_reductions[N_ens] = train_loss_reduction
        val_loss_reductions[N_ens] = val_loss_reduction
    end
    p = Plots.plot(title="Loss reduction versus N_ens", xlabel="N_ens", ylabel="Loss reduction (Final / Initial)", legend=false, lw=3)
    plot!(loss_reductions, label="training")
    plot!(val_loss_reductions, label="validation")
    Plots.savefig(p, joinpath(directory, "N_ens.pdf"))
end

## Observation noise level

function plot_observation_noise_level(calibration, validation, initial_parameters, directory; xrange=-2.0:0.2:3.0)
    loss_reductions = Dict()
    val_loss_reductions = Dict()
    for log_noise_level = xrange
        train_loss_reduction, val_loss_reduction = loss_reduction(calibration, validation, initial_parameters, (noise_level = 10.0^log_noise_level,))
        loss_reductions[log_noise_level] = train_loss_reduction
        val_loss_reductions[log_noise_level] = val_loss_reduction
    end
    p = Plots.plot(title="Loss Reduction versus Observation Noise Level", xlabel="log₁₀(Observation noise level)", ylabel="Loss reduction (Final / Initial)", legend=:topleft)
    plot!(loss_reductions, label="training", lw=4, color=:purple)
    plot!(val_loss_reductions, label="validation", lw=4, color=:blue)
    Plots.savefig(p, joinpath(directory, "obs_noise_level.pdf"))
    nl = argmin(val_loss_reductions)
    println("loss-minimizing obs noise level: $(nl)")
    return 10^nl
end

function plot_prior_variance_and_obs_noise_level(calibration, validation, initial_parameters, directory; vrange=0.40:0.025:0.90, nlrange=-2.5:0.1:0.5)
    Γθs = collect(vrange)
    Γys = 10 .^ collect(nlrange)
    losses = zeros((length(Γθs), length(Γys)))
    counter = 1
    countermax = length(Γθs)*length(Γys)
    for i in 1:length(Γθs)
        for j in 1:length(Γys)
            println("progress $(counter)/$(countermax)")
            Γθ = Γθs[i]
            Γy = Γys[j]
            train_loss_reduction, val_loss_reduction = loss_reduction(calibration, validation, initial_parameters, (stds_within_bounds=Γθ, noise_level=Γy, N_iter=10, N_ens=10, informed_priors=false))
            losses[i, j] = val_loss_reduction
            counter += 1
        end
    end
    p = Plots.heatmap(Γys, Γθs, losses, xlabel=L"\Gamma_y", ylabel=L"\Gamma_\theta", size=(250,250), yscale=:log10)
    Plots.savefig(p, joinpath(directory, "GammaHeatmap.pdf"))
    v = Γθs[argmin(losses)[1]]
    nl = Γys[argmin(losses)[2]]
    println("loss-minimizing Γθ: $(v)")
    println("loss-minimizing log10(Γy): $(log10(nl))")
    return v, nl
end


text:
module HaploSoup

export build_prefix_and_divergence_arrays
export build_encoding
encode_haplotype_matrix

include("encoding.jl")

end


text:
module Fortran

# const Label = String
#
# const TypeSpec = Union{}
#
# const Body = Union{}
#
# const MainProgram = Union{}
# const SubroutineSubprogram = Union{}
# const BlockDataSubprogram = Union{}
#
# struct Ident
#     ident::String
# end
#
# struct FunctionRange
#     functionParList::Vector{Ident}
#     body::Body
# end
#
# struct FunctionSubprogram
#     label::Union{Nothing, Label}
#     typeSpec::Union{Nothing, TypeSpec}
#     functionName::Ident
#     functionRange::FunctionRange
# end
#
# struct ProgramUnit
#     data::Union{MainProgram,FunctionSubprogram,SubroutineSubprogram,BlockDataSubprogram}
# end
#
# struct ExecutableProgram
#     data::Vector{ProgramUnit}
# end

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

clean_code(expr) = expr
function clean_code(expr::Expr)
    expr = Expr(expr.head, map(clean_code, filter(arg -> !(arg isa LineNumberNode), expr.args))...)
    # Remove line numbers.
    # Line numbers are usually wrong because they point to this file,
    # instead of the file where the code originates.
    if expr.head ≡ :block && length(expr.args) == 1
        expr = expr.args[1]
    end
    return expr
end
export clean_code

function feval end
export feval

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

struct FLogical
    value::Int32
    FLogical(b::Bool) = new(b)
end
export FLogical
Base.Bool(l::FLogical) = Bool(l.value)
Base.:%(l::FLogical, ::Type{Bool}) = l.value % Bool

const FInteger = Int32
export FInteger
const FReal = Float32
export FReal
const FDoublePrecision = Float64
export FDoublePrecision
const FComplex = Complex{FReal}
export FComplex
const FDoubleComplex = Complex{FDoublePrecision}
export FDoubleComplex

struct Fun{Domain<:Tuple,Codomain} <: Function
    fun::Function
end
export Fun
@inline (f::Fun{Domain,Codomain})(xs...) where {Domain,Codomain} = f.fun((xs::Domain)...)::Codomain

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

abstract type FExpr end
export FExpr

struct Const <: FExpr
    value::Any
    type::Type
    Const(value::T, ::Type{T}) where {T} = new(value, T)
end
export Const
Const(value) = Const(value, typeof(value))
feval(c::Const) = c.value

struct Var <: FExpr
    name::Symbol
    type::Type
end
export Var
feval(v::Var) = :($(v.name)::$(v.type))

struct Call <: FExpr
    fun::FExpr
    args::Vector{FExpr}
end
export Call
feval(c::Call) = :($(feval(c.fun))($((feval(arg) for arg in c.args)...)))

abstract type Stmt end
export Stmt

struct Assign <: Stmt
    lhs::Var
    rhs::FExpr
end
export Assign
feval(a::Assign) = quote
    $(a.lhs.name) = $(a.lhs.type)($(feval(a.rhs)))
end

struct Block <: Stmt
    stmts::Vector{Stmt}
end
export Block
feval(b::Block) = quote
    $((feval(stmt) for stmt in b.stmts)...)
end

struct Do <: Stmt
    var::Var
    lbnd::FExpr
    ubnd::FExpr
    step::FExpr
    body::Stmt
end
export Do
function feval(d::Do)
    name = d.var.name
    type = d.var.type
    lbnd = Symbol(name, ".lbnd")
    ubnd = Symbol(name, ".ubnd")
    step = Symbol(name, ".step")
    count = Symbol(name, ".count")
    idx = Symbol(name, ".idx")
    quote
        $lbnd = $(feval(d.lbnd))
        $ubnd = $(feval(d.ubnd))
        $step = $(feval(d.step))
        $step == 0 && error("Do loop step is zero")
        $count = FInteger(($ubnd - $lbnd) ÷ $step)
        $name = $type($lbnd)
        for $idx in FInteger(0):($count)
            $(feval(d.body))
            $name += $type($step)
        end
    end
end

struct DoWhile <: Stmt
    cond::FExpr
    body::Stmt
end
export DoWhile
function feval(d::DoWhile)
    quote
        while $(feval(d.cond))::FLogical % Bool
            $(feval(d.body))
        end
    end
end

struct If <: Stmt
    branches::Vector{Pair{FExpr,Stmt}}
    default::Stmt
end
export If
function feval(i::If)
    stmt = feval(i.default)
    for branch in Iterators.reverse(i.branches)
        stmt = quote
            if $(feval(branch[1]))::FLogical % Bool
                $(feval(branch[2]))
            else
                $stmt
            end
        end
    end
    return stmt
end

struct Print <: Stmt
    format::Any                 # TODO
    exprs::Vector{FExpr}
end
export Print
feval(p::Print) = quote
    println(" ", $((feval(expr) for expr in p.exprs)...))
end

struct FFunction
    name::Symbol
    type::Type
    args::Vector{Var}
    vars::Vector{Var}
    body::Stmt
end
export FFunction
feval(f::FFunction) = quote
    function $(f.name)($((:($(arg.name)::$(arg.type)) for arg in f.args)...))
        $(feval(f.body))
        return $(Symbol(f.name, ".result"))::$(f.type)
    end
end

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

feq(x::T, y::T) where {T} = FLogical(x == y)
fne(x::T, y::T) where {T} = FLogical(x ≠ y)
fle(x::T, y::T) where {T} = FLogical(x ≤ y)
flt(x::T, y::T) where {T} = FLogical(x < y)
fge(x::T, y::T) where {T} = FLogical(x ≥ y)
fgt(x::T, y::T) where {T} = FLogical(x > y)

const pos_integer = Const(Fun{Tuple{FInteger},FInteger}(+))
const neg_integer = Const(Fun{Tuple{FInteger},FInteger}(-))

const add_integer = Const(Fun{Tuple{FInteger,FInteger},FInteger}(+))
const sub_integer = Const(Fun{Tuple{FInteger,FInteger},FInteger}(-))
const mul_integer = Const(Fun{Tuple{FInteger,FInteger},FInteger}(*))
const div_integer = Const(Fun{Tuple{FInteger,FInteger},FInteger}(÷))
const mod_integer = Const(Fun{Tuple{FInteger,FInteger},FInteger}(%))

const eq_integer = Const(Fun{Tuple{FInteger,FInteger},FLogical}(feq))
const ne_integer = Const(Fun{Tuple{FInteger,FInteger},FLogical}(fne))
const le_integer = Const(Fun{Tuple{FInteger,FInteger},FLogical}(fle))
const lt_integer = Const(Fun{Tuple{FInteger,FInteger},FLogical}(flt))
const ge_integer = Const(Fun{Tuple{FInteger,FInteger},FLogical}(fge))
const gt_integer = Const(Fun{Tuple{FInteger,FInteger},FLogical}(fgt))

const pos_real = Const(Fun{Tuple{FReal},FReal}(+))
const neg_real = Const(Fun{Tuple{FReal},FReal}(-))

const add_real = Const(Fun{Tuple{FReal,FReal},FReal}(+))
const sub_real = Const(Fun{Tuple{FReal,FReal},FReal}(-))
const mul_real = Const(Fun{Tuple{FReal,FReal},FReal}(*))
const div_real = Const(Fun{Tuple{FReal,FReal},FReal}(/))

const eq_real = Const(Fun{Tuple{FReal,FReal},FLogical}(feq))
const ne_real = Const(Fun{Tuple{FReal,FReal},FLogical}(fne))
const le_real = Const(Fun{Tuple{FReal,FReal},FLogical}(fle))
const lt_real = Const(Fun{Tuple{FReal,FReal},FLogical}(flt))
const ge_real = Const(Fun{Tuple{FReal,FReal},FLogical}(fge))
const gt_real = Const(Fun{Tuple{FReal,FReal},FLogical}(fgt))

end


text:
if abspath(PROGRAM_FILE) == @__FILE__
  if length(ARGS) < 1
    return
  end
  open(Base.ARGS[1]) do f
    lastnums = [0,0,0]
    linecount = 0
    result = 0
    for line in eachline(f)
      linecount += 1
      n = parse(Int, line)
      if (sum(lastnums[2:end]) + n > sum(lastnums)) &
         (linecount > 3)
        result += 1
      end
      lastnums[1] = lastnums[2]
      lastnums[2] = lastnums[3]
      lastnums[3] = n
    end
    println(result)
  end
end


text:
using StaticArrays
using AbstractTrees
using DataStructures

function small_example()
    R, (x, y) = Singular.PolynomialRing(Singular.Fp(101), ["x", "y"])
    I = [x^2, x*y + y^2]
    order = SG.Grevlex(2)
    dat = SG.f5data(I, order=order)
    ctx = dat.ctx
    basis = SG.new_basis(ctx, 2)
    syz = SG.new_syz(ctx, 2)
    for i in 1:2
        SG.new_basis_elem!(ctx, basis, (SG.pos_type(ctx)(i), ctx.po.mo(R(1))))
    end
    R, (x, y), ctx, basis, syz
end

function contained_in(I1::sideal{MP},
                      I2::sideal{MP}) where {MP <: Singular.MPolyElem}

    gb = std(I2)
    all(p -> iszero(reduce(p, gb)), gens(I1))
end

function is_eq(I1::sideal{MP},
               I2::sideal{MP}) where {MP <: Singular.MPolyElem}

    contained_in(I1, I2) && contained_in(I2, I1)
end

function fancy_loop(I::Vector{MP};
                    interreduce = false) where {MP <: Singular.MPolyElem}

    R = parent(first(I))
    I_prime = Ideal(R, first(I))
    Js = typeof(I_prime)[]
    for f in I[2:end]
        f_id = Ideal(R, f)
        I_sat = saturation(I_prime, f_id)[1]
        if !(is_eq(I_prime, I_sat))
            push!(Js, saturation(I_prime, I_sat)[1])
        end
        I_prime = I_sat + f_id
    end

    for J in Js
        I_prime = saturation(I_prime, J)[1]
    end
    interreduce ? std(I_prime, complete_reduction = true) : I_prime
end

@testset "termorder" begin
    order = SG.Grevlex(5)
    v = @SVector [2,2,3,4,5]
    w = @SVector [3,1,2,5,5]

    @test SG.lt(order, w, v)
    @test SG.iscompatible(order, v, w)
    @test !(SG.divides(order, v, w))
end

# stolen from pierre
@testset "polynomials" begin
    char = 65521
    ctx = SG.polynomialctx(SG.Nmod32Γ(char), monomials = SG.monomialctx(order=SG.Grevlex(5)))

    # Conversion from and to AA
    R, x = SG.abstractalgebra(ctx)

    p1 = (1+x[1]+x[2]+x[3])^10
    p0 = ctx(p1)

    @test R(ctx, p0) == p1
    @test SG.monomial(p0, 1) == ctx.mo(x[1]^10)
    @test length(p0) == binomial(10+3,3)

    # LCM
    @test SG.lcm(ctx.mo, ctx.mo([1,0,0,2,0]), ctx.mo([0,1,0,0,0])) == ctx.mo([1,1,0,2,0])
end

@testset "sig polynomials" begin
    R, (x, y) = Singular.PolynomialRing(Singular.Fp(101), ["x", "y"])
    f, g = x + y, x^2 + x*y + y^2

    order = SG.Grevlex(2)
    char = 101
    ctx = SG.idxsigpolynomialctx(SG.Nmod32Γ(char), 2, order=order)

    sig1, sig2 = ctx(1, R(1)), ctx(2, R(1))
    ctx(sig1, f), ctx(sig2, g)
    m1 = ctx.po.mo(x)
    @test collect(keys(ctx.tbl)) == [sig1, sig2]
    @test R(ctx.po, ctx(sig1)[:poly]) == f
    @test R(ctx.po, ctx(m1, sig1)[:poly]) == x*f
end

@testset "f5 data" begin
    R, (x, y) = Singular.PolynomialRing(Singular.Fp(101), ["x", "y"])
    I = [x^2, y^2 + x*y]
    order = SG.Grevlex(2)
    dat = SG.f5data(I, order=order)
    sig1, sig2 = dat.ctx(1, R(1)), dat.ctx(2, R(1))
    @test SG.lt(dat.ctx, sig1, sig2)
end

@testset "pairs" begin
    R, (x, y), ctx, basis, syz = small_example()
    new_sig = ctx(2, x)
    ctx(new_sig, y^3)
    pairset = SG.pairset(ctx)
    SG.pairs!(ctx, pairset, new_sig, ctx.po.mo(y^3), basis, syz)
    @test isempty(pairset)
end

@testset "symbolic-pp" begin
    R, (x, y), ctx, basis, syz = small_example()
    pair_sig = (ctx.po.mo(x), ctx(2, R(1)))
    pair_sig_2 = (ctx.po.mo(y), ctx(1, R(1)))
    pairset = SG.mpairset(ctx, [pair_sig, pair_sig_2])
    SG.symbolic_pp!(ctx, pairset, basis, syz)
    test_sig_1 = (ctx.po.mo(y), ctx(1, R(1)))
    test_sig_2 = (ctx.po.mo(y), ctx(2, R(1)))
    test_pairset = SG.mpairset(ctx, [pair_sig, test_sig_1, test_sig_2])
    @test pairset == test_pairset
end

@testset "small-reduction" begin
    R, (x, y), ctx, basis, syz = small_example()
    pair_sig = (ctx.po.mo(x), ctx(2, R(1)))
    pairset = SG.mpairset(ctx)
    push!(pairset, pair_sig)
    two = SG.pos_type(ctx)(2)
    mons = SG.symbolic_pp!(ctx, pairset, basis, syz, are_pairs = false)
    rows = sort(collect(pairset), lt = (a, b) -> Base.Order.lt(SG.mpairordering(ctx), a, b))
    mat = SG.f5matrix(ctx, mons, rows, two, two, :f, trace_sig_tail_tags = [:f, :g])
    @test SG.mat_show(mat) == [1 0 0; 0 1 1; 1 1 0]
    SG.reduction!(mat)
    @test SG.mat_show(mat) == [1 0 0; 0 1 1; 0 0 100]
    p = SG.unindexpolynomial(mat.sigtail_mat.tbl,
                             mat.sigtail_mat.rows[pair_sig])
    @test R(ctx.po, p) == x - y
end

@testset "monomial hashing" begin
    order = SG.Grevlex(5)
    ctx = SG.monomialctx(exponents = Int64, order=SG.Grevlex(5))
    idx = SG.ixmonomialctx(ctx)

    v = @SVector [2,2,3,4,5]
    w = @SVector [3,1,2,5,5]

    m = SG.Monomial(v)
    n = SG.Monomial(w)
    i, j  = idx(m), idx(n)
    k = SG.mul(idx, i, j)
    @test SG.divides(idx, i, k)
    @test SG.divides(idx, j, k)
end

@testset "small groebner" begin
    R, (x, y) = Singular.PolynomialRing(Singular.Fp(101), ["x", "y"])
    I = [x^2, x*y + y^2]
    gb_2 = SG.f5(I, interreduction = false)
    gb = vcat(I, [-y^3])
    @test all(p -> p in gb, gb_2)
end

@testset "small groebner 2" begin
    R, (x, y, z, t) = Singular.PolynomialRing(Singular.Fp(7), ["x", "y", "z", "t"])
    I = [y*z - 2*t^2, x*y + t^2, x^2*z + 3*x*t^2 - 2*y*t^2]
    gb = SG.f5(I, interreduction = false)
    @test length(gb) == 7
end

@testset "cyclic 4" begin
    R, (x, y, z, w) = Singular.PolynomialRing(Singular.Fp(101), ["x", "y", "z", "w"])
    I = SG.cyclic([x,y,z,w])
    gb = SG.f5(I, interreduction = false, verbose = 1)
    @test SG.is_gb(gb)
end

@testset "cyclic 4 sigtails" begin
    R, (x, y, z, w) = Singular.PolynomialRing(Singular.Fp(101), ["x", "y", "z", "w"])
    I = SG.cyclic([x,y,z,w])
    dat = SG.f5setup(I, trace_sig_tail_tags = [:f, :g])
    G, H, pairs = SG.pairs_and_basis(dat, length(I))
    SG.f5core!(dat, G, H, pairs, interreduction = false)
    basis_sigs = [(i, g[1]) for i in keys(G) for g in G[i]]
    gb = [R(dat.ctx, s) for s in basis_sigs]
    projs = [R(dat.ctx.po, SG.project(dat.ctx, s)) for s in basis_sigs]
    gb_s = [std(Ideal(R, I[1:k])) for k in 1:3]
    for (p, g, sig) in zip(projs, gb, basis_sigs)
        isone(sig[1]) && continue
        k = sig[1]
        @test iszero(Singular.reduce(p*I[k] - g, gb_s[k - 1]))
    end
end

@testset "cyclic 6 sigtails" begin
    R, (x1, x2, x3, x4, x5, x6) = Singular.PolynomialRing(Singular.Fp(101), ["x$(i)" for i in 1:6])
    I = SG.cyclic([x1,x2,x3,x4,x5,x6])
    dat = SG.f5setup(I, trace_sig_tail_tags = [:f, :g])
    G, H, pairs = SG.pairs_and_basis(dat, length(I))
    SG.f5core!(dat, G, H, pairs, interreduction = false)
    basis_sigs = [(i, g[1]) for i in keys(G) for g in G[i]]
    gb = [R(dat.ctx, s) for s in basis_sigs]
    projs = [R(dat.ctx.po, SG.project(dat.ctx, s)) for s in basis_sigs]
    gb_s = [std(Ideal(R, I[1:k])) for k in 1:5]
    for (p, g, sig) in zip(projs, gb, basis_sigs)
        isone(sig[1]) && continue
        k = sig[1]
        @test iszero(Singular.reduce(p*I[k] - g, gb_s[k - 1]))
    end
end

@testset "cyclic 6" begin
    R, (x1, x2, x3, x4, x5, x6) = Singular.PolynomialRing(Singular.Fp(101), ["x$(i)" for i in 1:6])
    I = SG.cyclic([x1,x2,x3,x4,x5,x6])
    gb = SG.f5(I, interreduction = false)
    @test SG.is_gb(gb)
end

@testset "katsura 6 interreduction" begin
    R, (x1, x2, x3, x4, x5, x6) = Singular.PolynomialRing(Singular.Fp(65521), ["x$(i)" for i in 1:6])
    I = [x1+2*x2+2*x3+2*x4+2*x5+2*x6-1,
         x1^2+2*x2^2+2*x3^2+2*x4^2+2*x5^2+2*x6^2-x1,
         2*x1*x2+2*x2*x3+2*x3*x4+2*x4*x5+2*x5*x6-x2,
         x2^2+2*x1*x3+2*x2*x4+2*x3*x5+2*x4*x6-x3,
         2*x2*x3+2*x1*x4+2*x2*x5+2*x3*x6-x4,
         x3^2+2*x2*x4+2*x1*x5+2*x2*x6-x5]
    gb = SG.f5(I)
    @test length(gb) == 22
end

@testset "small-decomp" begin
    R, (x, y, z) = Singular.PolynomialRing(Singular.Fp(101), ["x", "y", "z"])
    I = [x*y, x*z]
    gb = decompose(I, interreduction = false)
    @test SG.is_gb(gb)
    @test is_eq(Ideal(R, gb), fancy_loop(I))
end

@testset "cyclic 4 decomp" begin
    R, (x, y, z, w) = Singular.PolynomialRing(Singular.Fp(101), ["x", "y", "z", "w"])
    I = SG.cyclic([x,y,z,w])
    gb = decompose(I, interreduction = false)
    @test SG.is_gb(gb)
    @test is_eq(Ideal(R, gb), fancy_loop(I))
end

@testset "cyclic 6 decomp" begin
    R, (x1, x2, x3, x4, x5, x6) = Singular.PolynomialRing(Singular.Fp(101), ["x$(i)" for i in 1:6])
    I = SG.cyclic([x1,x2,x3,x4,x5,x6])
    gb = decompose(I, interreduction = false)
    @test SG.is_gb(gb)
    @test is_eq(Ideal(R, gb), fancy_loop(I))
end

@testset "singhypcritpoints decomp" begin
    R, (x1, x2, x3, x4)  = Singular.PolynomialRing(Singular.Fp(65521), ["x1", "x2", "x3", "x4"])
    I = [2297232*x1^4+10816588*x1^3*x2-19696652*x1^3*x3+12014836*x1^3*x4+19038581*x1^2*x2^2+25140479*x1^2*x2*x3-124056256*x1^2*x2*x4-19352961*x1^2*x3^2+30086207*x1^2*x3*x4-15861565*x1^2*x4^2-64545276*x1*x2^3-3395230*x1*x2^2*x3+125153935*x1*x2^2*x4-10319588*x1*x2*x3^2+5839912*x1*x2*x3*x4+112915639*x1*x2*x4^2+12114032*x1*x3^3+32314614*x1*x3^2*x4-5646418*x1*x3*x4^2-3823868*x1*x4^3+40572828*x2^4-42051502*x2^3*x3-17361686*x2^3*x4+66618627*x2^2*x3^2-61020343*x2^2*x3*x4-66491064*x2^2*x4^2+16298428*x2*x3^3+46623254*x2*x3^2*x4-26720875*x2*x3*x4^2-16047619*x2*x4^3+15073900*x3^4+9644976*x3^3*x4-9790086*x3^2*x4^2-9372130*x3*x4^3+4661500*x4^4-39813328*x1^3-9400295*x1^2*x2-1289811*x1^2*x3+64756487*x1^2*x4-41226416*x1*x2^2+22255617*x1*x2*x3+115474646*x1*x2*x4+54583668*x1*x3^2+150465505*x1*x3*x4-19224417*x1*x4^2+52554279*x2^3-1318549*x2^2*x3-91874169*x2^2*x4+119183502*x2*x3^2-131012023*x2*x3*x4-88321124*x2*x4^2+41605616*x3^3+58052114*x3^2*x4-70772721*x3*x4^2+3123219*x4^3+99968607*x1^2-17533608*x1*x2+133696391*x1*x3-71771275*x1*x4+20484256*x2^2+10206682*x2*x3-58141098*x2*x4+72704712*x3^2-111938286*x3*x4-5946966*x4^2-46338000*x1+12230767*x2-23665371*x3+7665457*x4+9313394, 10816588*x1^3+38077162*x1^2*x2+25140479*x1^2*x3-124056256*x1^2*x4-193635828*x1*x2^2-6790460*x1*x2*x3+250307870*x1*x2*x4-10319588*x1*x3^2+5839912*x1*x3*x4+112915639*x1*x4^2+162291312*x2^3-126154506*x2^2*x3-52085058*x2^2*x4+133237254*x2*x3^2-122040686*x2*x3*x4-132982128*x2*x4^2+16298428*x3^3+46623254*x3^2*x4-26720875*x3*x4^2-16047619*x4^3-9400295*x1^2-82452832*x1*x2+22255617*x1*x3+115474646*x1*x4+157662837*x2^2-2637098*x2*x3-183748338*x2*x4+119183502*x3^2-131012023*x3*x4-88321124*x4^2-17533608*x1+40968512*x2+10206682*x3-58141098*x4+12230767, -19696652*x1^3+25140479*x1^2*x2-38705922*x1^2*x3+30086207*x1^2*x4-3395230*x1*x2^2-20639176*x1*x2*x3+5839912*x1*x2*x4+36342096*x1*x3^2+64629228*x1*x3*x4-5646418*x1*x4^2-42051502*x2^3+133237254*x2^2*x3-61020343*x2^2*x4+48895284*x2*x3^2+93246508*x2*x3*x4-26720875*x2*x4^2+60295600*x3^3+28934928*x3^2*x4-19580172*x3*x4^2-9372130*x4^3-1289811*x1^2+22255617*x1*x2+109167336*x1*x3+150465505*x1*x4-1318549*x2^2+238367004*x2*x3-131012023*x2*x4+124816848*x3^2+116104228*x3*x4-70772721*x4^2+133696391*x1+10206682*x2+145409424*x3-111938286*x4-23665371, 12014836*x1^3-124056256*x1^2*x2+30086207*x1^2*x3-31723130*x1^2*x4+125153935*x1*x2^2+5839912*x1*x2*x3+225831278*x1*x2*x4+32314614*x1*x3^2-11292836*x1*x3*x4-11471604*x1*x4^2-17361686*x2^3-61020343*x2^2*x3-132982128*x2^2*x4+46623254*x2*x3^2-53441750*x2*x3*x4-48142857*x2*x4^2+9644976*x3^3-19580172*x3^2*x4-28116390*x3*x4^2+18646000*x4^3+64756487*x1^2+115474646*x1*x2+150465505*x1*x3-38448834*x1*x4-91874169*x2^2-131012023*x2*x3-176642248*x2*x4+58052114*x3^2-141545442*x3*x4+9369657*x4^2-71771275*x1-58141098*x2-111938286*x3-11893932*x4+7665457]
    gb = decompose(I, interreduction = false)
    @test SG.is_gb(gb)
    @test is_eq(Ideal(R, gb), fancy_loop(I))
    gb = decompose(I)
    @test SG.is_gb(gb)
    @test length(gb) == length(gens(fancy_loop(I, interreduce = true)))
end


text:
using LinearAlgebraicRepresentation
Lar = LinearAlgebraicRepresentation

#Tetrahemihexahedron
"""
Return the Lar.Struct of tetrahemihexahedron
"""
function tetrahemihexahedron()
	V, (VV,EV,FV,CV) = Lar.simplex(3, true)
	tetra = V, EV,FV,CV
	twotetra = Lar.Struct([ tetra, Lar.r(0,π,0),tetra])
	V,EV,FV,CV = Lar.struct2lar(Lar.Struct([twotetra,Lar.r(0,0,π),twotetra]))
	return Lar.Struct([(V,collect(Set(EV)),FV,CV)])
end


text:
"""
    sixel_decode([T=RGB{N0f8}], src, [decoder]; kwargs...) -> img::IndirectArray

Decode the sixel format sequence provided by `src` and output as an indexed image.

# Arguments

- `T`: output eltype. By default it is `RGB{N0f8}`.
- `src`: the input sixel data source. It can be either an `IO`, `AbstractVector{UInt8}`, or `AbstractString`.
- `decoder::AbstractSixelDecoder`: the sixel decoder. Currently only `LibSixelDecoder` is available.

!!! warning
    `sixel_decode` for tiny images (including vectors) is undefined behavior; sixel format requires
    at least `6` pixels in the row direction. For example, you should not expect it returning a
    `Vector` data even if it is.

# Parameters

- `transpose::Bool`: whether we need to permute the image's width and height dimension before encoding.
  The default value is `false`.

# References

- [1] VT330/VT340 Programmer Reference Manual, Volume 1: Text Programming
- [2] VT330/VT340 Programmer Reference Manual, Volume 2: Graphics Programming
- [3] https://github.com/saitoha/libsixel
"""
function sixel_decode end

sixel_decode(src, dec=default_decoder(); kwargs...) = sixel_decode(RGB{N0f8}, src, dec; kwargs...)

sixel_decode(::Type{T}, data::AbstractString, dec=default_decoder(); kwargs...) where T =
    sixel_decode(T, collect(UInt8, convert(String, data)), dec; kwargs...)

function sixel_decode(::Type{T}, bytes::AbstractArray, dec=default_decoder(); transpose=false) where {T}
    bytes = convert(Vector{UInt8}, bytes)

    expected_size = read_sixel_size(bytes)
    index, values = dec(bytes)
    values = eltype(values) == T ? values : map(T, values)

    if dec isa LibSixelDecoder
        # Julia uses column-major order while libsixel uses row-major order,
        # thus transpose=true means no permutation.
        index = transpose ? index : PermutedDimsArray(index, (2, 1))
        expected_size = transpose ? (expected_size[2], expected_size[1]) : expected_size
    else
        throw(ArgumentError("Unsupported decoder type $(typeof(enc)). Please open an issue for this."))
    end

    check_size(size(index), expected_size)
    # We use IndirectArray to mark it as an indexed image so as to avoid unnecessary memory allocation
    # Users that expect a dense Array can always call `collect(rst)` or `convert(Array, rst)` on this.
    return IndirectArray(index, values)
end

function sixel_decode(::Type{T}, io::IO, dec=default_decoder(); kwargs...) where T
    # TODO: This is actually a duplicated check since the bytes method also checks it
    #       but this is a quite fast operation...
    expected_size = read_sixel_size(io)

    bytes = read(io)
    img = sixel_decode(T, bytes, dec; kwargs...)

    check_size(size(img), expected_size)
    return img
end

"""
    read_sixel_size(io::IO)
    read_sixel_size(bytes::Vector{UInt8})

Read the header from a sixel sequence `io`/`bytes` and return the size of the sixel image.
"""
function read_sixel_size(bytes::Vector{UInt8})
    # There's absolutely something wrong if the header content exceeds 50 bytes
    max_header_length = 50

    p_end = findfirst(isequal(UInt8('#')), bytes)
    p_end = isnothing(p_end) ? length(bytes) : p_end - 1
    buffer = view(bytes, 1:min(p_end, max_header_length))
    seps = findall(isequal(UInt8(';')), buffer)
    length(seps) == 3 || throw(ArgumentError("The input data is not recognizable as sixel sequence."))
    w = parse(Int, String(buffer[seps[2]+1:seps[3]-1]))
    h = parse(Int, String(buffer[seps[3]+1:end]))
    # (h, w) in column-major order convention
    return h, w
end

function read_sixel_size(io::IO)
    # Sixel sequence format always start with this
    #     \ePq"1;1;w;h#
    # where `#` indicates the first palette value
    buffer = UInt8[]
    p = position(io)
    try
        # There's absolutely something wrong if the header content exceeds 50 bytes
        i, max_header_length = 1, 50
        ch = read(io, Cuchar)
        # Cuchar('#') == 0x23
        while ch != 0x23 && i < max_header_length
            push!(buffer, ch)
            ch = read(io, Cuchar)
            i += 1
        end
    catch e
        e isa EOFError && throw(ArgumentError("The input data is not recognizable as sixel sequence."))
        rethrow(e)
    finally
        seek(io, p)
    end
    return read_sixel_size(buffer)
end

function check_size(actual_size, expected_size)
    valid = length(actual_size)==length(expected_size)
    valid = valid && all(zip(actual_size, expected_size)) do x
        # libsixel encode/decode for small images is undefined behavior to us
        minimum(x) <= 6 && return true
        x[1] == x[2]
    end
    if !valid
        @warn "Output size mismatch during decoding sixel sequences" actual_size expected_size
    end
end