text:
struct ConstantFrequency end
struct ReferenceBus end

function _get_frequency_state(d::DynamicWrapper{T}) where {T <: PSY.DynamicGenerator}
    return get_global_index(d)[:ω]
end

function _get_frequency_state(d::DynamicWrapper{T}) where {T <: PSY.DynamicInverter}
    return get_global_index(d)[:ω_oc]
end

function _get_frequency_state(d::DynamicWrapper{PSY.PeriodicVariableSource})
    return 0
end

function get_frequency_reference!(
    ::Type{ConstantFrequency},
    dynamic_injection_data::Vector,
    static_injection_data::Vector,
)
    ref_devices = filter(x -> get_bus_category(x) == SLACKBus, static_injection_data)
    if length(ref_devices) < 1
        dyn_devices = filter(x -> get_bus_category(x) == SLACKBus, dynamic_injection_data)
        if length(dyn_devices) < 1
            throw(
                IS.ConflictingInputsError(
                    "ConstantFrequency model requires at least one generation unit at the Bus BusTypes.REF",
                ),
            )
        end
    end
    return 0
end

function get_frequency_reference!(
    ::Type{ReferenceBus},
    dynamic_injection_data::Vector,
    static_injection_data::Vector,
)
    reference = -1
    static_ref_devices = filter(x -> get_bus_category(x) == SLACKBus, static_injection_data)
    if isempty(static_ref_devices)
        ref_devices = filter(x -> get_bus_category(x) == SLACKBus, dynamic_injection_data)
        if length(ref_devices) > 1
            ref_devices = filter(
                x -> (
                    haskey(get_global_index(x), :ω_oc) || haskey(get_global_index(x), :ω)
                ),
                ref_devices,
            )
            if isempty(ref_devices)
                throw(
                    IS.ConflictingInputsError(
                        "ReferenceBus model requires at least one generator with frequency controls.",
                    ),
                )
            end
            ref_devices = sort(ref_devices, by = x -> PSY.get_base_power(x))
        elseif length(ref_devices) < 1
            throw(
                IS.ConflictingInputsError(
                    "ReferenceBus model requires at least one bus of type BusTypes.REF with a DynamicInjection or Source device connected to it",
                ),
            )
        elseif length(ref_devices) != 1
            @assert false
        end
        reference = _get_frequency_state(first(ref_devices))
    else
        @info "The reference Bus has a Source connected to it. The frequency reference model will change to ConstantFrequency"
        reference = 0
    end
    return reference
end


text:
export
  getCliqueStatus,
  setCliqueStatus!

export
  stackCliqUpMsgsByVariable,
  getCliqDownMsgsAfterDownSolve

# likely to be deleted at some point

## =============================================================================
## Clique status accessors
## =============================================================================

"""
    $SIGNATURES

Return `::Symbol` status a particular clique is in, with specific regard to solution
or numerical initialization status:
- :needdownmsg
- UPSOLVED
- DOWNSOLVED
- INITIALIZED
- MARGINALIZED
- NULL

Notes:
- `NULL` represents the first uninitialized state of a cliq.
"""
getCliqueStatus(cliqdata::BayesTreeNodeData) = cliqdata.status
getCliqueStatus(cliq::TreeClique) = getCliqueStatus(getCliqueData(cliq))

"""
    $SIGNATURES

Set up initialization or solve status of this `cliq`.
"""
function setCliqueStatus!(cdat::BayesTreeNodeData, status::CliqStatus)
  cdat.status = status
end
setCliqueStatus!(cliq::TreeClique, status::CliqStatus) = setCliqueStatus!(getCliqueData(cliq), status)

## =============================================================================
## Regular up and down Message Registers/Channels, getters and setters
## =============================================================================

## =============================================================================
## Message channel put/take! + buffer message accessors
## =============================================================================

## -----------------------------------------------------------------------------
## UP
## -----------------------------------------------------------------------------
"""
$SIGNATURES
Get the message channel
"""
getMsgUpChannel(tree::MetaBayesTree, edge) = MetaGraphs.get_prop(tree.bt, edge, :upMsg)

"""
$SIGNATURES

Put a belief message on the up tree message channel `edge`. Blocks until a take! is performed by a different task.
"""
function putBeliefMessageUp!(tree::AbstractBayesTree, edge, beliefMsg::LikelihoodMessage)
  # Blocks until data is available.
  put!(getMsgUpChannel(tree, edge), beliefMsg)
  return beliefMsg
end

"""
$SIGNATURES

Remove and return belief message from the up tree message channel edge. Blocks until data is available.
"""
function takeBeliefMessageUp!(tree::AbstractBayesTree, edge)
  # Blocks until data is available.
  beliefMsg = take!(getMsgUpChannel(tree, edge))
  return beliefMsg
end

## -----------------------------------------------------------------------------
## DOWN
## -----------------------------------------------------------------------------
"""
$SIGNATURES
Get the message channel
"""
getMsgDwnChannel(tree::MetaBayesTree, edge) = MetaGraphs.get_prop(tree.bt, edge, :downMsg)

"""
    $SIGNATURES

Put a belief message on the down tree message channel edge. Blocks until a take! is performed by a different task.
"""
function putBeliefMessageDown!(tree::AbstractBayesTree, edge, beliefMsg::LikelihoodMessage)
  # Blocks until data is available.
  put!(getMsgDwnChannel(tree, edge), beliefMsg)
  return beliefMsg
end

"""
    $SIGNATURES

Remove and return a belief message from the down tree message channel edge. Blocks until data is available.
"""
function takeBeliefMessageDown!(tree::AbstractBayesTree, edge)
  # Blocks until data is available.
  beliefMsg = take!(getMsgDwnChannel(tree, edge))
  return beliefMsg
end

##==============================================================================
## Clique Message Buffers
##==============================================================================
function getMessageBuffer(btnd::BayesTreeNodeData)
  btnd.messages
end
getMessageBuffer(clique::TreeClique) = getCliqueData(clique).messages

# getMessageUpRx(clique::TreeClique) = getMessageBuffer(clique).upRx
# getMessageDownRx(clique::TreeClique) = getMessageBuffer(clique).downRx

# getMessageUpTx(clique::TreeClique) = getMessageBuffer(clique).upTx
# getMessageDownTx(clique::TreeClique) = getMessageBuffer(clique).downTx

ext: jl


text:
# Copyright (c) Microsoft Corporation.
# Licensed under the MIT license.

#Pkg.add("ReverseDiffSource")
#Pkg.update()
#Pkg.checkout("ReverseDiffSource","devl")
#Pkg.status()
using ReverseDiffSource

typealias Mat Matrix{Float64}

function function2expression(f::Function, sig::Tuple; order::Int=1, evalmod=Main, debug=false, allorders=true)
    fs = methods(f, map( typeof, sig ))
    length(fs) == 0 && error("ardiff: no function '$f' found for signature $sig")
    length(fs) > 1  && error("ardiff: several functions $f found for signature $sig")  # is that possible ?

    fdef  = fs[1].func.code
    fcode = Base.uncompressed_ast(fdef)
    fargs = fcode.args[1]  # function parameters

    cargs = [ (fargs[i], sig[i]) for i in 1:length(sig) ]
    expr = ReverseDiffSource.streamline(fcode.args[3])
    expr,cargs
end

function logsumexp(x::Mat)
  mx = maximum(x)
  log(sum(exp(x .- mx))) + mx
end

gmm_objective_code =
  :(begin
      n = size(x,2)
      d = size(x,1)
      CONSTANT = -n*d*0.5*log(2 * pi)

      sum_qs = similar(alphas)
      for ik in 1:k
        sum_qs[ik] = sum(icf[1:d,ik])
      end

      main_term = similar(alphas)
      slse = 0.
      for ix in 1:n
        for ik in 1:k
          Qxcentered = get_Q(d,icf[:,ik])*(x[:,ix] - means[:,ik])
          main_term[ik] = -0.5*sum(Qxcentered.^2)
        end
        slse += logsumexp(alphas + sum_qs + main_term)
      end
      return CONSTANT + slse - n*logsumexp(alphas)
  end)

k=3
d=2
n=3
icf_sz=div(d*(d + 1),2)
alphas = randn(1,k)
means = randn(d,k)
icf = randn(icf_sz,k)
x = randn(d,n)

logsumexp_code,cargs = function2expression(logsumexp, (alphas,))
logsumexp_d_code = rdiff(logsumexp_code; x=alphas, allorders=false)
@eval logsumexp_d(x) = $logsumexp_d_code
@deriv_rule logsumexp(x) x logsumexp_d(x)
@deriv_rule similar(x) x zeros(size(x))

function make_row(r::Int, diag_elem, L)
      elems=zeros(eltype(L),1,r-1)
      for i=1:r-1
        elems[i]=L[i]
      end
      return [elems diag_elem zeros(d-r)']
end

@deriv_rule make_row(r,diag_elem,L) r 0.
@deriv_rule make_row(r,diag_elem,L) diag_elem make_row(r,0.,ones(L)).*ds
@deriv_rule make_row(r,diag_elem,L) L make_row(r,1.,zeros(L)).*ds
make_row_code_d_code = rdiff(make_row,(1,2,icf[d+1:end,:]))

function ltri_unpack(D,LT)
  d=length(D)
  sum(D)+sum(LT)
  function make_row(r::Int, L)
    elems=zeros(eltype(L),1,r-1)
    for i=1:r-1
      elems[i]=L[i]
    end
    hcat(elems, D[r], zeros(1,d-r))
  end
  #row_start(r::Int) = div((r-1)*(r-2),2)
  #inds(r) = row_start(r)+(1:r-1)
  #vcat([ make_row(r, LT[inds(r)]) for r=1:d ]...)
end
ltri_unpack_code,cargs = function2expression(ltri_unpack, (icf[1:d,1],icf[d+1:end,1]))
ltri_unpack_d_code = rdiff(ltri_unpack_code; D=icf[1:d,1], LT=icf[d+1:end,1])
@eval logsumexp_d(x) = $ltri_unpack_d_code
@deriv_rule ltri_unpack(x) x ltri_unpack_d(x)

#gmm_objective_AD_code,cargs = function2expression(gmm_objective_AD, (alphas,x))
gmm_objective_AD_d_code = rdiff(gmm_objective_code;alphas=alphas,means=means,icf=icf)
@eval g(alphas,means,icf) = $gmm_objective_AD_d_code
g(alphas,means,icf)

include("common.jl")

# input should be 1 dimensional
function logsumexp(x)
  mx = maximum(x)
  log(sum(exp(x - mx))) + mx
end

@inbounds function gmm_objective(alphas,means,icf,x,wishart::Wishart)
  d = size(x,1)
  n = size(x,2)
  CONSTANT = -n*d*0.5*log(2 * pi)

  sum_qs = sum(icf[1:d,:],1)
  slse = sum(sum_qs)
  Qs = [get_Q(d,icf[:,ik]) for ik in 1:k]
  main_term = zeros(eltype(alphas),1,k)

  slse = 0.
  for ix=1:n
    for ik=1:k
      main_term[ik] = -0.5*sumabs2(Qs[ik] * (x[:,ix] - means[:,ik]))
    end
    slse += logsumexp(alphas + sum_qs + main_term)
  end

  CONSTANT + slse - n*logsumexp(alphas) + log_wishart_prior(wishart, sum_qs, Qs, icf)
end

# Read instance
dir_in = ARGS[1]
dir_out = ARGS[2]
fn = ARGS[3]
nruns_f = parse(Int,ARGS[4])
nruns_J = parse(Int,ARGS[5])
replicate_point = size(ARGS,1) >= 6 && ARGS[6] == "-rep"

fn_in = string(dir_in, fn)
fn_out = string(dir_out, fn)

alphas,means,icf,x,wishart = read_gmm_instance(string(fn_in,".txt"),replicate_point)
d = size(means,1)
k = size(means,2)
n = size(x,2)

# Objective
precompile(gmm_objective,(typeof(alphas),typeof(means),typeof(icf),typeof(x),Wishart))
tic()
for i in 1:nruns_f
  err = gmm_objective(alphas,means,icf,x,wishart)
end
tf = toq()/nruns_f;
@printf "tf: %g\n" tf
#@printf "err: %f\n" err

# Gradient helper
function wrapper_gmm_objective(packed)
  alphas,means,icf = unpack(d,k,packed)
  gmm_objective(alphas,means,icf,x,wishart)
end

# Gradient
g = ForwardDiff.gradient(wrapper_gmm_objective)
precompile(pack,(typeof(alphas),typeof(means),typeof(icf)))
precompile(g,(typeof(alphas),typeof(means),typeof(icf)))
J = zeros(1,length(alphas)+length(means)+length(icf))
tic()
for i in 1:nruns_J
  J = g(pack(alphas,means,icf))
end
tJ = toq()/nruns_J;
@printf "tJ: %g\n" tJ
#println("J:")
#println(J)

name = "Julia_F"

write_J(string(fn_out,"_J_",name,".txt"),J)
write_times(string(fn_out,"_times_",name,".txt"),tf,tJ)

max_stars_repo_name: dsyme/ADBench

max_issues_repo_name: dsyme/ADBench

max_forks_repo_name: dsyme/ADBench


text:
#
# This file is a part of ANISPROU. License is MIT
# Copyright (c) 2020 Andrea Zunino
#

"""
ANISPROU

Analysis of isothermal titration calorimetry (ITC) data on sodium dodecyl sulphate (SDS) mediated protein unfolding.

# Exports

$(EXPORTS)
"""
module ANISPROU

using DelimitedFiles
using Distributions
using ForwardDiff
#using ReverseDiff
#using HDF5
using JLD2
using LinearAlgebra
using Roots
using Cubature,QuadGK
using Optim
using PyPlot
PyPlot.rc("font", family="serif", size=12)
using WriteVTK
using DocStringExtensions
using Requires
using GLMakie

export ScaledBeta2DParams,BetaMix2D,ITCObsData
export readallexperiments
export solveinvprob
export plotobsdata,plotinitialguess,plotresults
export plotsingleexperiments
export saveresultVTK
export findcurvefeatures
export calcfreeSDSNbound
export area_enthalpy,volume_enthalpy,areasvsprotcon
export plotbindisotherm,plotfoundfeatures
export plotareavsprotcon,plotbetacomp1D
export plotparamlines
export plotsurface3D

include("utils.jl")
include("readITCdata.jl")
include("invertITCdata2D.jl")
include("betaintegrals.jl")
include("bindingisotherm.jl")
include("plotstuff.jl")

# 3D plotting
include("plot3d.jl")

# # 3D plotting
# function __init__()
#     @require GLMakie="e9467ef8-e4e7-5192-8a1a-b1aee30e663a"  begin
#         # @require AbstractPlotting="537997a7-5e4e-5d89-9595-2241ea00577e" begin
#         import .GLMakie
#         include("plot3d.jl")
#         export plotsurface3D
#        # end
#     end
# end

end # module


text:
__precompile__(true)
module uCSV

using DataFrames
using Dates

COLMAP{T} = Union{Dict{String, T}, Dict{Int, T}}

include("read.jl")
include("write.jl")
include("parsesource.jl")
include("parsefields.jl")
include("checkfield.jl")
include("helperfunctions.jl")

end # module


text:
@cache struct Feagin10Cache{uType,uNoUnitsType,rateType,TabType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  fsalfirst::rateType
  k2::rateType
  k3::rateType
  k4::rateType
  k5::rateType
  k6::rateType
  k7::rateType
  k8::rateType
  k9::rateType
  k10::rateType
  k11::rateType
  k12::rateType
  k13::rateType
  k14::rateType
  k15::rateType
  k16::rateType
  k17::rateType
  tmp::uType
  atmp::uNoUnitsType
  k::rateType
  tab::TabType
end

function alg_cache(alg::Feagin10,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Val{true})
  tab = Feagin10ConstantCache(constvalue(uBottomEltypeNoUnits),constvalue(tTypeNoUnits))
  k1 = zero(rate_prototype); k2 = zero(rate_prototype); k3 = zero(rate_prototype); k4 = zero(rate_prototype); k5 = zero(rate_prototype)
  k6 = zero(rate_prototype); k7 = zero(rate_prototype); k8 = zero(rate_prototype); k9 = zero(rate_prototype); k10 = zero(rate_prototype)
  k11 = zero(rate_prototype); k12 = zero(rate_prototype); k13 = zero(rate_prototype); k14 = zero(rate_prototype)
  k15 = zero(rate_prototype); k16 = zero(rate_prototype); k17 = zero(rate_prototype)
  tmp = zero(u); atmp = similar(u,uEltypeNoUnits); k = zero(rate_prototype)

  Feagin10Cache(u,uprev,k1,k2,k3,k4,k5,k6,k7,k8,k9,k10,k11,k12,k13,k14,k15,k16,k17,tmp,atmp,k,tab)
end

alg_cache(alg::Feagin10,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Val{false}) = Feagin10ConstantCache(constvalue(uBottomEltypeNoUnits),constvalue(tTypeNoUnits))

@cache struct Feagin12Cache{uType,uNoUnitsType,rateType,TabType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  fsalfirst::rateType
  k2::rateType
  k3::rateType
  k4::rateType
  k5::rateType
  k6::rateType
  k7::rateType
  k8::rateType
  k9::rateType
  k10::rateType
  k11::rateType
  k12::rateType
  k13::rateType
  k14::rateType
  k15::rateType
  k16::rateType
  k17::rateType
  k18::rateType
  k19::rateType
  k20::rateType
  k21::rateType
  k22::rateType
  k23::rateType
  k24::rateType
  k25::rateType
  tmp::uType
  atmp::uNoUnitsType
  k::rateType
  tab::TabType
end

function alg_cache(alg::Feagin12,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Val{true})
  tab = Feagin12ConstantCache(constvalue(uBottomEltypeNoUnits),constvalue(tTypeNoUnits))
  k1 = zero(rate_prototype); k2 = zero(rate_prototype); k3 = zero(rate_prototype); k4 = zero(rate_prototype); k5 = zero(rate_prototype)
  k6 = zero(rate_prototype); k7 = zero(rate_prototype); k8 = zero(rate_prototype); k9 = zero(rate_prototype); k10 = zero(rate_prototype)
  k11 = zero(rate_prototype); k12 = zero(rate_prototype); k13 = zero(rate_prototype); k14 = zero(rate_prototype)
  k15 = zero(rate_prototype); k16 = zero(rate_prototype); k17 = zero(rate_prototype); k18 = zero(rate_prototype)
  k19 = zero(rate_prototype); k20 = zero(rate_prototype); k21 = zero(rate_prototype); k22 = zero(rate_prototype)
  k23 = zero(rate_prototype); k24 = zero(rate_prototype); k25 = zero(rate_prototype)
  tmp = zero(u); atmp = similar(u,uEltypeNoUnits); k = zero(rate_prototype)

  Feagin12Cache(u,uprev,k1,k2,k3,k4,k5,k6,k7,k8,k9,k10,k11,k12,k13,k14,k15,k16,k17,k18,k19,k20,k21,k22,k23,k24,k25,tmp,atmp,k,tab)
end

alg_cache(alg::Feagin12,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Val{false}) = Feagin12ConstantCache(constvalue(uBottomEltypeNoUnits),constvalue(tTypeNoUnits))

@cache struct Feagin14Cache{uType,uNoUnitsType,rateType,TabType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  fsalfirst::rateType
  k2::rateType
  k3::rateType
  k4::rateType
  k5::rateType
  k6::rateType
  k7::rateType
  k8::rateType
  k9::rateType
  k10::rateType
  k11::rateType
  k12::rateType
  k13::rateType
  k14::rateType
  k15::rateType
  k16::rateType
  k17::rateType
  k18::rateType
  k19::rateType
  k20::rateType
  k21::rateType
  k22::rateType
  k23::rateType
  k24::rateType
  k25::rateType
  k26::rateType
  k27::rateType
  k28::rateType
  k29::rateType
  k30::rateType
  k31::rateType
  k32::rateType
  k33::rateType
  k34::rateType
  k35::rateType
  tmp::uType
  atmp::uNoUnitsType
  k::rateType
  tab::TabType
end

function alg_cache(alg::Feagin14,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Val{true})
  tab = Feagin14ConstantCache(constvalue(uBottomEltypeNoUnits),constvalue(tTypeNoUnits))
  k1 = zero(rate_prototype); k2 = zero(rate_prototype); k3 = zero(rate_prototype); k4 = zero(rate_prototype); k5 = zero(rate_prototype)
  k6 = zero(rate_prototype); k7 = zero(rate_prototype); k8 = zero(rate_prototype); k9 = zero(rate_prototype); k10 = zero(rate_prototype)
  k11 = zero(rate_prototype); k12 = zero(rate_prototype); k13 = zero(rate_prototype); k14 = zero(rate_prototype)
  k15 = zero(rate_prototype); k16 = zero(rate_prototype); k17 = zero(rate_prototype); k18 = zero(rate_prototype)
  k19 = zero(rate_prototype); k20 = zero(rate_prototype); k21 = zero(rate_prototype); k22 = zero(rate_prototype)
  k23 = zero(rate_prototype); k24 = zero(rate_prototype); k25 = zero(rate_prototype)
  k26 = zero(rate_prototype); k27 = zero(rate_prototype); k28 = zero(rate_prototype)
  k29 = zero(rate_prototype); k30 = zero(rate_prototype); k31 = zero(rate_prototype); k32 = zero(rate_prototype)
  k33 = zero(rate_prototype); k34 = zero(rate_prototype); k35 = zero(rate_prototype)
  tmp = zero(u); atmp = similar(u,uEltypeNoUnits); k = zero(rate_prototype)

  Feagin14Cache(u,uprev,k1,k2,k3,k4,k5,k6,k7,k8,k9,k10,k11,k12,k13,k14,k15,k16,
                k17,k18,k19,k20,k21,k22,k23,k24,k25,k26,k27,k28,k29,k30,
                k31,k32,k33,k34,k35,tmp,atmp,k,tab)
end

alg_cache(alg::Feagin14,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Val{false}) = Feagin14ConstantCache(constvalue(uBottomEltypeNoUnits),constvalue(tTypeNoUnits))


text:
@generated function promote_rule(::Type{V1}, ::Type{V2}) where {V1 <: Variable, V2 <: Variable}
    if name(V1) < name(V2)
        :(Monomial{(V1(), V2())})
    else
        :(Monomial{(V2(), V1())})
    end
end

promote_rule(::Type{M}, ::Type{V}) where {V <: Variable, M <: Monomial} = promote_rule(Monomial{(V(),), (1,)}, M)
promote_rule(V::Type{<:Variable}, M::Type{<:Monomial}) = promote_rule(M, V)

function _promote_monomial(::Type{<:Monomial{V1}}, ::Type{<:Monomial{V2}}) where {V1, V2}
    varnames = union(name.(V1), name.(V2))
    sort!(varnames)
    vars = Expr(:tuple, [Expr(:call, Expr(:curly, :Variable, Expr(:quote, n))) for n in varnames]...)
    :(Monomial{$vars})
end

@generated function promote_rule(::Type{M1}, ::Type{M2}) where {M1 <: Monomial, M2 <: Monomial}
    _promote_monomial(M1, M2)
end

promote_rule(::Type{Term{T, M}}, ::Type{V}) where {V <: Variable, T, M <: Monomial} = Term{T, promote_type(V, M)}
promote_rule(V::Type{<:Variable}, T::Type{<:Term}) = promote_rule(T, V)

function promote_rule(::Type{Term{T, M2}}, ::Type{M1}) where {M1 <: Monomial, T, M2 <: Monomial}
    Term{T, promote_type(M1, M2)}
end
promote_rule(M::Type{<:Monomial}, T::Type{<:Term}) = promote_rule(T, M)

function _promote_term(::Type{Term{T1, M1}}, ::Type{Term{T2, M2}}) where {T1, M1, T2, M2}
    Term{promote_type(T1, T2), promote_type(M1, M2)}
end

promote_rule(T1::Type{<:Term}, T2::Type{<:Term}) = _promote_term(T1, T2)

function promote_rule(::Type{V}, ::Type{S}) where {S, V <: Variable}
    Term{S, Monomial{(V(),)}}
end

function promote_rule(::Type{M}, ::Type{S}) where {S, M <: Monomial}
    Term{S, M}
end

function _promote_any_term(::Type{S}, ::Type{Term{T, M}}) where {S, T, M <: Monomial}
    Term{promote_type(S, T), M}
end

promote_rule(t::Type{<:Term}, s) = _promote_any_term(s, t)


text:
using NNlib, Test, Statistics, Random
using ChainRulesCore, ChainRulesTestUtils
using Base.Broadcast: broadcasted
import FiniteDifferences
import Zygote
using Zygote: gradient
using StableRNGs
using CUDA

if VERSION < v"1.6"
    @info "skipping doctests, on Julia $VERSION"
else
    using Documenter
    DocMeta.setdocmeta!(NNlib, :DocTestSetup, :(using NNlib, UnicodePlots); recursive=true)
    @testset "Doctests" begin
        doctest(NNlib, manual=false)
    end
end

const rng = StableRNG(123)

include("test_utils.jl")

@testset "Activation Functions" begin
    include("activations.jl")
end

@testset "Batched Multiplication" begin
    include("batchedmul.jl")
end

@testset "Convolution" begin
    include("conv.jl")
    include("conv_bias_act.jl")
end

@testset "Inference" begin
    include("inference.jl")
end

@testset "Pooling" begin
    include("pooling.jl")
end

@testset "Padding" begin
    include("padding.jl")
end

@testset "Softmax" begin
    include("softmax.jl")
end

@testset "Upsampling" begin
    include("upsample.jl")
end

@testset "Gather" begin
    include("gather.jl")
end

@testset "Scatter" begin
    include("scatter.jl")
end

@testset "Utilities" begin
    include("utils.jl")
end

@testset "Grid Sampling" begin
    include("sampling.jl")
end

@testset "Functions" begin
    include("functions.jl")
end

if VERSION >= v"1.6" && CUDA.functional()
    if get(ENV, "NNLIB_TEST_CUDA", "false") == "true"
        import Pkg
        using NNlibCUDA
        @testset "CUDA" begin
            Pkg.test("NNlibCUDA")
        end
    else
        @info "Skipping CUDA tests, set NNLIB_TEST_CUDA=true to run them"
    end
else
    @info "Insufficient version or CUDA not found; Skipping CUDA tests"
end


text:
struct Response{S<:IO}
    response::HTTP.Response
    io::S  # I/O type must support `seek`
    dict::Ref{AbstractDict}  # Caches the parsed dict response
end

Response(response::HTTP.Response, io::IO) = Response(response, io, Ref{AbstractDict}())

function Base.getproperty(r::Response, f::Symbol)
    if f === :io || f === :response || f === :dict
        # Direct access to the fields stored in `Response`
        return getfield(r, f)
    elseif f === :body
        # As we're streaming the requests we'll fake the `HTTP.Response` body field
        return _rewind(read, r.io)::Vector{UInt8}
    else
        # Pretend like we're an `HTTP.Response` type for any other field access
        return getproperty(r.response, f)
    end
end

function mime_type(r::Response)
    # Parse response data according to mimetype...
    # Using "application/octet-stream" as the fallback MIME type as recommended by:
    # https://developer.mozilla.org/en-US/docs/Web/HTTP/Basics_of_HTTP/MIME_types/Common_types
    mime = HTTP.header(r.response, "Content-Type", "application/octet-stream")

    # When the MIME type is not specified attempt to determine type from first bytes of the
    # stream.
    magic_bytes = if isempty(HTTP.header(r.response, "Content-Type", ""))
        _rewind(r.io) do io
            read(io, 5)
        end
    else
        UInt8[]
    end

    T = if occursin(r"/xml", mime) || magic_bytes == b"<?xml"
        MIME"application/xml"
    elseif occursin(r"/x-amz-json-1\.[01]$", mime) || endswith(mime, "json")
        MIME"application/json"
    elseif startswith(mime, "text/")
        MIME"text/plain"
    else
        MIME{Symbol(mime)}  # Unhandled MIME type
    end

    return T
end

function Base.parse(f::Function, r::Response, mime::MIME=mime_type(r)())
    result = _rewind(r.io) do io
        f(io, mime)
    end
    return result
end

Base.parse(r::Response, mime::MIME=mime_type(r)()) = parse(_read, r, mime)

function _read(io::IO, ::MIME"application/xml")
    xml = parse_xml(read(io, String))
    root = XMLDict.root(xml.x)  # Drop XML declaration
    return xml_dict(root, LittleDict{Union{String,Symbol},Any})
end

function _read(io::IO, ::MIME"application/json")
    # Note: Using JSON instead of JSON3 since it does not support OrderedDict/LittleDict
    return JSON.parse(io; dicttype=LittleDict{String,Any})
end

_read(io::IO, ::MIME"text/plain") = read(io, String)
_read(io::IO, ::MIME) = read(io)

# Dict-like access
function _dict(r::Response)
    if !isassigned(r.dict)
        dict = parse(r)::AbstractDict
        r.dict[] = dict
    end

    return r.dict[]
end

Base.getindex(r::Response, key::Union{AbstractString,Symbol}) = getindex(_dict(r), key)
Base.haskey(r::Response, key::Union{AbstractString,Symbol}) = haskey(_dict(r), key)
Base.keys(r::Response) = keys(_dict(r))
Base.values(r::Response) = values(_dict(r))

function Base.iterate(r::Response)
    iter = parse(r)
    x = iterate(iter)
    x === nothing && return nothing
    el, s = x
    return (el, (iter, s))
end

function Base.iterate(r::Response, state)
    iter, s = state
    x = iterate(iter, s)
    x === nothing && return nothing
    el, s = x
    return (el, (iter, s))
end

Base.String(r::Response) = _rewind(io -> read(io, String), r.io)

function Base.show(io::IO, m::MIME"text/plain", r::Response)
    println(io, "$(Response): $(mime_type(r)()) interpreted as:")

    content = parse(r)
    show(io, m, content)

    return nothing
end

function _rewind(f::Function, io::IO)
    pos = position(io)
    seekstart(io)
    result = try
        f(io)
    finally
        seek(io, pos)
    end
    return result
end


text:
__precompile__()

module RowEchelon
using LinearAlgebra

export rref, rref!

function rref!(A::Matrix{T}, ɛ=T <: Union{Rational,Integer} ? 0 : eps(norm(A,Inf))) where T
    nr, nc = size(A)
    i = j = 1
    while i <= nr && j <= nc
        (m, mi) = findmax(abs.(A[i:nr,j]))
        mi = mi+i - 1
        if m <= ɛ
            if ɛ > 0
                A[i:nr,j] .= zero(T)
            end
            j += 1
        else
            for k=j:nc
                A[i, k], A[mi, k] = A[mi, k], A[i, k]
            end
            d = A[i,j]
            for k = j:nc
                A[i,k] /= d
            end
            for k = 1:nr
                if k != i
                    d = A[k,j]
                    for l = j:nc
                        A[k,l] -= d*A[i,l]
                    end
                end
            end
            i += 1
            j += 1
        end
    end
    A
end

rrefconv(::Type{T}, A::Matrix) where {T} = rref!(copyto!(similar(A, T), A))

"""
    rref(A)
Compute the reduced row echelon form of the matrix A.
Since this algorithm is sensitive to numerical imprecision,
* Complex numbers are converted to ComplexF64
* Integer, Float16 and Float32 numbers are converted to Float64
* Rational are kept unchanged

```jldoctest
julia> rref([ 1  2 -1  -4;
              2  3 -1 -11;
             -2  0 -3  22])
3×4 Array{Float64,2}:
 1.0  0.0  0.0  -8.0
 0.0  1.0  0.0   1.0
 0.0  0.0  1.0  -2.0

julia> rref([16  2  3  13;
              5 11 10   8;
              9  7  6  12;
              4 14 15   1])
4×4 Array{Float64,2}:
 1.0  0.0  0.0   1.0
 0.0  1.0  0.0   3.0
 0.0  0.0  1.0  -3.0
 0.0  0.0  0.0   0.0

julia> rref([ 1  2  0   3;
              2  4  0   7])
2×4 Array{Float64,2}:
 1.0  2.0  0.0  0.0
 0.0  0.0  0.0  1.0
```
"""
rref(A::Matrix{T}) where {T} = rref!(copy(A))
rref(A::Matrix{T}) where {T <: Complex} = rrefconv(ComplexF64, A)
rref(A::Matrix{ComplexF64}) = rref!(copy(A))
rref(A::Matrix{T}) where {T <: Union{Integer, Float16, Float32}} = rrefconv(Float64, A)
rref(A::AbstractMatrix) = rref(Matrix(A))
rref(A::Adjoint) = rref(Matrix(A))

include("RowEchelon_with_pivots.jl")

end # module


text:
using DifferentialEquations, ModelingToolkit
using HoJMTK.SimpleFluid

R = 1.0
C = 1.0
P = 1.0
@named resistor = FluidResistor(R=R)
@named compliance = VesselCompliance(C=C)
@named source = Pump(P=P)
@named ground = FluidGround()

frc_eqs = [
          connect(source.i, resistor.i)
          connect(resistor.o, compliance.i)
          connect(compliance.o, source.o, ground.g)
         ]

@variables t
@named _frc_model = ODESystem(frc_eqs, t)
@named frc_model = compose(_frc_model, [resistor, compliance, source, ground])
sys = structural_simplify(frc_model)
u0 = [
        compliance.p => 0.0
     ]
prob = ODAEProblem(sys, u0, (0, 10.0))
sol = solve(prob, Tsit5())

using Plots
plot(sol)


text:
run(`$(joinpath(JULIA_HOME, "julia")) $(joinpath(dirname(@__FILE__), "real.jl"))`)
run(`$(joinpath(JULIA_HOME, "julia")) $(joinpath(dirname(@__FILE__), "complex.jl"))`)


text:
using Adapt
using CUDA: CuArray
using Oceananigans.Fields: ReducedField, fill_halo_regions!
using Oceananigans.Architectures: arch_array

abstract type AbstractGridFittedBoundary <: AbstractImmersedBoundary end

#####
##### GridFittedBoundary
#####

struct GridFittedBoundary{S} <: AbstractGridFittedBoundary
    mask :: S
end

@inline is_immersed(i, j, k, underlying_grid, ib::GridFittedBoundary) = ib.mask(node(c, c, c, i, j, k, underlying_grid)...)

#####
##### GridFittedBottom
#####

"""
    GridFittedBottom(bottom)

Return an immersed boundary...
"""
struct GridFittedBottom{B} <: AbstractGridFittedBoundary
    bottom :: B
end

@inline function is_immersed(i, j, k, underlying_grid, ib::GridFittedBottom)
    x, y, z = node(c, c, c, i, j, k, underlying_grid)
    return z < ib.bottom(x, y)
end

@inline function is_immersed(i, j, k, underlying_grid, ib::GridFittedBottom{<:AbstractArray})
    x, y, z = node(c, c, c, i, j, k, underlying_grid)
    return @inbounds z < ib.bottom[i, j]
end

const ArrayGridFittedBottom = GridFittedBottom{<:Array}
const CuArrayGridFittedBottom = GridFittedBottom{<:CuArray}

function ImmersedBoundaryGrid(grid, ib::Union{ArrayGridFittedBottom, CuArrayGridFittedBottom})
    # Wrap bathymetry in an OffsetArray with halos
    arch = grid.architecture
    bottom_field = ReducedField(Center, Center, Nothing, arch, grid; dims=3)
    bottom_data = arch_array(arch, ib.bottom)
    bottom_field .= bottom_data
    fill_halo_regions!(bottom_field, arch)
    offset_bottom_array = dropdims(bottom_field.data, dims=3)
    new_ib = GridFittedBottom(offset_bottom_array)
    return ImmersedBoundaryGrid(grid, new_ib)
end

const GFBIBG = ImmersedBoundaryGrid{<:Any, <:Any, <:Any, <:Any, <:Any, <:GridFittedBottom}

@inline Δzᵃᵃᶜ(i, j, k, ibg::GFBIBG) = ifelse(is_immersed(i, j, k, ibg.grid, ibg.immersed_boundary),
                                             zero(eltype(ibg.grid)),
                                             Δzᵃᵃᶜ(i, j, k, ibg.grid))

@inline Δzᶠᶜᶜ(i, j, k, ibg::GFBIBG) =  ifelse(solid_node(Face(), Center(), Center(), i  , j, k, ibg),
                                             zero(eltype(ibg)),
                                             Δzᵃᵃᶜ(i, j, k, ibg.grid))

@inline Δzᶜᶠᶜ(i, j, k, ibg::GFBIBG) = ifelse(solid_node(Center(), Face(), Center(), i  , j, k, ibg),
                                             zero(eltype(ibg)),
                                             Δzᵃᵃᶜ(i, j, k, ibg.grid))

@inline Δzᵃᵃᶠ(i, j, k, ibg::GFBIBG) = ifelse(is_immersed(i, j, k, ibg.grid, ibg.immersed_boundary),
                                             zero(eltype(ibg.grid)),
                                             Δzᵃᵃᶠ(i, j, k, ibg.grid))

Adapt.adapt_structure(to, ib::GridFittedBottom) = GridFittedBottom(adapt(to, ib.bottom))


text:
#=
old_eval:
- Julia version: 1.3.1
- Author: Berk
- Date: 2020-11-05
=#

function evaluate(constraint::JuMP.ConstraintRef, data::Union{Dict, DataFrame})
    constr_obj = constraint_object(constraint)
    rhs_const = get_constant(constr_obj.set)
    clean_data = data_to_DataFrame(data);
    vals = [JuMP.value(constr_obj.func, i -> get(clean_data, string(i), Inf)[j]) for j=1:size(clean_data,1)]
    if any(isinf.(vals))
        throw(OCTException(string("Constraint ", constraint, " returned an infinite value.
                                   Make sure you have declared all bbl.vars!")))
    end
    if isnothing(rhs_const)
        vals = [-1*distance_to_set(vals[i], constr_obj.set) for i=1:length(vals)]
    else
        coeff=1
        if constr_obj.set isa MOI.LessThan
            coeff=-1
        end
        vals = coeff.*(vals .- rhs_const)
    end
    length(vals) == 1 && return vals[1]
    return vals
end

function evaluate(constraint::JuMP.ConstraintRef, data::Union{Dict, DataFrame},
                           model::JuMP.Model)
   return evaluate(constraint, data)
end

"""
    evaluate(constraint::Expr, data::Union{Dict, DataFrame}, model::JuMP.Model)

    Evaluates constraint violation on data in an expression.
    Returns Array of values.

    Note that the keys of the Dict have to be uniform.
"""
function evaluate(constraint::Function, data::Union{Dict, DataFrame}, vars::Array)
    args = []
    clean_data = data_to_DataFrame(data)
    arrays = [[j=1:length()] for i=1:size(clean_data,1)]
    for var_element in vars
    n_vals, n_vars = size(clean_data)
    for (outer, vars) in outer_var_mapping
        if vars isa JuMP.VariableRef
            new_data[outer] = clean_data[!, string(vars)]
        elseif vars isa Array{JuMP.VariableRef}
            new_data[outer] = Inf*ones(size(clean_data,1), length(vars))
            for i=1:length(vars)
                if haskey(clean_data, string(vars[i]))
                    new_data[outer][:,i] = clean_data[!, string(vars[i])]
                end
            end
        else
            throw(ArgumentError(string(constraint, "had an error in evaluation.")))
        end
    end
    new_data = chop_dict(new_data) # makes sure to "scalarize" scalar variable inputs.
    if n_vals == 1
        return substitute_expr(constraint, new_data[1])
    else
        return [substitute_expr(constraint, new_data[i]) for i=1:length(new_data)]
    end
end


text:
# Learn how to copy files and directories

if Sys.iswindows()
    println("Windows")
    src = raw"C:\Users\Yi-Kun.Yang\tmp\org.allnix.julia"
    tgt = raw"C:\TEMP\org.allnix.julia"
    cmd = `cmd /c robocopy /move /s $(src) $(tgt) `

    cerr = IOBuffer()
    cout = IOBuffer()
    proc = run(pipeline(cmd, stdout=cout, stderr=cerr), wait=false)

    wait(proc)
    if proc.exitcode > 1
        @error "exitcode = $(proc.exitcode)"
        err = String(take!(cerr))
        @error err
    end
else
    println("Unsupported OS: $(Sys.MACHINE)")
end


text:
using CircularArrays
using Base.Test

# @test 1 == 1

# # getindex test
# x = reshape([i for i in 1:15],3,5)
# carray = CircularArray(x)
# advance_head!(carray,2)
# @test getindex(carray,1:2,2:5) == [10 13 1 4; 11 14 2 5];

# # shorthand test
# @test carray[1:2,2:5] == [10 13 1 4; 11 14 2 5];

# # setindex! test
# c2 = CircularArray(reshape([i for i in 1:28],4,7))
# advance_head!(c2,3)
# setindex!(c2,[3:6],1:4,3)
# @test c2.data ==
# [1  5   9  13  17  3  25;
#  2  6  10  14  18  4  26;
#  3  7  11  15  19  5  27;
#  4  8  12  16  20  6  28]

# different slicing dimension test
cdim = 1
X = falses(4,3,5)
C = CircularArray(X,cdim)
advance_head!(C,1)
C[1,1:3,1:5] = trues(3,5)
@test reshape(C.data[2,:,:],3,5) == trues(3,5)


text:
abstract type AbstractStopCondition end

function is_terminated(sc::AbstractStopCondition, args...; kwrags...)
    error("Defined this method for type: $(typeof(sc))")
end

struct DistanceStopCondition <: AbstractStopCondition
    eps::Real
    p::Real
    function DistanceStopCondition(eps=1e-1, p=2)
        @assert p > 1
        new(p, eps)
    end
end

function is_terminated(sc::DistanceStopCondition, w, w_new)
    @unpack eps, p = sc
    norm(w_new - w, p) < eps
end


text:
function SetWireInValue(rhd::FPGA, ep, val, mask = 0xffffffff)
    er=ccall(Libdl.dlsym(rhd.lib,:okFrontPanel_SetWireInValue),Cint,(Ptr{Void},Int,Culong,Culong),rhd.board,ep, val, mask)
end

function UpdateWireIns(rhd::FPGA)
    ccall(Libdl.dlsym(rhd.lib,:okFrontPanel_UpdateWireIns),Void,(Ptr{Void},),rhd.board)
    nothing
end

function UpdateWireOuts(rhd::FPGA)
    ccall(Libdl.dlsym(rhd.lib,:okFrontPanel_UpdateWireOuts),Void,(Ptr{Void},),rhd.board)
    nothing
end

function ActivateTriggerIn(rhd::FPGA,epAddr::UInt8,bit::Int)
    er=ccall(Libdl.dlsym(rhd.lib,:okFrontPanel_ActivateTriggerIn),Cint,(Ptr{Void},Int32,Int32),rhd.board,epAddr,bit)
end

function GetWireOutValue(rhd::FPGA,epAddr::UInt8)
    value = ccall(Libdl.dlsym(rhd.lib,:okFrontPanel_GetWireOutValue),Culong,(Ptr{Void},Int32),rhd.board,epAddr)
end

function ReadFromPipeOut(rhd::FPGA,epAddr::UInt8, length, data)
    ccall(Libdl.dlsym(rhd.lib,:okFrontPanel_ReadFromPipeOut),Clong,(Ptr{Void},Int32,Clong,Ptr{UInt8}),rhd.board,epAddr,length,data)
end

function ReadFromBlockPipeOut(fpga::FPGA,epAddr::UInt8,length,data)
ccall(Libdl.dlsym(fpga.lib,:okFrontPanel_ReadFromBlockPipeOut),Clong,(Ptr{Void},Int32,Int32,Clong,Ptr{UInt8}),fpga.board,epAddr,USB3_BLOCK_SIZE,length,data)
end

function ReadUsbBuffer(fpga::FPGA)
    ReadFromPipeOut(fpga,PipeOutData, convert(Clong, fpga.numBytesPerBlock), fpga.usbBuffer)
    nothing
end


text:
using COVID
using Test

@testset "COVID.jl" begin
    params = ParameterSetting()
    switch_pts = [(0.3, 5), (0.9, 5)]
    res1 = run_model(params, .4, 20)
    res2 = run_model(params, switch_pts)
    res3 = run_model(params, 0.03, 6, 20)
    @test typeof(res1) == Array{Int64,2}
    @test typeof(res2) == Array{Int64,2}
    @test typeof(res3) == Array{Float64,2}
end


text:
# -*- coding: utf-8 -*-
using ModelingToolkit
using MacroTools
using SemanticModels
using SemanticModels.PetriModels
import MacroTools: postwalk
using Test
using Petri

# +
@variables S, I, R

# generated by AMIDOL, basic SIR Model
sir = model(PetriModel, Petri.Model([S, I, R],
                 [(S+I, 2I), (I,R)]))

# +
@variables I′, R′

# rule1 = SI <- SI -> SII
si = model(PetriModel, Petri.Model([S, I, R],
                 [(S+I, 2I)]))

sii = model(PetriModel, Petri.Model([S, I, I′, R],
                  [(S+I,  2I ),
                   (S+I′, 2I′)]
                 ))

rule1 = PetriModels.PetriSpan(si, si, sii)

# +
# rule2 = I <- I -> IR
i = model(PetriModel, Petri.Model([I′], Tuple{Operation,Operation}[]))

ir = model(PetriModel, Petri.Model([I′, R′], [(I′, R′)]))

rule2 = PetriModels.PetriSpan(i, i, ir)

# +
# rule3 = IRI′R′ <- II′ -> II′R
irir = model(PetriModel, Petri.Model([I, I′, R, R′], [(I, R), (I′, R′)]))

ii = model(PetriModel, Petri.Model([I, I′], Tuple{Operation,Operation}[]))

iir = model(PetriModel, Petri.Model([I, I′, R], [(I, R), (I′, R)]))

rule3 = PetriModels.PetriSpan(irir, ii, iir)

# +
# Apply rules to each model in succession

siir = PetriModels.solve(PetriModels.DPOProblem(rule1, sir))

siirr = PetriModels.solve(PetriModels.DPOProblem(rule2, siir))

c′ = PetriModels.dropdown(irir, ii, siirr)
siir′ = PetriModels.solve(PetriModels.DPOProblem(rule3, c′))

println(sir.model.Δ)
println()
println(siir′.model.Δ)


text:
struct EnumerableUnique{T,TKEY,S,Q<:Function} <: Enumerable
    source::S
    f::Q
end

function unique(source::Enumerable, f::Function, f_expr::Expr)
    T = eltype(source)
    S = typeof(source)
    TKEY = Base._return_type(f, Tuple{T,})
    return EnumerableUnique{T,TKEY,S,typeof(f)}(source, f)
end

Base.eltype(::Type{EnumerableUnique{T,TKEY,S,Q}}) where {T,TKEY,S,Q} = T

function Base.iterate(iter::EnumerableUnique{T,TKEY,S,Q}) where {T,TKEY,S,Q}
    ret = iterate(iter.source)

    ret===nothing && return nothing

    observed_keys = Set{TKEY}()

    first_element = ret[1]
    source_state = ret[2]

    key_first_element = iter.f(first_element)

    push!(observed_keys, key_first_element)

    return first_element, (observed_keys=observed_keys, source_state=source_state)
end

function Base.iterate(iter::EnumerableUnique{T,TKEY,S,Q}, state) where {T,TKEY,S,Q}
    ret = iterate(iter.source, state.source_state)

    ret===nothing && return nothing

    while true
        current_element = ret[1]
        key_current_element=iter.f(current_element)
        if key_current_element in state.observed_keys
            ret = iterate(iter.source, ret[2])
            ret===nothing && return nothing
        else
            push!(state.observed_keys, key_current_element)

            return current_element, (observed_keys=state.observed_keys, source_state=ret[2])
        end
    end
end


text:
struct ThreadLocalCache{T}
    cache::Vector{Vector{T}}
    size::Int
end

function ThreadLocalCache{T}(; size::Integer = 4) where {T}
    cache = [empty!(Vector{T}(undef, size)) for _ in 1:Threads.nthreads()]
    return ThreadLocalCache(cache, size)
end

function _maybepop!(cache::ThreadLocalCache)
    buffer = cache.cache[Threads.threadid()]
    if isempty(buffer)
        return nothing
    else
        return pop!(buffer)
    end
end

function trypush!(cache::ThreadLocalCache{T}, x::T) where {T}
    buffer = cache.cache[Threads.threadid()]
    if length(buffer) ≥ cache.size
        return false
    else
        push!(buffer, x)
        return true
    end
end


text:
function compute_dists²(x::AbstractMatrix, y::AbstractMatrix)
    y = y'
    return sum(x.^2; dims=2) .+ sum(y.^2; dims=1) .- 2 .* (x * y)
end

@testset "setconv.jl" begin
    for density in [true, false]
        scale = 0.1f0
        n_context = 5
        n_target = 10
        dimensionality = 2
        batch_size = 3
        num_channels = 2

        # Generate context and target set.
        xc = randn(Float32, n_context, dimensionality, batch_size)
        yc = randn(Float32, n_context, num_channels, batch_size)
        xt = randn(Float32, n_target, dimensionality, batch_size)

        # Compute with layer.
        layer = set_conv(num_channels, scale; density=density)
        _, out = NeuralProcesses.code(layer, xc, yc, xt)

        # Brute-force the calculation.
        batches = []
        for i in 1:batch_size
            channels = []

            # Compute weights.
            dists² = compute_dists²(xt[:, :, i], xc[:, :, i]) ./ scale.^2
            weights = NeuralProcesses.rbf(dists²)

            # Prepend density channel.
            density && push!(channels, weights * ones(Float32, n_context))

            # Compute other channels.
            for j in 1:num_channels
                channel = weights * yc[:, j, i]
                # Normalise by density channel.
                density && (channel ./= channels[1] .+ 1f-8)
                push!(channels, channel)
            end

            push!(batches, hcat(channels...))
        end
        ref = cat(batches..., dims=3)

        # Check that the brute-force calculation lines up with the layer.
        @test out ≈ ref
    end
end


text:
using DoctorDocstrings
using Test
using TerminalMenus
using Printf
using TypedTables
using REPL

    module testme
        struct Foo end

        f() = nothing

        """
        hello

        Examples
        =======

        ```jldoctest
        julia> 3+3
        6
        ```
        """
        g() = nothing

        h() = 1

        """
        yolo swag
        """
        i() = 0

        """

        j(type, str; base)

        parse is a number ....

        !!! note

        """
        j() = 0
        export f, g, h, i, j
    end

@testset "DoctorDocstrings.jl" begin
    # Write your tests here.

    @test hasdocs(testme.g)
    @test hasexamples(testme.g)
    @test hasdocs(testme.i)
    @test !hasexamples(testme.i)

   # list = list_no_docs(testme)
   # @test length(list) == 24
   # @test list[2,3] == "-"
   # @test list[3,3] == "Yes"
   # @test all(list[4:7, 3] .== "-")

@test length(list_no_docs(testme)) == 7
@test length(diagnosedocs(testme)) == 3

end


text:
# Copyright 2020 Google LLC
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

const TRIVIAL_LP_MODEL = FirstOrderLp.QuadraticProgrammingProblem(
  [0.0, 1.0], # variable_lower_bound
  [1.0, 2.0], # variable_upper_bound
  zeros(2, 2), # objective_matrix
  [2.0, -1.0], # objective_vector
  0.0, # objective_constant
  [-1.0 -1.0], # constraint_matrix
  [-3.0], # right_hand_side
  0, # num_equalities
)

const TRIVIAL_QP_MODEL = FirstOrderLp.QuadraticProgrammingProblem(
  [0.0, 1.0], # variable_lower_bound
  [1.0, 2.0], # variable_upper_bound
  [2.0 2.0; 2.0 4.0], # objective_matrix
  [2.0, -1.0], # objective_vector
  0.0, # objective_constant
  [-1.0 -1.0], # constraint_matrix
  [-3.0], # right_hand_side
  0, # num_equalities
)

@testset "Read .mps file (LP)" begin
  qp = FirstOrderLp.qps_reader_to_standard_form(
    joinpath(dirname(@__FILE__), "trivial_lp_model.mps"),
  )
  test_fields_equal(qp, TRIVIAL_LP_MODEL)
end

@testset "Read .mps file (QP)" begin
  qp = FirstOrderLp.qps_reader_to_standard_form(
    joinpath(dirname(@__FILE__), "trivial_qp_model.mps"),
  )
  test_fields_equal(qp, TRIVIAL_QP_MODEL)
end

@testset "Read .mps.gz file (QP)" begin
  mps_contents =
    read(open(joinpath(dirname(@__FILE__), "trivial_qp_model.mps")), String)
  gzfile = tempname() * ".mps.gz"
  try
    fd = GZip.open(gzfile, "w")
    write(fd, mps_contents)
    close(fd)
    qp = FirstOrderLp.qps_reader_to_standard_form(gzfile)
    test_fields_equal(qp, TRIVIAL_QP_MODEL)
  finally
    rm(gzfile, force = true)
  end
end

@testset "Two-sided rows to slacks" begin
  qp = FirstOrderLp.TwoSidedQpProblem(
    [-Inf, -Inf],  # variable_lower_bound
    [Inf, Inf],  # variable_upper_bound
    [-3.0, -2.0],  # constraint_lower_bound
    [1.0, Inf],  # constraint_upper_bound
    [1.0 1.0; 1.0 1.0],  # constraint_matrix
    2.0,  # objective_offset
    [0.0, 1.0],  # objective_vector
    diagm([1.0, 3.0]),  # objective_matrix
  )
  FirstOrderLp.two_sided_rows_to_slacks(qp)
  test_fields_equal(
    qp,
    FirstOrderLp.TwoSidedQpProblem(
      [-Inf, -Inf, -3.0],  # variable_lower_bound
      [Inf, Inf, 1.0],  # variable_upper_bound
      [0.0, -2.0],  # constraint_lower_bound
      [0.0, Inf],  # constraint_upper_bound
      [
        1.0 1.0 -1.0
        1.0 1.0 0.0
      ],  # constraint_matrix
      2.0,  # objective_offset
      [0.0, 1.0, 0.0],  # objective_vector
      diagm([1.0, 3.0, 0.0]),  # objective_matrix
    ),
  )
end


text:
cd(@__DIR__)
using AdvancedPS, Random, Distributions

Random.seed!(50514)

priors = (
    μ=(Normal(0, 10),),
    σ=(Truncated(Cauchy(0, 1), 0.0, Inf),)
)

bounds = ((-Inf,Inf),(0.0,Inf))

data = rand(Normal(0.0, 1.0), 50)

function loglike(μ, σ, data)
    return sum(logpdf.(Normal(μ, σ), data))
end

loglike(θ) = loglike(θ..., data)

model = DEModel(priors=priors, model=loglike)

de = DE(bounds=bounds, burnin=1000, priors=priors, progress=true)
n_iter = 2000
chains = psample(model, de, n_iter)
println(chains)


text:
include(joinpath(pwd(), "src/constraints/contact.jl"))
# include(joinpath(pwd(), "src/constraints/contact_no_slip.jl"))

# Simulate contact model one step
function step_contact(model::Model{<: Integration, FixedTime}, x1, u1, w1, h;
        tol_c = 1.0e-5, tol_opt = 1.0e-5, tol_s = 1.0e-4, nlp = :ipopt)
    # Horizon
    T = 2

    # Objective
    obj = PenaltyObjective(1.0e5, model.m)

    # Constraints
    con_dynamics = dynamics_constraints(model, T; w = [w1])
    con_contact = contact_constraints(model, T)
    con = multiple_constraints([con_dynamics, con_contact])

    # Bounds
    _uu = Inf * ones(model.m)
    _uu[model.idx_u] .= u1
    _ul = zeros(model.m)
    _ul[model.idx_u] .= u1
    ul, uu = control_bounds(model, T, _ul, _uu)

    xl, xu = state_bounds(model, T, x1 = x1)

    # Problem
    prob = trajectory_optimization_problem(model,
                   obj,
                   T,
                   h = h,
                   xl = xl,
                   xu = xu,
                   ul = ul,
                   uu = uu,
                   con = con,
                   dynamics = false)

    # Trajectory initialization
    x0 = [x1 for t = 1:T]
    u0 = [[u1; 1.0e-5 * rand(model.m - model.nu)] for t = 1:T-1] # random controls

    # Pack trajectories into vector
    z0 = pack(x0, u0, prob)

    @time z, info = solve(prob, copy(z0),
        tol = tol_opt, c_tol = tol_c, mapl = 0, nlp = nlp)

    @assert check_slack(z, prob) < tol_s
    x, u = unpack(z, prob)

    return x[end]
end


text:
# Submodule for filtering design

###################################################
# Function that calls filters from Wavelets.jl
"""
...
filt_gen(WT.filter::Wavelets.WT) generates de desired filter with the next possible
choices:

* Haar: WT.Haar

* Coiflet: coif_n_ (n in 2:2:8)

* Daubechies: db_n_ (n in 1:Inf)

* Symlet: sym_n_ (n in 4:10)

* Battle: batt_n_ (n in 2:2:6)

* Beylkin: beyl

* Vaidyanathan: vaid

* CDF: cdf_n_ (n in (9,7))

#### Examples
#### filt_gen(WT.haar)
...
"""
function filt_gen(filt1)
	if filt1=="scaling_shearlet"
		filter = [0.0104933261758410,-0.0263483047033631,-0.0517766952966370,0.276348304703363,0.582566738241592,0.276348304703363,-0.0517766952966369,-0.0263483047033631,0.0104933261758408]
	elseif filt1=="directional_shearlet"
		filter = [-0.0000000e+00   0.0000000e+00  -3.0861315e-07   0.0000000e+00  -3.7033578e-07   0.0000000e+00  -4.8143652e-07   0.0000000e+00  -8.3942778e-07   0.0000000e+00  -4.8143652e-07   0.0000000e+00  -3.7033578e-07   0.0000000e+00  -3.0861315e-07   0.0000000e+00  -0.0000000e+00
  -0.0000000e+00   6.1722631e-07  -0.0000000e+00  -6.0488178e-06  -0.0000000e+00  -7.6782952e-06  -0.0000000e+00  -1.2171703e-05  -0.0000000e+00  -1.2171703e-05  -0.0000000e+00  -7.6782952e-06  -0.0000000e+00  -6.0488178e-06  -0.0000000e+00   6.1722631e-07  -0.0000000e+00
  -3.0861315e-07   0.0000000e+00   1.2838307e-05   0.0000000e+00  -5.9401860e-05   0.0000000e+00  -8.7893026e-05   0.0000000e+00  -1.0969346e-04   0.0000000e+00  -8.7893026e-05   0.0000000e+00  -5.9401860e-05   0.0000000e+00   1.2838307e-05   0.0000000e+00  -3.0861315e-07
  -0.0000000e+00  -6.0488178e-06  -0.0000000e+00   1.3512318e-04  -0.0000000e+00  -4.3544081e-04  -0.0000000e+00  -5.7848918e-04  -0.0000000e+00  -5.7848918e-04  -0.0000000e+00  -4.3544081e-04  -0.0000000e+00   1.3512318e-04  -0.0000000e+00  -6.0488178e-06  -0.0000000e+00
  -3.7033578e-07   0.0000000e+00  -5.9401860e-05   0.0000000e+00   1.0726899e-03   6.1035156e-04  -1.9615452e-03   3.6621094e-04  -2.2362850e-03   3.6621094e-04  -1.9615452e-03   6.1035156e-04   1.0726899e-03   0.0000000e+00  -5.9401860e-05   0.0000000e+00  -3.7033578e-07
  -0.0000000e+00  -7.6782952e-06  -0.0000000e+00  -4.3544081e-04  -6.1035156e-04   5.3247620e-03   6.3476562e-03  -5.3114299e-03   4.1503906e-03  -5.3114299e-03   6.3476562e-03   5.3247620e-03  -6.1035156e-04  -4.3544081e-04  -0.0000000e+00  -7.6782952e-06  -0.0000000e+00
  -4.8143652e-07   0.0000000e+00  -8.7893026e-05   0.0000000e+00  -1.9615452e-03  -6.3476562e-03   1.6444291e-02   3.3691406e-02  -8.7910061e-03   3.3691406e-02   1.6444291e-02  -6.3476562e-03  -1.9615452e-03   0.0000000e+00  -8.7893026e-05   0.0000000e+00  -4.8143652e-07
  -0.0000000e+00  -1.2171703e-05  -0.0000000e+00  -5.7848918e-04  -3.6621094e-04  -5.3114299e-03  -3.3691406e-02   3.3130363e-02   1.7749023e-01   3.3130363e-02  -3.3691406e-02  -5.3114299e-03  -3.6621094e-04  -5.7848918e-04  -0.0000000e+00  -1.2171703e-05  -0.0000000e+00
  -8.3942778e-07   0.0000000e+00  -1.0969346e-04   0.0000000e+00  -2.2362850e-03  -4.1503906e-03  -8.7910061e-03  -1.7749023e-01   5.9486541e-01  -1.7749023e-01  -8.7910061e-03  -4.1503906e-03  -2.2362850e-03   0.0000000e+00  -1.0969346e-04   0.0000000e+00  -8.3942778e-07
  -0.0000000e+00  -1.2171703e-05  -0.0000000e+00  -5.7848918e-04  -3.6621094e-04  -5.3114299e-03  -3.3691406e-02   3.3130363e-02   1.7749023e-01   3.3130363e-02  -3.3691406e-02  -5.3114299e-03  -3.6621094e-04  -5.7848918e-04  -0.0000000e+00  -1.2171703e-05  -0.0000000e+00
  -4.8143652e-07   0.0000000e+00  -8.7893026e-05   0.0000000e+00  -1.9615452e-03  -6.3476562e-03   1.6444291e-02   3.3691406e-02  -8.7910061e-03   3.3691406e-02   1.6444291e-02  -6.3476562e-03  -1.9615452e-03   0.0000000e+00  -8.7893026e-05   0.0000000e+00  -4.8143652e-07
  -0.0000000e+00  -7.6782952e-06  -0.0000000e+00  -4.3544081e-04  -6.1035156e-04   5.3247620e-03   6.3476562e-03  -5.3114299e-03   4.1503906e-03  -5.3114299e-03   6.3476562e-03   5.3247620e-03  -6.1035156e-04  -4.3544081e-04  -0.0000000e+00  -7.6782952e-06  -0.0000000e+00
  -3.7033578e-07   0.0000000e+00  -5.9401860e-05   0.0000000e+00   1.0726899e-03   6.1035156e-04  -1.9615452e-03   3.6621094e-04  -2.2362850e-03   3.6621094e-04  -1.9615452e-03   6.1035156e-04   1.0726899e-03   0.0000000e+00  -5.9401860e-05   0.0000000e+00  -3.7033578e-07
  -0.0000000e+00  -6.0488178e-06  -0.0000000e+00   1.3512318e-04  -0.0000000e+00  -4.3544081e-04  -0.0000000e+00  -5.7848918e-04  -0.0000000e+00  -5.7848918e-04  -0.0000000e+00  -4.3544081e-04  -0.0000000e+00   1.3512318e-04  -0.0000000e+00  -6.0488178e-06  -0.0000000e+00
  -3.0861315e-07   0.0000000e+00   1.2838307e-05   0.0000000e+00  -5.9401860e-05   0.0000000e+00  -8.7893026e-05   0.0000000e+00  -1.0969346e-04   0.0000000e+00  -8.7893026e-05   0.0000000e+00  -5.9401860e-05   0.0000000e+00   1.2838307e-05   0.0000000e+00  -3.0861315e-07
  -0.0000000e+00   6.1722631e-07  -0.0000000e+00  -6.0488178e-06  -0.0000000e+00  -7.6782952e-06  -0.0000000e+00  -1.2171703e-05  -0.0000000e+00  -1.2171703e-05  -0.0000000e+00  -7.6782952e-06  -0.0000000e+00  -6.0488178e-06  -0.0000000e+00   6.1722631e-07  -0.0000000e+00
  -0.0000000e+00   0.0000000e+00  -3.0861315e-07   0.0000000e+00  -3.7033578e-07   0.0000000e+00  -4.8143652e-07   0.0000000e+00  -8.3942778e-07   0.0000000e+00  -4.8143652e-07   0.0000000e+00  -3.7033578e-07   0.0000000e+00  -3.0861315e-07   0.0000000e+00  -0.0000000e+00]
	elseif filt1=="directional_shearlet2"
		filter = [0.0000000e+00  -0.0000000e+00   0.0000000e+00  -0.0000000e+00   1.6717973e-03  -0.0000000e+00   0.0000000e+00  -0.0000000e+00   0.0000000e+00
   0.0000000e+00  -0.0000000e+00   0.0000000e+00  -6.6871894e-03  -2.1080148e-03  -6.6871894e-03   0.0000000e+00  -0.0000000e+00   0.0000000e+00
   0.0000000e+00  -0.0000000e+00   1.0030784e-02   6.3240444e-03  -1.9555817e-02   6.3240444e-03   1.0030784e-02  -0.0000000e+00   0.0000000e+00
   0.0000000e+00  -6.6871894e-03  -6.3240444e-03   5.2486012e-02   1.3975610e-01   5.2486012e-02  -6.3240444e-03  -6.6871894e-03   0.0000000e+00
   1.6717973e-03   2.1080148e-03  -1.9555817e-02  -1.3975610e-01   6.8785947e-01  -1.3975610e-01  -1.9555817e-02   2.1080148e-03   1.6717973e-03
   0.0000000e+00  -6.6871894e-03  -6.3240444e-03   5.2486012e-02   1.3975610e-01   5.2486012e-02  -6.3240444e-03  -6.6871894e-03   0.0000000e+00
   0.0000000e+00  -0.0000000e+00   1.0030784e-02   6.3240444e-03  -1.9555817e-02   6.3240444e-03   1.0030784e-02  -0.0000000e+00   0.0000000e+00
   0.0000000e+00  -0.0000000e+00   0.0000000e+00  -6.6871894e-03  -2.1080148e-03  -6.6871894e-03   0.0000000e+00  -0.0000000e+00   0.0000000e+00
   0.0000000e+00  -0.0000000e+00   0.0000000e+00  -0.0000000e+00   1.6717973e-03  -0.0000000e+00   0.0000000e+00  -0.0000000e+00   0.0000000e+00]
	elseif filt1=="directional_shearlet3"
		filter = [-6.0515365e-02   0.0000000e+00   1.2103073e-01   0.0000000e+00  -6.0515365e-02   0.0000000e+00
   4.6875000e-02   7.8125000e-02  -4.6875000e-01   4.6875000e-01  -7.8125000e-02  -4.6875000e-02
  -0.0000000e+00  -6.0515365e-02  -0.0000000e+00   1.2103073e-01  -0.0000000e+00  -6.0515365e-02]
  elseif filt1=="Coiflet1"
		filter = [3.8580778e-02,-1.2696913e-01,-7.7161555e-02,6.0749164e-01,7.4568756e-01,2.2658427e-01]
	else
		#Call the function and add a zero
		filter = [0;wavelet(filt1).qmf]
	end
    filter
end


text:
"""
    @enter ex

Step into `ex`.
"""
macro enter(ex)
  Main.Atom.enter(__module__, ex)
end

"""
    @run ex

Run `ex` with the Interpreter and break when a breakpoint is encountered.
"""
macro run(ex)
  Main.Atom.enter(__module__, ex; initial_continue = true)
end

"""
    @trace ex

Analyse `ex` for type instabilities using [Traceur.jl](https://github.com/JunoLab/Traceur.jl).

!!! warning
    You need to first run `using Traceur` before using this macro.
"""
macro trace(ex, args...)
  return if isdefined(Main.Atom,  :trace)
    Main.Atom.trace(ex, args...)
  else
    notification("`Juno.@trace`"; kind = :Warning, options = (
      dismissable = true,
      description = """
      You haven't loaded [Traceur package](https://github.com/JunoLab/Traceur.jl) into this running session.
      Run `using Traceur` first and use this macro again.
      """
    ))
    :($(esc(ex)))
  end
end

function connect(args...; kws...)
  activate()
  eval(:(Main.Atom.connect($args...; $kws...)))
  return
end


text:
using DynamicalSystems
using DifferentialEquations
using PyPlot
using Interpolations
using LinearAlgebra
using Statistics
using DelimitedFiles

function pyragas_rhs(state, history, params, time)
    chaotic_system, feedback_amplitude_vector, feedback_delay = params
    delayed_state = history(params, time - feedback_delay)
    F = feedback_amplitude_vector .* (delayed_state .- state)
    return chaotic_system.f(state, chaotic_system.p, time) .+ F
end

function pyragas_control(chaotic_system, feedback_amplitude_vector, feedback_delay, tmax;
        dt=0.1, system_transit=0.0, dense=false, ode_algorithm=AutoTsit5(Rosenbrock23()),
        solve_kwargs=(abstol=1e-6, reltol=1e-3, maxiters=10^10))
    ode_problem = ODEProblem(chaotic_system, (- feedback_delay - system_transit, 0.0))
    initial_solution = solve(ode_problem, ode_algorithm; dense=true, solve_kwargs...)
    initial_history(p, t) = initial_solution(t)
    initial_state = initial_solution.u[end]
    params = (chaotic_system, feedback_amplitude_vector, feedback_delay)
    problem = DDEProblem(pyragas_rhs, initial_state, initial_history, (0.0, tmax), params; constant_lags=[feedback_delay])
    dde_algorithm = MethodOfSteps(ode_algorithm)
    return solve(problem, dde_algorithm; solve_kwargs, dense=dense, saveat=dense ? [] : dt)
end

const sys = Systems.roessler()
const system_transit = 1000.0
const v = SVector(0.0, 1.0, 0.0)
const default_K = 0.2
const default_tau = 17.5

function plot_control_example()
    K = default_K
    tau = default_tau
    tmax = 500.0
    t_chaotic = trajectory(sys, tmax; Ttr=system_transit)
    feedback_amplitude_vector = K.*v
    t_control = pyragas_control(sys, feedback_amplitude_vector, tau, tmax; system_transit=system_transit)
    dt = t_control.t[2] - t_control.t[1]
    start = ceil(Int, tau/dt)
    transit_end = 1 + ceil(Int, 480.0/dt)
    F = [dot(feedback_amplitude_vector .* (t_control(t - tau) .- t_control(t)), v) for t in t_control.t[start:end]]

    figure()
    subplot(2,1,1)
    plot(t_chaotic[:, 1], t_chaotic[:, 2], lw=0.5, label="chaos")
    plot(t_control[1, transit_end:end], t_control[2, transit_end:end], lw=2.0, label="control")
    legend(framealpha=1.0)
    xlabel("x")
    ylabel("y")
    subplot(2,1,2)
    plot(t_control.t[start:end], F)
    xlabel("t")
    ylabel("F")
    suptitle("Pyragas Control of the Rossler Oscillator")
    tight_layout()
    savefig("pyragas_control_of_roessler_attractor.pdf")
end

const disp_disp_tau_filename = "pyragas_control_displacement_dispersions_tau.dat"
const disp_disp_K_filename = "pyragas_control_displacement_dispersions_K.dat"

function compute_dispersions_tau()
    tmax = 5000.0
    nsamples = 1000
    control_transit = 1000.0
    taus = range(4.0, 20.0, length=500)
    @time displacement_dispersions = map(taus) do tau
        sol = pyragas_control(sys, default_K.*v, tau, tmax; system_transit=system_transit, dense=true)
        mean((sol(t - tau)[2] - sol(t)[2])^2 for t = tau + control_transit:tmax/nsamples:tmax)
    end
    writedlm(disp_disp_tau_filename, [taus displacement_dispersions])
end

const disp_tau = 5.9

function compute_dispersions_K()
    tmax = 5000.0
    nsamples = 1000
    control_transit = 1000.0
    Ks = range(0.0, 1.1, length=500)
    tau = disp_tau
    @time displacement_dispersions = map(Ks) do K
        sol = pyragas_control(sys, K.*v, tau, tmax; system_transit=system_transit, dense=true)
        mean((sol(t - tau)[2] - sol(t)[2])^2 for t = tau + control_transit:tmax/nsamples:tmax)
    end
    writedlm(disp_disp_K_filename, [Ks displacement_dispersions])
end

function plot_dispersions()
    figure()

    data = readdlm(disp_disp_tau_filename)
    taus, displacement_dispersions = data[:, 1], data[:, 2]
    subplot(2,1,1)
    plot(taus, displacement_dispersions, ".")
    xlabel(raw"$\tau$")
    ylabel(raw"$\langle D(t)^2 \rangle$")
    yscale("log")
    ylim(1e-11, 1e2)
    title("K = $default_K")

    data = readdlm(disp_disp_K_filename)
    Ks, displacement_dispersions = data[:, 1], data[:, 2]
    subplot(2,1,2)
    plot(Ks, displacement_dispersions, ".")
    xlabel("K")
    ylabel(raw"$\langle D(t)^2 \rangle$")
    yscale("log")
    # ylim(1e-11, 1e2)
    title("\$\\tau\$ = $disp_tau")

    tight_layout()
    savefig("pyragas_control_displacement_dispersion_plot.pdf")
end

const max_lyap_filename = "pyragas_max_lyap.dat"

function pyragas_map!(next_history_vector, _initial_history_vector, params, _t)
    chaotic_system, feedback_amplitude_vector, feedback_delay, N = params

    initial_history_vector = [SVector{N}(_initial_history_vector[i:i+N-1]) for i = 1:N:length(_initial_history_vector) - N]

    dt = feedback_delay / (length(initial_history_vector) - 1)
    initial_state = initial_history_vector[end]
    interp = LinearInterpolation(range(-feedback_delay, 0.0, length=length(initial_history_vector)), initial_history_vector)
    initial_history(p, t) = interp(t)
    problem = DDEProblem(pyragas_rhs, initial_state, initial_history, (0.0, feedback_delay + dt), params; constant_lags=[feedback_delay])

    dde_algorithm = MethodOfSteps(AutoTsit5(Rosenbrock23()))
    sol = solve(problem, dde_algorithm; saveat=dt)

    for (i, u) in enumerate(@view(sol.u[2:end]))
        copy!(@view(next_history_vector[(i - 1)*N + 1 : i*N]), u)
    end
    return nothing
end

function max_lyapunov(K, tau; N=80)
    dt = tau / (N - 1)
    sol = pyragas_control(sys, K.*v, tau, tau; system_transit=1000.0, dt=dt)
    ds = DiscreteDynamicalSystem(pyragas_map!, vcat(sol.u...), [sys, K.*v, tau, length(sys.u0)])
    return lyapunov(ds, 100; Ttr=round(Int, 500.0/(tau + dt)), show_progress=true)
end

function compute_max_lyapunov_exponent()
    Ks = range(0.0, 1.1, length=10) # 10 .^ range(-2, 2, length=10)
    @time lambda_maxs = map(K -> max_lyapunov(K, disp_tau), Ks)
    writedlm(max_lyap_filename, [Ks lambda_maxs])
end

function plot_max_lyapunov_exponent()
    data = readdlm(max_lyap_filename)
    Ks, lambda_maxs = data[:, 1], data[:, 2]
    figure()
    plot(Ks, lambda_maxs, ".")
    axhline(0.0, color="black")
    xlabel("K")
    ylabel(raw"$\lambda_\mathrm{max}$")
    #xscale("log")
    savefig("pyragas_max_lyapunov_exponent.pdf")
end

function main()
    plot_control_example()
    compute_dispersions_tau()
    compute_dispersions_K()
    plot_dispersions()
    compute_max_lyapunov_exponent()
    plot_max_lyapunov_exponent()
end

function old_max_lyapunov(K, tau; debug=false)
    tmax = 5000.0
    ks = 1:40
    steps = 1000.0:1.0:tmax
    sol = pyragas_control(sys, K.*v, tau, tmax; system_transit=1000.0, dense=true)
    data = Dataset(sol.(steps))
    Es = lyapunov_from_data(data, ks)
    starts, slops = linear_regions(ks, Es; tol=1.0)
    if debug
        for i = 1:length(starts) - 1
            interval = ks[i] : ks[i + 1]
            plot(interval, Es[interval])
        end
    end
    return slops[1]
end

max_stars_repo_name: cosmo-jana/numerics-physics-stuff

max_issues_repo_name: cosmo-jana/numerics-physics-stuff

max_forks_repo_name: cosmo-jana/numerics-physics-stuff


text:
lines = readlines("input.txt")
octomap = parse.(Int, reduce(hcat, collect.(lines)))
idx = CartesianIndices(octomap)
minval = typemin(Int)

o = CartesianIndex(1,1)
V = CartesianIndices((-1:1, -1:1))
padded = fill(minval, size(octomap) .+ 2)
padded[idx .+ o] = octomap[idx]

Σ = 0
for step = 1:100
    padded .+= 1
    flashmap = fill(false, size(padded))
    while (yxs = findall(padded .>= 10)) |> length > 0
        flashmap[yxs] .= true
        padded[yxs] .= minval # Prevent any further flashing for this step
        [padded[yx .+ V] .+= 1 for yx = yxs]
    end
    padded[findall(flashmap)] .= 0
    global Σ += sum(flashmap)
end

println(Σ)


text:
export BackLinks

struct BackLinks
  postlinks::String

  function BackLinks()
    default_postlinks = """

    """
    return new(default_postlinks)
  end

  BackLinks(default_pl::String) = new(default_pl)
end

function process_html(B::BackLinks,
                      html::AbstractString,
                      fileinfo::FileInfo;
                      args...)
  folders = split(fileinfo["folderstring"],"/")[2:end]
  filename = fileinfo["filename"]

  bhtml = ""
  if length(folders) > 0 || filename!="index.md" #<-- don't show for main page
    bhtml *= "
	\n"
    bhtml *= "main/"
    tfold = "/"
    for fold in folders[1:end-1]
      tfold *= fold * "/"
      bhtml *= "$fold/"
    end
    if filename!="index.md"
      if length(folders) > 0
          fold = folders[end]
          tfold *= fold * "/"
          bhtml *= "$fold/"
      end
      bhtml *= fileinfo["basename"]
    else
      if length(folders) > 0
        bhtml*= "$(folders[end])/"
      end
    end
    bhtml *= "
"
  end
  return bhtml*B.postlinks*html
end


text:
using Images
using ImageView
using Boltzmann

function nonzero_indexes{T <: Number}(mat::Matrix{T})
    idxs = (Int, Int)[]
    for i=1:size(mat, 1)
        for j=1:size(mat, 2)
            if mat[i, j] != 0
                push!(idxs, (i, j))
            end
        end
    end
    return idxs
end

function collect_nonzeros{T <: Number}(mat::Matrix{T}, idxs::Array{(Int, Int)})
    arr = zeros(T, length(idxs))
    k = 1
    for (i, j) in idxs
        arr[k] = mat[i, j]
        k += 1
    end
    return arr
end

function map_nonzeros{T <: Number}(matsize::(Int, Int), arr::Array{T, 1},
                      idxs::Array{(Int, Int)})
    mat = zeros(T, matsize)
    k = 1
    for (i, j) in idxs
        mat[i, j] = arr[k]
        k += 1
    end
    return mat
end

IMSIZE = (96, 96)

using Images

function imresize!(resized, original)
    scale1 = (size(original,1)-1)/(size(resized,1)-0.999f0)
    scale2 = (size(original,2)-1)/(size(resized,2)-0.999f0)
    for jr = 0:size(resized,2)-1
        jo = scale2*jr
        ijo = itrunc(jo)
        fjo = jo - oftype(jo, ijo)
        @inbounds for ir = 0:size(resized,1)-1
            io = scale1*ir
            iio = itrunc(io)
            fio = io - oftype(io, iio)
            tmp = (1-fio)*((1-fjo)*original[iio+1,ijo+1] +
                           fjo*original[iio+1,ijo+2]) +
            fio*((1-fjo)*original[iio+2,ijo+1] +
                   fjo*original[iio+2,ijo+2])
            resized[ir+1,jr+1] = convertsafely(eltype(resized), tmp)
        end
    end
    resized
end
imresize(original, new_size) =
    imresize!(similar(original,
                      new_size), original)

convertsafely{T<:FloatingPoint}(::Type{T}, val) = convert(T, val)
convertsafely{T<:Integer}(::Type{T}, val::Integer) = convert(T, val)
convertsafely{T<:Integer}(::Type{T}, val::FloatingPoint) =
    itrunc(T, val+oftype(val, 0.5))

# Create dataset consiting of nonzero pixels of face images
# Returns:
#   dat : Matrix(# of nonzeros x # of filenames) - data matrix
#   nonzero_idxs : Vector((i, j)) - nonzero pixel coordinates
function facedata(datadir="../../data/CK/faces_aligned", imsize=IMSIZE)
    filenames = readdir(datadir)
    refmat = convert(Array, imread(datadir * "/" * filenames[1]))
    refmat = imresize(refmat, imsize)
    # refmat = refmat ./ 256  -- imresize already does it
    nonzero_idxs = nonzero_indexes(refmat)
    dat = zeros(length(nonzero_idxs), length(filenames))
    for (i, fname) in enumerate(filenames)
        if i % 500 == 0 println(i, " ", fname) end
        mat = convert(Array, imread(datadir * "/" * fname))
        mat = imresize(mat, imsize)
        # mat = mat ./ 256
        dat[:, i] = collect_nonzeros(mat, nonzero_idxs)
    end
    return dat, nonzero_idxs
end

function wplot(W, nzs, padding=10)
    h, w = IMSIZE
    n = size(W, 1)
    rows = int(floor(sqrt(n)))
    cols = int(ceil(n / rows))
    halfpad = div(padding, 2)
    dat = zeros(rows * (h + padding), cols * (w + padding))
    for i=1:n
        wt = W[i, :]
        wt = reshape(wt, length(wt))
        wim = map_nonzeros(IMSIZE, wt, nzs)
        wim = wim ./ (maximum(wim) - minimum(wim))
        r = div(i - 1, cols) + 1
        c = rem(i - 1, cols) + 1
        dat[(r-1)*(h+padding)+halfpad+1 : r*(h+padding)-halfpad,
            (c-1)*(w+padding)+halfpad+1 : c*(w+padding)-halfpad] = wim
    end
    view(dat)
    return dat
end

# REPL

using HDF5, JLD, PyPlot

flatten(a) = reshape(a, prod(size(a)))
ihist(a) = plt.hist(flatten(a), bins=100)

function run(n_hid=1024)
    dat, nzs = facedata()
    n_feat, n_samples = size(dat)
    model = GRBM(n_feat, n_hid, sigma=0.001)
    for i=1:10
        println("meta-iteration #", i)
        fit(model, dat, n_iter=10, n_gibbs=3, lr=0.01)
        println("Sleeping for 20 seconds to cool down")
        sleep(20)
    end
    min_images = min(n_hid, 64)
    for i=1:min(3, div(n_hid, min_images))
        wplot(model.weights[(i-1)*min_images+1:i*min_images, :], nzs)
    end
    ## rest = rem(n_hid, min_images)
    ## if rest > 0
    ##     wplot(model.weights[end-rest:end, :], nzs)
    ## end
    # wplot(model.weights[div(n_hid)*min_images+1:end, :], nzs)
    return model, nzs
end

# some results
# n_hid=192, n_gibbs=3, n_meta=10 ~ 7 faces of good quality
# n_hid=192, n_gibbs=5, n_meta=10 ~ 4 faces of good quality
# n_hid=256, n_gibbs=1, n_meta=10 ~ no faces at all
# n_hid=256, n_gibbs=3, n_meta=10 ~ 10 faces of almost good quality
# n_hid=256, n_gibbs=5, n_meta=10 ~ 11 faces of almost good quality
# n_hid=256, n_gibbs=10, n_meta=10 ~ 11 faces of almost good quality
# n_hid=256, n_gibbs=5, n_meta=20 ~ 12 faces of almost good quality
# n_hid=128, n_gibbs=3, n_meta=10 ~ all faces of bad quality
# n_hid=192, n_gibbs=3, n_meta=10 ~ 3 faces of a modate quality
# n_hid=256, n_gibbs=3, n_meta=10 ~ 10 faces of moderate quality
# n_hid=320, n_gibbs=3, n_meta=10 ~ 16 faces of almost good qualiy
# n_hid=384, n_gibbs=3, n_meta=10 ~ 17 faces of almost good qualiy
# n_hid=512, n_gibbs=3, n_meta=10 ~ 22 faces of different quality
# n_hid=768, n_gibbs=3, n_meta=10 ~ 25 faces of different quality
# n_hid=1024, n_gibbs=3, n_meta=10 ~ 21 face of different quality
# n_hid=768, n_gibbs=3, n_meta=10, sigma=0.001 ~ almost all faces!
# n_hid=512, n_gibbs=3, n_meta=10, sigma=0.001 ~ all good, but many similar
# n_hid=1024, n_gibbs=3, n_meta=10, sigma=0.001 ~ all good
# n_hid=128, n_gibbs=3, n_meta=10, lr=0.01 ~ 1 good, others similar and bad
# n_hid=128, n_gibbs=10, n_meta=10, lr=0.01 ~ 1 good, others similar and bad
# brbm: n_hid=768, n_gibbs=3, n_meta=10, lr=0.01 ~ 20 good
# grbm: n_hid=1024, n_gibbs=3, n_meta=10 ~ half almost good
# grbm: n_hid=768, n_gibbs=3, n_meta=10 ~ half almost good
# grbm: n_hid=1024, n_gibbs=10, n_meta=10 ~ half almost good
# grbm: n_hid=1024, n_gibbs=10, n_meta=3, lr=0.01 ~ half almost good
#! grbm: n_hid=1024, n_gibbs=3, n_meta=10, sigma=0.001, lr=0.01 ~ all really good!
# brbm: n_hid=1024, n_gibbs=3, n_meta=10, sigma=0.001, lr=0.01 ~ almost the same

# grbm (w/ momentum): n_hid=1024, n_gibbs=2, n_meta=5, sigma=0.001, lr=0.01 ~ 2243

function load_model(filename="session.jld")
    d = load(filename)
    return d["m"], d["nzs"]
end

## if isinteractive()
##     im = imread("../../data/CK/faces_aligned/S005_001_00000010.png")
##     mat = convert(Array, im)
## end

max_stars_repo_name: dfdx/masque

max_issues_repo_name: dfdx/masque

max_forks_repo_name: dfdx/masque


text:
@info("Testing basic random number generators...")
for brng_type in brng_types
    brng = BasicRandomNumberGenerator(brng_type, rand(UInt))
    brng = BasicRandomNumberGenerator(brng_type, [rand(Cuint) for i in 1:10])
end


text:
using LinearAlgebra
using Fermi.DIIS
# Functions specific to CTF
#
############ PREPROCESSING ################
"""
    Fermi.CoupledCluster.RCCSD{T}(Alg::CTF)

Compute a RCCSD wave function using the Compiled time factorization algorithm (CTF)
"""
function BCCD{T}(guess::BCCD{Tb},Alg::CTF) where { T <: AbstractFloat,
                                                    Tb <: AbstractFloat }
    refwfn = Fermi.HartreeFock.RHF()

    drop_occ = Fermi.CurrentOptions["drop_occ"]
    drop_vir = Fermi.CurrentOptions["drop_vir"]

    @output "Transforming Integrals..."
    tint = @elapsed begin
        refwfn.ints["OOOO"]
        refwfn.ints["OOOV"]
        refwfn.ints["OOVV"]
        refwfn.ints["OVOV"]
        refwfn.ints["OVVV"]
        refwfn.ints["VVVV"]
    end
    @output " done in {} s" tint
    BCCD{T}(refwfn, guess, Alg)
end

function BCCD{T}(refwfn::RHF, guess::BCCD{Tb}, Alg::CTF) where { T <: AbstractFloat,
                                                                  Tb <: AbstractFloat }
    ints = refwfn.ints
    d = [i - a for i = diag(ints["FOO"]), a = diag(ints["FVV"])]
    D = [i + j - a - b for i = diag(ints["FOO"]), j = diag(ints["FOO"]), a = diag(ints["FVV"]), b = diag(ints["FVV"])]
    newT1 = ints["FOV"]./d
    newT2 = ints["OOVV"]./D
    o_small,v_small = size(guess.T1)
    o = 1:o_small
    v = 1:v_small
    newT1[o,v] .= Fermi.data(guess.T1)
    newT2[o,o,v,v] .= Fermi.data(guess.T2)
    BCCD{T}(refwfn, ints, newT1, newT2, Alg)
end
############ END PREPROCESSING ################
#
#
#
############# KERNEL FUNCTIONS ###################
function print_bcc_alg(alg::CTF)
    @output "Computing BCC with CTF algorithm.\n"
end

function compute_integrals(ints,alg::CTF)
    ints["OOOO"]
    ints["OOOV"]
    ints["OOVV"]
    ints["OVOV"]
    ints["OVVV"]
    ints["VVVV"]
end

function delete_integrals(ints,alg::CTF)
    delete!(ints.cache,"OOOO")
    delete!(ints.cache,"OOOV")
    delete!(ints.cache,"OOVV")
    delete!(ints.cache,"OVOV")
    delete!(ints.cache,"OVVV")
    delete!(ints.cache,"VVVV")
    delete!(ints.cache,"F")
    delete!(ints.cache,"FOO")
    delete!(ints.cache,"FVV")
    delete!(ints.cache,"FOV")
end

function compute_oovv(ints,alg::CTF)
    ints["OOVV"]
end

function update_T1(T2::Array{T,4}, newT1::Array{T,2}, foo, fov, fvv, ints::IntegralHelper, alg::CTF) where T <: AbstractFloat
    fill!(newT1, 0.0)
    Voooo, Vooov, Voovv, Vovov, Vovvv, Vvvvv = ints["OOOO"], ints["OOOV"], ints["OOVV"], ints["OVOV"], ints["OVVV"], ints["VVVV"]
    @tensoropt (i=>x, j=>x, k=>x, l=>x, a=>10x, b=>10x, c=>10x, d=>10x) begin
        newT1[i,a] += fov[i,a]
        newT1[i,a] += 2.0*fov[k,c]*T2[i,k,a,c]
        newT1[i,a] -= fov[k,c]*T2[k,i,a,c]
        newT1[i,a] -= T2[k,i,c,d]*Vovvv[k,a,d,c]
        newT1[i,a] += 2.0*T2[i,k,c,d]*Vovvv[k,a,d,c]
        newT1[i,a] += -2.0*T2[k,l,a,c]*Vooov[k,l,i,c]
        newT1[i,a] += T2[l,k,a,c]*Vooov[k,l,i,c]
    end
end

function update_T2(T2::Array{T,4},newT2::Array{T,4},foo,fov,fvv,ints::IntegralHelper, alg::CTF) where T <: AbstractFloat
    fill!(newT2, 0.0)
    Voooo, Vooov, Voovv, Vovov, Vovvv, Vvvvv = ints["OOOO"], ints["OOOV"], ints["OOVV"], ints["OVOV"], ints["OVVV"], ints["VVVV"]

    @tensoropt (i=>x, j=>x, k=>x, l=>x, a=>10x, b=>10x, c=>10x, d=>10x) begin
        newT2[i,j,a,b] += Voovv[i,j,a,b]
        newT2[i,j,a,b] += T2[i,j,c,d]*Vvvvv[c,d,a,b]
        newT2[i,j,a,b] += T2[k,l,a,b]*Voooo[i,j,k,l]
        newT2[i,j,a,b] += T2[k,l,a,c]*T2[i,j,d,b]*Voovv[k,l,c,d]
        newT2[i,j,a,b] += -2.0*T2[i,k,a,c]*T2[l,j,b,d]*Voovv[k,l,c,d]
        newT2[i,j,a,b] += -2.0*T2[l,k,a,c]*T2[i,j,d,b]*Voovv[k,l,c,d]
        newT2[i,j,a,b] += T2[k,i,a,c]*T2[l,j,d,b]*Voovv[l,k,c,d]
        newT2[i,j,a,b] += T2[i,k,a,c]*T2[l,j,b,d]*Voovv[l,k,c,d]
        newT2[i,j,a,b] += -2.0*T2[i,k,a,c]*T2[j,l,b,d]*Voovv[l,k,c,d]
        newT2[i,j,a,b] += T2[k,i,a,c]*T2[l,j,b,d]*Voovv[k,l,c,d]
        newT2[i,j,a,b] += -2.0*T2[k,i,a,c]*T2[j,l,b,d]*Voovv[k,l,c,d]
        newT2[i,j,a,b] += T2[i,j,a,c]*T2[l,k,b,d]*Voovv[k,l,c,d]
        newT2[i,j,a,b] += -2.0*T2[i,j,a,c]*T2[k,l,b,d]*Voovv[k,l,c,d]
        newT2[i,j,a,b] += T2[k,j,a,c]*T2[i,l,d,b]*Voovv[l,k,c,d]
        newT2[i,j,a,b] += 4.0*T2[i,k,a,c]*T2[j,l,b,d]*Voovv[k,l,c,d]
        newT2[i,j,a,b] += T2[i,j,d,c]*T2[l,k,a,b]*Voovv[k,l,c,d]
        P_OoVv[i,j,a,b] := -1.0*foo[i,k]*T2[k,j,a,b]
        P_OoVv[i,j,a,b] += fvv[c,a]*T2[i,j,c,b]
        P_OoVv[i,j,a,b] += -1.0*T2[k,i,a,c]*Voovv[k,j,c,b]
        P_OoVv[i,j,a,b] += 2.0*T2[i,k,a,c]*Voovv[k,j,c,b]
        P_OoVv[i,j,a,b] += -1.0*T2[i,k,a,c]*Vovov[j,c,k,b]
        P_OoVv[i,j,a,b] += -1.0*T2[k,j,a,c]*Vovov[i,c,k,b]
        P_OoVv[i,j,a,b] += T2[j,k,c,d]*T2[i,l,a,b]*Voovv[k,l,c,d]
        P_OoVv[i,j,a,b] += -2.0*T2[i,k,d,c]*T2[l,j,a,b]*Voovv[k,l,c,d]

        newT2[i,j,a,b] += P_OoVv[i,j,a,b] + P_OoVv[j,i,b,a]
    end
end
#
#
############# END KERNEL FUNCTIONS ###################


text:
using OceanColorData
using Test

@testset "OceanColorData.jl" begin

wvbd_out=[412, 443, 490, 510, 555, 670]
wvbd_in=[400,425,450,475,500,525,550,575,600,625,650,675,700]
rirr_in=1e-3*[23.7641,26.5037,27.9743,30.4914,28.1356,
    21.9385,18.6545,13.5100,5.6338,3.9272,2.9621,2.1865,1.8015]
rrs_ref=1e-3*[4.4099, 4.8533, 5.1247, 4.5137, 3.0864, 0.4064]
chla_ref=0.6101219875535429
fcm_ref=[5.683734365402806e-26, 4.958850399758996e-64, 2.0058235898296587e-43,
  2.399985306394786e-23, 1.0914931954631329e-15, 3.610712128260008e-8,
  0.06660122938880934, 0.4566726319832639, 0.001153891939362191, 5.901424234631198e-11,
  1.8433658038301877e-10, 0.4133547888920851, 0.027225883646784382, 0.005595266829573927]

rrs_out=RemotelySensedReflectance(rirr_in,wvbd_in,wvbd_out)
chla_out=RrsToChla(rrs_out; Eqn="OC4")
(M,Sinv)=Jackson2017()
fcm_out=FuzzyClassification(M,Sinv,rrs_out)

    @test isapprox(rrs_out,rrs_ref,atol=1e-5)
    @test isapprox(chla_out,chla_ref,atol=1e-5)
    @test isapprox(fcm_out,fcm_ref,atol=1e-5)
end


text:
using Documenter, AbstractStorage

makedocs(
    sitename="AbstractStorage",
    modules=[AbstractStorage],
    pages = [
        "index.md",
        "links.md",
        "reference.md"
        ]
)

deploydocs(
    repo = "git@github.com:ChevronETC/AbstractStorage.jl.git"
)


text:
using GeoPhyInv
using Statistics
using Test

medium = Medium(:elastic_homo2D, 5)

medium_new = similar(medium)
update!(medium_new, [:vp, :rho, :vs], rectangle = [[-500, -500], [500, 500]], perc = 5.0) # perturb velocity

ageom = AGeom(medium.mgrid, :surf, SSrcs(3), Recs(100)); # surface seismic

tgrid = range(0.0, stop = 2.0, length = 2500) # generate a temporal grid
wav = ricker(15.0, tgrid, tpeak = 0.25); # Choose a source wavelet
srcwav = SrcWav(tgrid, ageom, [:vz]) # initialize
update!(srcwav, [:vz], wav) # distribute to all supersources

pa = SeisForwExpt(
    FdtdElastic(),
    medium = medium,
    ageom = ageom,
    srcwav = srcwav,
    tgrid = tgrid,
    rfields = [:vz],
    verbose = true,
);

@time update!(pa);
d1 = copy(pa[:data][:vz])

update!(pa, medium_new)

@time update!(pa);
d2 = copy(pa[:data][:vz])

@test d1 ≠ d2


text:
Qx, x = PolynomialRing(FlintQQ, "x")
K2, a2 = NumberField(x^3 - 2, "a1")
K3, (a3,) = NumberField([x^3 - 2], "a2")
@testset "Fractional ideals for $K1" for (K1, a1) in [(K2, a2), (K3, a3)]
  O1 = Order(K1, Hecke.FakeFmpqMat(FlintZZ[1 0 0; 0 2 0; 0 0 4], one(FlintZZ)))

  i = ideal(O1, O1(2*a1))

  @testset "Construction" begin
    I = @inferred fractional_ideal(O1, i)
    @test basis_matrix(I) == Hecke.FakeFmpqMat(FlintZZ[16 0 0; 0 1 0; 0 0 1], fmpz(1))
    @test basis_mat_inv(I) == Hecke.FakeFmpqMat(FlintZZ[1 0 0; 0 16 0; 0 0 16], fmpz(16))

    J = @inferred fractional_ideal(O1, i, 2)
    @test basis_matrix(J) == Hecke.FakeFmpqMat(FlintZZ[16 0 0; 0 1 0; 0 0 1], fmpz(2))
    @test basis_matrix(J) == Hecke.FakeFmpqMat(FlintZZ[16 0 0; 0 1 0; 0 0 1], fmpz(2))
    @test basis_matrix(J) == Hecke.FakeFmpqMat(FlintZZ[16 0 0; 0 1 0; 0 0 1], fmpz(2))
    @test basis_mat_inv(J) == Hecke.FakeFmpqMat(FlintZZ[1 0 0; 0 16 0; 0 0 16], fmpz(8))
    @test basis_mat_inv(J) == Hecke.FakeFmpqMat(FlintZZ[1 0 0; 0 16 0; 0 0 16], fmpz(8))

    K = @inferred fractional_ideal(O1, FlintZZ[16 0 0; 0 1 0; 0 0 1], fmpz(2))
    @test K.basis_matrix == Hecke.FakeFmpqMat(FlintZZ[16 0 0; 0 1 0; 0 0 1], fmpz(2))
    @test basis_matrix(K) == Hecke.FakeFmpqMat(FlintZZ[16 0 0; 0 1 0; 0 0 1], fmpz(2))
    @test basis_mat_inv(K) == Hecke.FakeFmpqMat(FlintZZ[1 0 0; 0 16 0; 0 0 16], fmpz(8))
    @test basis_mat_inv(K) == Hecke.FakeFmpqMat(FlintZZ[1 0 0; 0 16 0; 0 0 16], fmpz(8))

    L = @inferred fractional_ideal(O1, Hecke.FakeFmpqMat(FlintZZ[16 0 0; 0 1 0; 0 0 1], fmpz(2)))
    @test L.basis_matrix == Hecke.FakeFmpqMat(FlintZZ[16 0 0; 0 1 0; 0 0 1], fmpz(2))
    @test basis_matrix(L) == Hecke.FakeFmpqMat(FlintZZ[16 0 0; 0 1 0; 0 0 1], fmpz(2))
    @test basis_mat_inv(L) == Hecke.FakeFmpqMat(FlintZZ[1 0 0; 0 16 0; 0 0 16], fmpz(8))
    @test basis_mat_inv(L) == Hecke.FakeFmpqMat(FlintZZ[1 0 0; 0 16 0; 0 0 16], fmpz(8))

    M = @inferred fractional_ideal(O1, 2*a1//2)
    @test M.basis_matrix == Hecke.FakeFmpqMat(FlintZZ[16 0 0; 0 1 0; 0 0 1], fmpz(2))
    @test basis_matrix(M) == Hecke.FakeFmpqMat(FlintZZ[16 0 0; 0 1 0; 0 0 1], fmpz(2))
    @test basis_mat_inv(M) == Hecke.FakeFmpqMat(FlintZZ[1 0 0; 0 16 0; 0 0 16], fmpz(8))
    @test basis_mat_inv(M) == Hecke.FakeFmpqMat(FlintZZ[1 0 0; 0 16 0; 0 0 16], fmpz(8))

    @test J == K
    @test K == L
    @test L == M
    @test M == J
  end

  J = fractional_ideal(O1, i, 2)
  K = fractional_ideal(O1, FlintZZ[16 0 0; 0 1 0; 0 0 1], fmpz(2))
  L = fractional_ideal(O1, Hecke.FakeFmpqMat(FlintZZ[16 0 0; 0 1 0; 0 0 1], fmpz(2)))
  M = fractional_ideal(O1, 2*a1//2)

  @testset "Basis" begin
    b = @inferred basis(J)
    @test b == [ K1(8), a1, 2*a1^2 ]
    b = @inferred basis(K)
    @test b == [ K1(8), a1, 2*a1^2 ]
    b = @inferred basis(L)
    @test b == [ K1(8), a1, 2*a1^2 ]
    b = @inferred basis(M)
    @test b == [ K1(8), a1, 2*a1^2 ]
  end

  @testset "Norm" begin
    b = @inferred norm(J)
    @test b == FlintQQ(16, 2^3)
  end

  @testset "Ring of multipliers" begin
    i = ideal(O1, FlintZZ[2 0 0; 0 1 0; 0 0 1])
    N = @inferred ring_of_multipliers(i)

    @test basis_matrix(N) == Hecke.FakeFmpqMat(FlintZZ[1 0 0; 0 2 0; 0 0 2], fmpz(1))
  end

  @testset "Denominator" begin
    # This was once a bug found by Johannes
    @testset begin
      R, x = PolynomialRing(FlintQQ, "x")
      K, a = NumberField(x, "a")
      O = maximal_order(K)
      I = Hecke.NfOrdFracIdl(ideal(O, O(2)), fmpz(2))
      @test denominator(I) == fmpz(2)
      basis_matrix(I)
      @test denominator(I) == fmpz(2)
    end
  end
end


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

"""
    Scale(; low=0.25, high=0.75)

Applies the scale transform to all columns of the table.
The scale transform of the column `x` is defined by `(x .- xl) ./ (xh - xl)`,
where `xl = quantile(x, low)` and `xh = quantile(x, high)`.

# Examples

```julia
Scale()
Scale(low=0, high=1)
Scale(low=0.3, high=0.7)
```

## Notes

* The `low` and `high` values are restricted to the interval [0, 1].
"""
struct Scale{T<:Real} <: Colwise
  low::T
  high::T

  function Scale(low::T, high::T) where {T<:Real}
    @assert 0 ≤ low ≤ high ≤ 1 "invalid quantiles"
    new{T}(low, high)
  end
end

Scale(low::Real, high::Real) = Scale(promote(low, high)...)

Scale(; low=0.25, high=0.75) = Scale(low, high)

assertions(::Type{<:Scale}) = [assert_continuous]

isrevertible(::Type{<:Scale}) = true

function colcache(transform::Scale, x)
  levels = (transform.low, transform.high)
  xl, xh = quantile(x, levels)
  xl == xh && ((xl, xh) = (zero(xl), one(xh)))
  (xl=xl, xh=xh)
end

colapply(::Scale, x, c)  = @. (x - c.xl) / (c.xh - c.xl)

colrevert(::Scale, y, c) = @. (c.xh - c.xl) * y + c.xl

"""
    MinMax()

The transform that is equivalent to `Scale(low=0, high=1)`.

See also [`Scale`](@ref).
"""
MinMax() = Scale(low=0, high=1)

"""
    Interquartile()

The transform that is equivalent to `Scale(low=0.25, high=0.75)`.

See also [`Scale`](@ref).
"""
Interquartile() = Scale(low=0.25, high=0.75)