text:
function binary_op(op, W, @nospecialize(T))
  ty = LLVM_TYPES[T]
  if isone(W)
    V = T
  else
    ty = "<$W x $ty>"
    V = NTuple{W,VecElement{T}}
  end
  instrs = "%res = $op $ty %0, %1\nret $ty %res"
  call = :($LLVMCALL($instrs, $V, Tuple{$V,$V}, data(v1), data(v2)))
  W > 1 && (call = Expr(:call, :Vec, call))
  Expr(:block, Expr(:meta, :inline), call)
end

# Integer
for (op, f) ∈ [("add", :+), ("sub", :-), ("mul", :*), ("shl", :<<)]
  ff = Symbol('v', op)
  fnsw = Symbol(ff, "_nsw")
  fnuw = Symbol(ff, "_nuw")
  fnw = Symbol(ff, "_nw")
  ff_fast = Symbol(ff, :_fast)
  @eval begin
    # @inline $ff(a,b) = $ff_fast(a,b)
    @inline $ff(a::T, b::T) where {T<:Union{FloatingTypes,IntegerTypesHW,AbstractSIMD}} =
      $ff_fast(a, b)
    # @generated $ff_fast(v1::Vec{W,T}, v2::Vec{W,T}) where {W,T<:IntegerTypesHW} = binary_op($op * (T <: Signed ? " nsw" : " nuw"), W, T)
    @generated $ff_fast(v1::Vec{W,T}, v2::Vec{W,T}) where {W,T<:IntegerTypesHW} =
      binary_op($op, W, T)
    @generated $fnsw(v1::Vec{W,T}, v2::Vec{W,T}) where {W,T<:IntegerTypesHW} =
      binary_op($(op * " nsw"), W, T)
    @generated $fnuw(v1::Vec{W,T}, v2::Vec{W,T}) where {W,T<:IntegerTypesHW} =
      binary_op($(op * " nuw"), W, T)
    @generated $fnw(v1::Vec{W,T}, v2::Vec{W,T}) where {W,T<:IntegerTypesHW} =
      binary_op($(op * " nsw nuw"), W, T)
    # @generated Base.$f(v1::Vec{W,T}, v2::Vec{W,T}) where {W,T<:IntegerTypesHW} = binary_op($op, W, T)
    @inline Base.$f(v1::Vec{W,T}, v2::Vec{W,T}) where {W,T<:IntegerTypesHW} =
      $ff_fast(v1, v2)
    # @generated $ff(v1::Vec{W,T}, v2::Vec{W,T}) where {W,T<:IntegerTypesHW} = binary_op($op, W, T)
    @inline $ff_fast(x, y) = $f(x, y)

    # @generated $ff_fast(v1::T, v2::T) where {T<:IntegerTypesHW} = binary_op($op * (T <: Signed ? " nsw" : " nuw"), 1, T)
    # @generated $ff_fast(v1::T, v2::T) where {T<:IntegerTypesHW} = binary_op($op, 1, T)
    @generated $fnsw(v1::T, v2::T) where {T<:IntegerTypesHW} =
      binary_op($(op * " nsw"), 1, T)
    @generated $fnuw(v1::T, v2::T) where {T<:IntegerTypesHW} =
      binary_op($(op * " nuw"), 1, T)
    @generated $fnw(v1::T, v2::T) where {T<:IntegerTypesHW} =
      binary_op($(op * " nsw nuw"), 1, T)
  end
end
@inline vadd_fast(v1::T, v2::T) where {T<:IntegerTypesHW} = Base.add_int(v1, v2)
@inline vsub_fast(v1::T, v2::T) where {T<:IntegerTypesHW} = Base.sub_int(v1, v2)
@inline vmul_fast(v1::T, v2::T) where {T<:IntegerTypesHW} = Base.mul_int(v1, v2)
@inline vshl_fast(v1::T, v2::T) where {T<:IntegerTypesHW} = Base.shl_int(v1, v2)
for (op, f) ∈ [("div", :÷), ("rem", :%)]
  ff = Symbol('v', op) #_ff = Symbol(:_, ff)
  sbf = Symbol('s', op, :_int)
  ubf = Symbol('u', op, :_int)
  sbf128 = Sys.WORD_SIZE == 32 ? f : sbf
  ubf128 = Sys.WORD_SIZE == 32 ? f : ubf
  @eval begin
    @generated $ff(v1::Vec{W,T}, v2::Vec{W,T}) where {W,T<:Integer} =
      binary_op((T <: Signed ? 's' : 'u') * $op, W, T)
    @inline $ff(a::I, b::I) where {I<:SignedHW} = Base.$sbf(a, b)
    @inline $ff(a::I, b::I) where {I<:UnsignedHW} = Base.$ubf(a, b)
    @inline $ff(a::Int128, b::Int128) = Base.$sbf128(a, b)
    @inline $ff(a::UInt128, b::UInt128) = Base.$ubf128(a, b)
    # @generated $_ff(v1::T, v2::T) where {T<:Integer} = binary_op((T <: Signed ? 's' : 'u') * $op, 1, T)
    # @inline $ff(v1::T, v2::T) where {T<:IntegerTypesHW} = $_ff(v1, v2)
  end
end
# for (op,f) ∈ [("div",:÷),("rem",:%)]
#   @eval begin
#     @generated Base.$f(v1::Vec{W,T}, v2::Vec{W,T}) where {W,T<:Integer} = binary_op((T <: Signed ? 's' : 'u') * $op, W, T)
#     @generated Base.$f(v1::T, v2::T) where {T<:IntegerTypesHW} = binary_op((T <: Signed ? 's' : 'u') * $op, 1, T)
#   end
# end
@inline vcld(x, y) = vadd(vdiv(vsub(x, one(x)), y), one(x))
@inline function vdivrem(x, y)
  d = vdiv(x, y)
  r = vsub(x, vmul(d, y))
  d, r
end
for (op, sub) ∈ [
  ("ashr", :SignedHW),
  ("lshr", :UnsignedHW),
  ("lshr", :IntegerTypesHW),
  ("and", :IntegerTypesHW),
  ("or", :IntegerTypesHW),
  ("xor", :IntegerTypesHW),
]
  ff = sub === :UnsignedHW ? :vashr : Symbol('v', op)
  @eval begin
    @generated $ff(v1::Vec{W,T}, v2::Vec{W,T}) where {W,T<:$sub} = binary_op($op, W, T)
    @generated $ff(v1::T, v2::T) where {T<:$sub} = binary_op($op, 1, T)
  end
end

for (op, f) ∈ [("fadd", :vadd), ("fsub", :vsub), ("fmul", :vmul), ("fdiv", :vfdiv)]#,("frem",:vrem)]
  ff = Symbol(f, :_fast)
  fop_fast = f === :vfdiv ? "fdiv fast" : op * ' ' * fast_flags(true)
  fop_contract = op * ' ' * fast_flags(false)
  @eval begin
    @generated $f(v1::Vec{W,T}, v2::Vec{W,T}) where {W,T<:Union{Float32,Float64}} =
      binary_op($fop_contract, W, T)
    @generated $ff(v1::Vec{W,T}, v2::Vec{W,T}) where {W,T<:Union{Float32,Float64}} =
      binary_op($fop_fast, W, T)
    @inline $f(v1::Vec{W,Float16}, v2::Vec{W,Float16}) where {W} =
      $f(convert(Vec{W,Float32}, v1), convert(Vec{W,Float32}, v2))
    @inline $ff(v1::Vec{W,Float16}, v2::Vec{W,Float16}) where {W} =
      $ff(convert(Vec{W,Float32}, v1), convert(Vec{W,Float32}, v2))
  end
end
@inline vsub(a::T, b::T) where {T<:Union{Float32,Float64}} = Base.sub_float(a, b)
@inline vadd(a::T, b::T) where {T<:Union{Float32,Float64}} = Base.add_float(a, b)
@inline vmul(a::T, b::T) where {T<:Union{Float32,Float64}} = Base.mul_float(a, b)
@inline vsub_fast(a::T, b::T) where {T<:Union{Float32,Float64}} = Base.sub_float_fast(a, b)
@inline vadd_fast(a::T, b::T) where {T<:Union{Float32,Float64}} = Base.add_float_fast(a, b)
@inline vmul_fast(a::T, b::T) where {T<:Union{Float32,Float64}} = Base.mul_float_fast(a, b)

@inline vdiv(v1::AbstractSIMD{W,T}, v2::AbstractSIMD{W,T}) where {W,T<:FloatingTypes} =
  trunc(vfdiv_fast(v1, v2))
@inline vdiv_fast(v1::AbstractSIMD{W,T}, v2::AbstractSIMD{W,T}) where {W,T<:FloatingTypes} =
  trunc(vfdiv_fast(v1, v2))
@inline vdiv_fast(v1::T, v2::T) where {T<:FloatingTypes} =
  trunc(Base.FastMath.div_float_fast(v1, v2))
@inline vdiv_fast(v1::T, v2::T) where {T<:Number} = v1 ÷ v2
@inline vdiv(v1::T, v2::T) where {T<:Number} = v1 ÷ v2
@inline vdiv(v1::T, v2::T) where {T<:FloatingTypes} = vdiv_fast(v1, v2)
@inline vrem(a, b) = vfnmadd(vdiv_fast(a, b), b, a)
@inline vrem_fast(a, b) = vfnmadd(vdiv_fast(a, b), b, a)
# @inline vdiv_fast(v1::AbstractSIMD{W,T}, v2::AbstractSIMD{W,T}) where {W,T<:IntegerTypesHW} = trunc(T, vfloat_fast(v1) / vfloat_fast(v2))
@inline vdiv_fast(
  v1::AbstractSIMD{W,T},
  v2::AbstractSIMD{W,T},
) where {W,T<:IntegerTypesHW} = trunc(T, vfloat(v1) / vfloat(v2))
@inline function vdiv_fast(v1, v2)
  v3, v4 = promote_div(v1, v2)
  vdiv_fast(v3, v4)
end
@inline vdiv_fast(v1::VecUnroll{N,1,T,T}, s::T) where {N,T<:SignedHW} =
  VecUnroll(fmap(Base.sdiv_int, data(v1), s))
@inline vdiv_fast(v1::VecUnroll{N,1,T,T}, s::T) where {N,T<:UnsignedHW} =
  VecUnroll(fmap(Base.udiv_int, data(v1), s))
@inline vdiv_fast(v1::VecUnroll{N,1,T,T}, v2::VecUnroll{N,1,T,T}) where {N,T<:SignedHW} =
  VecUnroll(fmap(Base.sdiv_int, data(v1), data(v2)))
@inline vdiv_fast(v1::VecUnroll{N,1,T,T}, v2::VecUnroll{N,1,T,T}) where {N,T<:UnsignedHW} =
  VecUnroll(fmap(Base.udiv_int, data(v1), data(v2)))

@inline vfdiv(a::AbstractSIMDVector{W}, b::AbstractSIMDVector{W}) where {W} =
  vfdiv(vfloat(a), vfloat(b))
# @inline vfdiv_fast(a::AbstractSIMDVector{W}, b::AbstractSIMDVector{W}) where {W} = vfdiv_fast(vfloat_fast(a), vfloat_fast(b))
@inline vfdiv_fast(a::AbstractSIMDVector{W}, b::AbstractSIMDVector{W}) where {W} =
  vfdiv_fast(vfloat(a), vfloat(b))
@inline vfdiv(a, b) = a / b
@inline vfdiv_fast(a, b) = Base.FastMath.div_fast(a, b)

for (f, op) ∈ ((:vand, :and_int), (:vor, :or_int), (:vxor, :xor_int))
  @eval @inline $f(b1::Bool, b2::Bool) = Base.$op(b1, b2)
end

for f ∈ [:vadd, :vsub, :vmul]
  for s ∈ [Symbol(""), :_fast, :_nsw, :_nuw, :_nw]
    fs = Symbol(f, s)
    @eval begin
      @inline function $fs(
        a::Union{FloatingTypes,IntegerTypesHW,AbstractSIMD},
        b::Union{FloatingTypes,IntegerTypesHW,AbstractSIMD},
      )
        c, d = promote(a, b)
        $fs(c, d)
      end
    end
  end
end
# @inline vsub(a::T, b::T) where {T<:Base.BitInteger} = Base.sub_int(a, b)
for (vf, bf) ∈ [
  (:vadd, :add_int),
  (:vsub, :sub_int),
  (:vmul, :mul_int),
  (:vadd_fast, :add_int),
  (:vsub_fast, :sub_int),
  (:vmul_fast, :mul_int),
  (:vadd_nsw, :add_int),
  (:vsub_nsw, :sub_int),
  (:vmul_nsw, :mul_int),
  (:vadd_nuw, :add_int),
  (:vsub_nuw, :sub_int),
  (:vmul_nuw, :mul_int),
  (:vadd_nw, :add_int),
  (:vsub_nw, :sub_int),
  (:vmul_nw, :mul_int),
]
  @eval begin
    @inline $vf(a::Int128, b::Int128) = Base.$bf(a, b)
    @inline $vf(a::UInt128, b::UInt128) = Base.$bf(a, b)
  end
end
# @inline vrem(a::Float32, b::Float32) = Base.rem_float_fast(a, b)
# @inline vrem(a::Float64, b::Float64) = Base.rem_float_fast(a, b)

@inline function Base.FastMath.add_fast(a::AbstractSIMD, b::AbstractSIMD, c::AbstractSIMD)
  Base.FastMath.add_fast(Base.FastMath.add_fast(a, b), c)
end
@inline function Base.FastMath.add_fast(a::T, b::T, c::T) where {T<:AbstractSIMD}
  Base.FastMath.add_fast(Base.FastMath.add_fast(a, b), c)
end

@inline function Base.FastMath.add_fast(
  a::AbstractSIMD,
  b::AbstractSIMD,
  c::AbstractSIMD,
  d::AbstractSIMD,
)
  x = Base.FastMath.add_fast(a, b)
  y = Base.FastMath.add_fast(c, d)
  Base.FastMath.add_fast(x, y)
end
@inline function Base.FastMath.add_fast(a::T, b::T, c::T, d::T) where {T<:AbstractSIMD}
  x = Base.FastMath.add_fast(a, b)
  y = Base.FastMath.add_fast(c, d)
  Base.FastMath.add_fast(x, y)
end

@inline function Base.FastMath.add_fast(
  a::AbstractSIMD,
  b::AbstractSIMD,
  c::AbstractSIMD,
  d::AbstractSIMD,
  e::AbstractSIMD,
  f::Vararg{Number,K},
) where {K}
  x = Base.FastMath.add_fast(a, b)
  y = Base.FastMath.add_fast(c, d)
  Base.FastMath.add_fast(Base.FastMath.add_fast(x, y), e, f...)
end
@inline function Base.FastMath.add_fast(
  a::T,
  b::T,
  c::T,
  d::T,
  e::T,
  f::Vararg{T,K},
) where {T<:AbstractSIMD,K}
  x = Base.FastMath.add_fast(a, b)
  y = Base.FastMath.add_fast(c, d)
  Base.FastMath.add_fast(Base.FastMath.add_fast(x, y), e, f...)
end


text:
# THIS FILE IS MACHINE-GENERATED
# See `../generate_osi_list` for more information.
const OSI_LICENSES = Set(["0BSD",
                          "AAL",
                          "AFL-1.1",
                          "AFL-1.2",
                          "AFL-2.0",
                          "AFL-2.1",
                          "AFL-3.0",
                          "AGPL-3.0",
                          "AGPL-3.0-only",
                          "AGPL-3.0-or-later",
                          "APL-1.0",
                          "APSL-1.0",
                          "APSL-1.1",
                          "APSL-1.2",
                          "APSL-2.0",
                          "Apache-1.1",
                          "Apache-2.0",
                          "Artistic-1.0",
                          "Artistic-1.0-Perl",
                          "Artistic-1.0-cl8",
                          "Artistic-2.0",
                          "BSD-1-Clause",
                          "BSD-2-Clause",
                          "BSD-2-Clause-Patent",
                          "BSD-3-Clause",
                          "BSD-3-Clause-LBNL",
                          "BSL-1.0",
                          "CAL-1.0",
                          "CAL-1.0-Combined-Work-Exception",
                          "CATOSL-1.1",
                          "CDDL-1.0",
                          "CECILL-2.1",
                          "CNRI-Python",
                          "CPAL-1.0",
                          "CPL-1.0",
                          "CUA-OPL-1.0",
                          "ECL-1.0",
                          "ECL-2.0",
                          "EFL-1.0",
                          "EFL-2.0",
                          "EPL-1.0",
                          "EPL-2.0",
                          "EUDatagrid",
                          "EUPL-1.1",
                          "EUPL-1.2",
                          "Entessa",
                          "Fair",
                          "Frameworx-1.0",
                          "GPL-2.0",
                          "GPL-2.0+",
                          "GPL-2.0-only",
                          "GPL-2.0-or-later",
                          "GPL-3.0",
                          "GPL-3.0+",
                          "GPL-3.0-only",
                          "GPL-3.0-or-later",
                          "GPL-3.0-with-GCC-exception",
                          "HPND",
                          "IPA",
                          "IPL-1.0",
                          "ISC",
                          "Intel",
                          "LGPL-2.0",
                          "LGPL-2.0+",
                          "LGPL-2.0-only",
                          "LGPL-2.0-or-later",
                          "LGPL-2.1",
                          "LGPL-2.1+",
                          "LGPL-2.1-only",
                          "LGPL-2.1-or-later",
                          "LGPL-3.0",
                          "LGPL-3.0+",
                          "LGPL-3.0-only",
                          "LGPL-3.0-or-later",
                          "LPL-1.0",
                          "LPL-1.02",
                          "LPPL-1.3c",
                          "LiLiQ-P-1.1",
                          "LiLiQ-R-1.1",
                          "LiLiQ-Rplus-1.1",
                          "MIT",
                          "MPL-1.0",
                          "MPL-1.1",
                          "MPL-2.0",
                          "MPL-2.0-no-copyleft-exception",
                          "MS-PL",
                          "MS-RL",
                          "MirOS",
                          "Motosoto",
                          "MulanPSL-2.0",
                          "Multics",
                          "NASA-1.3",
                          "NCSA",
                          "NGPL",
                          "NPOSL-3.0",
                          "NTP",
                          "Naumen",
                          "Nokia",
                          "OCLC-2.0",
                          "OFL-1.1",
                          "OFL-1.1-RFN",
                          "OFL-1.1-no-RFN",
                          "OGTSL",
                          "OSET-PL-2.1",
                          "OSL-1.0",
                          "OSL-2.0",
                          "OSL-2.1",
                          "OSL-3.0",
                          "PHP-3.0",
                          "PHP-3.01",
                          "PostgreSQL",
                          "Python-2.0",
                          "QPL-1.0",
                          "RPL-1.1",
                          "RPL-1.5",
                          "RPSL-1.0",
                          "RSCPL",
                          "SISSL",
                          "SPL-1.0",
                          "SimPL-2.0",
                          "Sleepycat",
                          "UCL-1.0",
                          "UPL-1.0",
                          "Unlicense",
                          "VSL-1.0",
                          "W3C",
                          "Watcom-1.0",
                          "Xnet",
                          "ZPL-2.0",
                          "Zlib"])

ext: jl


text:
FsumSq(x)=([sum(x.*x)],nothing) # must always return an array for the value
FsumSq_inplace(value,auxvalue,x)=copy!(value,sum(x.*x))

function FsumSq(x...)
    tmp=0.0
    for i in 1:length(x)
        tmp+=sum(x[i].*x[i])
    end
    return ([tmp],nothing) # must always return an array for the value
end

function FsumSq_inplace(value,auxvalue,x...) # inplace
    tmp=0.0
    for i in 1:length(x)
        tmp+=sum(x[i].*x[i])
    end
    copy!(value,tmp)
end

DsumSq(derivativeIDX,f_c,faux_c,grad_c,grad_n,x_n...)=axpy!(2.0,grad_c.*x_n[derivativeIDX],grad_n)

if PROC=="GPU"
    function FsumSq(x::CudaArray...)
        tmp=CudaArray(zeros(1))
        for i in 1:length(x)
            CUBLAS.gemv!('T',1.0,flatten(Float64,x[i]),vec(x[i]),1.0,tmp)
        end
    return (tmp,nothing) # must always return an array for the value
    end

    function FsumSq_inplace(value,auxvalue,x::CudaArray...) # inplace
        fill!(value,0.0)
        for i in 1:length(x)
            CUBLAS.gemv!('T',1.0,flatten(Float64,x[i]),vec(x[i]),1.0,value);
        end
    end

    function DsumSq(derivativeIDX,f_c,faux_c,grad_c,grad_n,x_n::CudaArray...)
        alphaaxpy!(2.0,grad_c,x_n[derivativeIDX],grad_n)
    end

end

Inplace[FsumSq]=FsumSq_inplace
Derivative[FsumSq]=DsumSq

sumSq(n::ADnode)=ADnode(FsumSq,n)

function sumSq(n::ArrayADnode)
    return ADnode(FsumSq,n)
end

sumSquare=sumSq

export sumSq
export sumSquare

max_stars_repo_name: davidbarber/AutoDiff

max_issues_repo_name: davidbarber/AutoDiff

max_forks_repo_name: davidbarber/AutoDiff


text:
# ------------ Runtime ------------ #

(d::Type{D})(args::Diffed...) where D <: Distributions.Distribution = d(map(unwrap, args)...)
(d::Type{D})(arg1::Diffed, others...) where D <: Distributions.Distribution = d(unwrap(arg1), others...)

# ------------ Abstract ------------ #

@abstract DiffInterpreter (::Distribution)(args...) = propagate(args...)
@abstract DiffInterpreter (::Type{D})(args...) where D <: Distributions.Distribution = propagate(args...)

@abstract DiffInterpreter rand(::D) where D <: Distribution = Change
@abstract DiffInterpreter rand(::Symbol, args...) = propagate(args...)
@abstract DiffInterpreter rand(::Random._GLOBAL_RNG, ::D) where D <: Distribution = Change
@abstract DiffInterpreter rand(::Symbol, call::Function, args...) = propagate(args...)

@abstract DiffInterpreter trace(::D) where D <: Distribution = Change
@abstract DiffInterpreter trace(::Symbol, args...) = consume(args...)
@abstract DiffInterpreter trace(::Random._GLOBAL_RNG, ::D) where D <: Distribution = Change
@abstract DiffInterpreter trace(::Symbol, call::Function, args...) = consume(args...)
@abstract DiffInterpreter trace(args...) = consume(args...)


text:
module Express

using ..GraphObjects: Figure

using PyCall

abstract type PxModule end
abstract type PxClass end
export px

struct Px end
px = Px()

const plotly = PyNULL()
const express = PyNULL()

const px_functions = [
    :scatter,
    :scatter_3d,
    :scatter_polar,
    :scatter_ternary,
    :scatter_mapbox,
    :scatter_geo,
    :scatter_matrix,
    :density_contour,
    :density_heatmap,
    :density_mapbox,
    :line,
    :line_3d,
    :line_polar,
    :line_ternary,
    :line_mapbox,
    :line_geo,
    :parallel_coordinates,
    :parallel_categories,
    :area,
    :bar,
    :timeline,
    :bar_polar,
    :violin,
    :box,
    :strip,
    :histogram,
    :ecdf,
    :choropleth,
    :choropleth_mapbox,
    :pie,
    :sunburst,
    :treemap,
    :icicle,
    :funnel,
    :funnel_area,
    :imshow,
    :set_mapbox_access_token,
    :get_trendline_results,
]

sym2obj = Dict{Symbol,Union{Function,DataType}}()

for func in px_functions
    @eval begin
        function $(func)(args...; kwargs...)
            ret = express.$(func)(args...; kwargs...)
            if pybuiltin("isinstance")(ret, plotly.graph_objs._figure.Figure)
                Figure(ret)
            end
        end
        sym2obj[nameof($func)] = $func
    end
end

const px_modules = [:data, :colors, :trendline_functions]

for m in px_modules
    @eval begin
        struct $(m) <: PxModule
            pyobj::PyObject
            $(m)() = new(getproperty(express, nameof($m)))
        end

        PyObject(t::$(m)) = t.pyobj

        function Base.propertynames(t::$(m))
            propertynames(getfield(t, :pyobj))
        end

        function Base.getproperty(t::$(m), s::Symbol)
            if s in fieldnames($(m))
                return getfield(t, s)
            else
                return getproperty(getfield(t, :pyobj), s)
            end
        end
        sym2obj[nameof($m)] = $(m)
    end
end

const px_classes = [:IdentityMap, :Constant, :Range]
for class in px_classes
    @eval begin
        struct $(class) <: PxClass
            pyobj::PyObject
            function $(class)(args..., ; kwargs...)
                new(express.$(class)(args...; kwargs...))
            end
        end

        PyObject(t::$(class)) = t.pyobj

        function Base.propertynames(t::$(class))
            propertynames(getfield(t, :pyobj))
        end

        function Base.getproperty(t::$(class), s::Symbol)
            if s in fieldnames($(class))
                return getfield(t, s)
            else
                return getproperty(getfield(t, :pyobj), s)
            end
        end
        sym2obj[nameof($class)] = $(class)
    end
end

function Base.getproperty(px::Px, s::Symbol)
    if s in fieldnames(Px)
        return getfield(px, s)
    else
        if s in px_modules
            return sym2obj[s]().pyobj
        else
            return sym2obj[s]
        end
    end
end
Base.propertynames(px::Px) = vcat(px_functions, px_modules, px_classes)

function __init__()
    pyimport_conda("pandas", "pandas")
    copy!(plotly, pyimport_conda("plotly", "plotly", "plotly"))
    copy!(express, pyimport_conda("plotly.express", "plotly", "plotly"))
end
# Write your package code here.

end


text:
# Implementing Forward Mode
# =========================
#
# In the [intro notebook](./intro.ipynb) we covered forward-mode differentiation
# thoroughly, but we don't yet have a real implementation that can work on our
# programs. Implementing AD effectively and efficiently is a field of its own,
# and we'll need to learn a few more tricks to get off the ground.

include("utils.jl");

# Up to now, we have differentiated things by creating a new Wengert list which
# contains parts of the original expression.

y = Wengert(:(5sin(log(x))))
derive(y, :x)

# We're now going to explicitly split our lists into two pieces: the original
# expression, and a new one which only calculates derivatives (but might refer
# back to values from the first). For example:

y = Wengert(:(5sin(log(x))))
#-
dy = derive(y, :x, out = Wengert(variable = :dy))
#-
Expr(dy)

# If we want to distinguish them, we can call `y` the *primal code* and `dy`
# the *tangent code*. Nothing fundamental has changed here, but it's useful
# to organise things this way.
#
# Almost all of the subtlety in differentiating programs comes from a
# mathematically trivial question: in what order do we evaluate the statements
# of the Wengert list? We have discussed the
# [forward/reverse](./backandforth.ipynb) distinction, but even once that choice
# is made, we have plenty of flexibility, and those choices can affect efficiency.
#
# For example, imagine if we straightforwardly evaluate `y` followed by `dy`. If
# we only cared about the final output of `y`, this would be no problem at all,
# but in general `dy` also needs to re-use variables like `y1` (or possibly any
# $y_i$). If our primal Wengert list has, say, a billion instructions, we end up
# having to store a billion intermediate $y_i$ before we run our tangent code.
#
# Alternatively, one can imagine running each instruction of the tangent code as
# early as possible; as soon as we run `y1 = log(x)`, for example, we know we
# can run `dy2 = cos(y1)` also. Then our final, combined program would look
# something like this:

# ```julia
# y0 = x
# dy = 1
# y1 = log(y0)
# dy = dy/y0
# y2 = cos(y1)
# dy = dy*sin(y1)
#   ...
# ```

# Now we can throw out `y1` soon after creating it, and we no longer have to
# store those billion intermediate results.
#
# The ability to do this is a very general property of forward differentiation;
# once we run $a = f(b)$, we can then run $\frac{da}{dx} = \frac{da}{db}
# \frac{db}{dx}$ using only `a` and `b`. It's really just a case of replacing
# basic instructions like `cos` with versions that calculate both the primal and
# tangent at once.

# Dual Numbers
# ------------
#
# Finally, the trick that we've been building up to: making our programming
# language do this all for us! Almost all common languages – with the notable
# exception of C – provide good support for *operator overloading*, which allows
# us to do exactly this replacement.
#
# To start with, we'll make a container that holds both a $y$ and a
# $\frac{dy}{dx}$, called a *dual number*.

struct Dual{T<:Real} <: Real
  x::T
  ϵ::T
end

Dual(1, 2)
#-
Dual(1.0,2.0)

# Let's print it nicely.

Base.show(io::IO, d::Dual) = print(io, d.x, " + ", d.ϵ, "ϵ")

Dual(1, 2)

# And add some of our rules for differentiation. The rules have the same basic
# pattern-matching structure as the ones we originally applied to our Wengert
# list, just with different notation.

import Base: +, -, *, /
a::Dual + b::Dual = Dual(a.x + b.x, a.ϵ + b.ϵ)
a::Dual - b::Dual = Dual(a.x - b.x, a.ϵ - b.ϵ)
a::Dual * b::Dual = Dual(a.x * b.x, b.x * a.ϵ + a.x * b.ϵ)
a::Dual / b::Dual = Dual(a.x * b.x, b.x * a.ϵ - a.x * b.ϵ)

Base.sin(d::Dual) = Dual(sin(d.x), d.ϵ * cos(d.x))
Base.cos(d::Dual) = Dual(cos(d.x), - d.ϵ * sin(d.x))

Dual(2, 2) * Dual(3, 4)

# Finally, we'll hook into Julia's number promotion system; this isn't essential
# to understand, but just makes it easier to work with Duals since we can now
# freely mix them with other number types.

Base.convert(::Type{Dual{T}}, x::Dual) where T = Dual(convert(T, x.x), convert(T, x.ϵ))
Base.convert(::Type{Dual{T}}, x::Real) where T = Dual(convert(T, x), zero(T))
Base.promote_rule(::Type{Dual{T}}, ::Type{R}) where {T,R} = Dual{promote_type(T,R)}

Dual(1, 2) * 3

# We already have enough to start taking derivatives of some simple functions.
# If we pass a dual number into a function, the $\epsilon$ component represents
# the derivative.

f(x) = x / (1 + x*x)

f(Dual(5., 1.))

# We can make a utility which allows us to differentiate any function.

D(f, x) = f(Dual(x, one(x))).ϵ

D(f, 5.)

# Dual numbers seem pretty scarcely related to all the Wengert list stuff we
# were talking about earlier. But we need take a closer look at how this is
# working. To start with, look at Julia's internal representation of `f`.

@code_typed f(1.0)

# This is just a Wengert list! Though the naming is a little different – `mul_float`
# rather than the more general `*` and so on – it's still essentially the same
# data structure we were working with earlier. Moreover, you'll recognise
# the code for the derivative, too!

@code_typed D(f, 1.0)

# This code is again the same as the Wengert list derivative we worked out at
# the very beginning of this handbook. The order of operations is just a little
# different, and there's the odd missing or new instruction due to the different
# set of optimisations that Julia applies. Still, we have not escaped our fundamental
# symbolic differentiation algorithm, just tricked the compiler into doing most
# of the work for us.

derive(Wengert(:(sin(cos(x)))), :x)
#-
@code_typed D(x -> sin(cos(x)), 0.5)

# What of data structures, control flow, function calls? Although these things
# are all present in Julia's internal "Wengert list", they end up being the same
# in the tangent program as in the primal; so an operator overloading approach
# need not deal with them explicitly to do the right thing. This won't be true
# when we come to talk more about reverse mode, which demands a more complex
# approach.

# Perturbation Confusion
# ----------------------
#
# Actually, that's not quite true. Operator-overloading-based forward mode
# *almost always* does the right thing, but it is not flawless. This more
# advanced section will talk about nested differentiation and the nasty bug that
# can come with it.
#
# We can differentiate any function we want, as long as we have the right
# primitive definitions for it. For example, the derivative of $\sin(x)$ is
# $\cos(x)$.

D(sin, 0.5), cos(0.5)

# We can also differentiate the differentiation operator itself. We'll find that
# the second derivative of $\sin(x)$ is $-\sin(x)$.

D(x -> D(sin, x), 0.5), -sin(0.5)

# This worked because we ended up nesting dual numbers. If we create a dual number
# whose $\epsilon$ component is another dual number, then we end up tracking the
# derivative of the derivative.

# The issue comes about when we close over a variable that *is itself* being
# differentiated.

D(x -> x*D(y -> x+y, 1), 1) # == 1

# The derivative $\frac{d}{dy} (x + y) = 1$, so this is equivalent to
# $\frac{d}{dx}x$, which should also be $1$. So where did this go wrong? The
# problem is that when we closed over $x$, we didn't just get a numeric value
# but a dual number with $\epsilon = 1$. When we then calculated $x + y$, both
# epsilons were added as if $\frac{dx}{dy} = 1$ (effectively $x = y$). If we had
# written this down, the answer would be correct.

D(x -> x*D(y -> y+y, 1), 1)

# I leave this second example as an excercise to the reader. Needless to say,
# this has caught out many an AD implementor.

D(x -> x*D(y -> x*y, 1), 4) # == 8

# More on Dual Numbers
# --------------------
#
# The above discussion presented dual numbers as essentially being a trick for
# applying the chain rule. I wanted to take the opportunity to present an
# alternative viewpoint, which might be appealing if, like me, you have any
# training in physics.
#
# Complex arithmetic involves a new number, $i$, which behaves like no other:
# specifically, because $i^2 = -1$. We'll introduce a number called $\epsilon$,
# which is a bit like $i$ except that $\epsilon^2 = 0$; this is effectively a
# way of saying the $\epsilon$ is a very small number. The relevance of this
# comes from the original definition of differentiation, which also requires
# $\epsilon$ to be very small.
#
# $$
# \frac{d}{dx} f(x) = \lim_{\epsilon \to 0} \frac{f(x+\epsilon)-f(x)}{\epsilon}
# $$
#
# We can see how our definition of $\epsilon$ works out by applying it to
# $f(x+\epsilon)$; let's say that $f(x) = sin(x^2)$.
#
# \begin{align}
# f(x + \epsilon) &= \sin((x + \epsilon)^2) \\
#                 &= \sin(x^2 + 2x\epsilon + \epsilon^2) \\
#                 &= \sin(x^2 + 2x\epsilon) \\
#                 &= \sin(x^2)\cos(2x\epsilon) + \cos(x^2)\sin(2x\epsilon) \\
#                 &= \sin(x^2) + 2x\cos(x^2)\epsilon \\
# \end{align}
#
# A few things have happened here. Firstly, we directly expand $(x+\epsilon)^2$
# and remove the $\epsilon^2$ term. We expand $sin(a+b)$ and then apply a *small
# angle approximation*: for small $\theta$, $\sin(\theta) \approx \theta$ and
# $\cos(\theta) \approx 1$. (This sounds pretty hand-wavy, but does follow from
# our original definition of $\epsilon$ if we look at the Taylor expansion of
# both functions). Finally we can plug this into our derivative rule.
#
# \begin{align}
# \frac{d}{dx} f(x) &= \frac{f(x+\epsilon)-f(x)}{\epsilon} \\
#                   &= 2x\cos(x^2)
# \end{align}
#
# This is, in my opinion, a rather nice way to derive functions by hand.
#
# This also leads to another nice trick, and a third way to look at forward-mode
# AD; if we replace $x + \epsilon$ with $x + \epsilon i$ then we still have
# $(\epsilon i)^2 = 0$. If $\epsilon$ is a small real number (say
# $1\times10^{-10}$), this is still true within floating point error, so our
# derivative still works out.

ϵ = 1e-10im
x = 0.5

f(x) = sin(x^2)

(f(x+ϵ) - f(x)) / ϵ
#-
2x*cos(x^2)

# So complex numbers can be used to get exact derivatives! This is very efficient
# and can be written using only one call to `f`.

imag(f(x+ϵ)) / imag(ϵ)

# Another way of looking at this is that we are doing bog-standard numerical
# differentiation, but the use of complex numbers avoids the typical problem
# with that technique (i.e. that a small perturbation ends up being overwhelmed
# by floating point error). The dual number is then a slight variation which
# makes the limit $\epsilon \rightarrow 0$ exact, rather than approximate.
# Forward mode AD can be described as "just" a clever implementation of
# numerical differentiation. Both numerical and forward derivatives propagate
# a perturbation $\epsilon$ using the same basic rules, and they have the
# same algorithmic properties.


text:
module Kdb
import Q._k: k, khp
const OS = Dict(:Darwin => 'm', :Linux => 'l', :Windows => 'w')
const QHOME = get(Dict(ENV), "QHOME") do
    joinpath(homedir(), "q")
end
const STARTUP_CODE = """
p:1024
while[0~@[system;"p ",string p;0];p+:1]
-1 string p
\\1 /dev/null
"""

"""
Returns the path to the kdb+ binary.
"""
function kdb_binary()
    os = OS[Sys.KERNEL]
    binary = joinpath(QHOME, "$(os)64", "q")
    if isfile(binary)
        return binary
    end
    joinpath(QHOME, "$(os)32", "q")
end

q_command() = `$(kdb_binary())`

"""
Start a slave kdb+ process and return a comm handle.
"""
function start()
    p = spawn(q_command(), (Pipe(), Pipe(), STDERR))
    write(p.in, STARTUP_CODE)
    port = readline(p.out)
    close(p.out)
    parse(Int, port), p
end

function stop(handle, process)
    k(handle, "exit 0")
    wait(process)
    process.exitcode
end

function __init__()
    global QBIN
    QBIN = kdb_binary()
end
end  # module Kdb


text:
using Stheno: GPC, EQ

@timedtestset "product" begin
    @timedtestset "GP mul errors" begin
        gpc = GPC()
        f1, f2 = GP(EQ(), gpc), GP(EQ(), gpc)
        @test_throws ArgumentError f1 * f2
    end
    @timedtestset "multiply by constant" begin
        rng, N, N′, D = MersenneTwister(123456), 3, 5, 2
        X, X′ = ColVecs(randn(rng, D, N)), ColVecs(randn(rng, D, N′))
        g1, c, c′ = GP(1, EQ(), GPC()), -4.3, 2.1
        g2, g2′ = c * g1, g1 * c′
        g3, g3′ = c * g2, g2′ * c′
        g4, g4′ = c * g3, g3′ * c′

        @test mean(g2(X)) == c .* mean(g1(X))
        @test mean(g3(X)) == c .* mean(g2(X))
        @test mean(g4(X)) == c .* mean(g3(X))
        @test mean(g2′(X)) == mean(g1(X)) .* c′
        @test mean(g3′(X)) == mean(g2′(X)) .* c′
        @test mean(g4′(X)) == mean(g3′(X)) .* c′

        @test cov(g2(X)) ≈ c^2 .* cov(g1(X))
        @test cov(g3(X)) ≈ c^2 .* cov(g2(X))
        @test cov(g4(X)) ≈ c^2 .* cov(g3(X))
        @test cov(g2′(X)) ≈ c′^2 .* cov(g1(X))
        @test cov(g3′(X)) ≈ c′^2 .* cov(g2′(X))
        @test cov(g4′(X)) ≈ c′^2 .* cov(g3′(X))

        @test cov(g2(X), g1(X′)) ≈ c .* cov(g1(X), g1(X′))
        @test cov(g3(X), g1(X′)) ≈ c^2 .* cov(g1(X), g1(X′))
        @test cov(g4(X), g1(X′)) ≈ c^3 .* cov(g1(X), g1(X′))
        @test cov(g2′(X), g1(X′)) ≈ c′ .* cov(g1(X), g1(X′))
        @test cov(g3′(X), g1(X′)) ≈ c′^2 .* cov(g1(X), g1(X′))
        @test cov(g4′(X), g1(X′)) ≈ c′^3 .* cov(g1(X), g1(X′))

        @test cov(g2(X), g2′(X′)) ≈ (c * c′) .* cov(g1(X), g1(X′))
        @test cov(g3(X), g3′(X′)) ≈ (c * c′)^2 .* cov(g1(X), g1(X′))
        @test cov(g4(X), g4′(X′)) ≈ (c * c′)^3 .* cov(g1(X), g1(X′))

        @test cov(g2(X), g3′(X′)) ≈ (c * c′^2) .* cov(g1(X), g1(X′))
        @test cov(g3(X), g2′(X′)) ≈ (c^2 * c′) .* cov(g1(X), g1(X′))
        @test cov(g2(X), g4′(X′)) ≈ (c * c′^3) .* cov(g1(X), g1(X′))
        @test cov(g4(X), g2′(X′)) ≈ (c^3 * c′) .* cov(g1(X), g1(X′))

        @timedtestset "Consistency Tests" begin
            rng, P, Q = MersenneTwister(123456), 3, 5
            x0 = collect(range(-1.0, 1.0; length=P))
            x1 = collect(range(-0.5, 1.5; length=Q))
            x2, x3 = randn(rng, Q), randn(rng, P)

            gpc = GPC()
            f1 = GP(cos, EQ(), gpc)
            f2 = 5 * f1
            abstractgp_interface_tests(f2, f1, x0, x1, x2, x3)
        end
        @timedtestset "Diff Tests" begin
            standard_1D_tests(
                MersenneTwister(123456),
                Dict(:σ=>2.3),
                θ->begin
                    f = GP(0.5, EQ(), GPC())
                    return θ[:σ] * f, f
                end,
                X, X′,
            )
        end
    end
    @timedtestset "multiply by function" begin
        rng, N, N′, D = MersenneTwister(123456), 3, 5, 2
        X, X′ = ColVecs(randn(rng, D, N)), ColVecs(randn(rng, D, N′))
        g1, f, f′ = GP(1, EQ(), GPC()), x->sum(sin, x), x->sum(cos, x)
        g2, g2′ = f * g1, g1 * f′
        g3, g3′ = f * g2, g2′ * f′
        g4, g4′ = f * g3, g3′ * f′

        @test mean(g2(X)) == map(f, X) .* mean(g1(X))
        @test mean(g3(X)) == map(f, X) .* mean(g2(X))
        @test mean(g4(X)) == map(f, X) .* mean(g3(X))
        @test mean(g2′(X)) == mean(g1(X)) .* map(f′, X)
        @test mean(g3′(X)) == mean(g2′(X)) .* map(f′, X)
        @test mean(g4′(X)) == mean(g3′(X)) .* map(f′, X)

        fX, f′X = map(f, X), map(f′, X)
        @test cov(g2(X)) ≈ fX .* cov(g1(X)) .* fX'
        @test cov(g3(X)) ≈ fX .* cov(g2(X)) .* fX'
        @test cov(g4(X)) ≈ fX .* cov(g3(X)) .* fX'
        @test cov(g2′(X)) ≈ f′X .* cov(g1(X)) .* f′X'
        @test cov(g3′(X)) ≈ f′X .* cov(g2′(X)) .* f′X'
        @test cov(g4′(X)) ≈ f′X .* cov(g3′(X)) .* f′X'

        fX′, f′X′ = map(f, X′), map(f′, X′)
        @test cov(g2(X), g1(X′)) ≈ fX .* cov(g1(X), g1(X′))
        @test cov(g3(X), g1(X′)) ≈ fX.^2 .* cov(g1(X), g1(X′))
        @test cov(g4(X), g1(X′)) ≈ fX.^3 .* cov(g1(X), g1(X′))
        @test cov(g2′(X), g1(X′)) ≈ f′X .* cov(g1(X), g1(X′))
        @test cov(g3′(X), g1(X′)) ≈ (f′X.^2) .* cov(g1(X), g1(X′))
        @test cov(g4′(X), g1(X′)) ≈ (f′X.^3) .* cov(g1(X), g1(X′))

        @test cov(g1(X′), g2(X)) == cov(g2(X), g1(X′))'
        @test cov(g1(X′), g3(X)) == cov(g3(X), g1(X′))'
        @test cov(g1(X′), g4(X)) == cov(g4(X), g1(X′))'

        @test cov(g2(X), g2′(X′)) ≈ fX .* cov(g1(X), g1(X′)) .* f′X′'
        @test cov(g3(X), g3′(X′)) ≈ fX.^2 .* cov(g1(X), g1(X′)) .* (f′X′.^2)'
        @test cov(g4(X), g4′(X′)) ≈ fX.^3 .* cov(g1(X), g1(X′)) .* (f′X′.^3)'

        @test cov(g2(X), g3′(X′)) ≈ fX .* cov(g1(X), g1(X′)) .* (f′X′').^2
        @test cov(g3(X), g2′(X′)) ≈ fX.^2 .* cov(g1(X), g1(X′)) .* (f′X′')
        @test cov(g2(X), g4′(X′)) ≈ fX .* cov(g1(X), g1(X′)) .* (f′X′').^3
        @test cov(g4(X), g2′(X′)) ≈ fX.^3 .* cov(g1(X), g1(X′)) .* (f′X′')

        @timedtestset "Consistency Tests" begin
            rng, P, Q = MersenneTwister(123456), 3, 5
            x0 = collect(range(-1.0, 1.0; length=P))
            x1 = collect(range(-0.5, 1.5; length=Q))
            x2, x3 = randn(rng, Q), randn(rng, P)

            gpc = GPC()
            f1 = GP(cos, EQ(), gpc)
            f2 = sin * f1
            abstractgp_interface_tests(f2, f1, x0, x1, x2, x3)
        end
        @timedtestset "Diff Tests" begin
            standard_1D_tests(
                MersenneTwister(123456),
                Dict(:σ=>2.3, :c=>1.3),
                θ->begin
                    f = GP(θ[:c], EQ(), GPC())
                    return (x->θ[:σ] * x) * f, f
                end,
                collect(range(-2.0, 2.0; length=N)),
                collect(range(-1.5, 2.5; length=N′)),
            )
        end
    end
end


text:
# This file was generated, do not modify it. # hide
@load FillImputer
filler = machine(FillImputer(), datac)
fit!(filler)
datac = transform(filler, datac)

y, X = unpack(datac, ==(:outcome), name->true);
X = coerce(X, autotype(X, :discrete_to_continuous));


text:
println("Imports")
using AdventOfCode2018
using AdventOfCode2018.Day5

println("getting input")
input = strip(get_aoc_input(5))

println("start stuff")
println("Solution 1: ", part_1(input))
println("Solution 2: ", part_2(input))


text:
# Add cadd to cbase
function add_constraint!(cbase::Grid2DConstraints, cadd::Grid2DConstraints)

    union!(cbase.vertex_constraints, cadd.vertex_constraints)
    union!(cbase.edge_constraints, cadd.edge_constraints)

end

function overlap_between_constraints(cbase::Grid2DConstraints, cother::Grid2DConstraints)

    vertex_intersect = intersect(cbase.vertex_constraints, cother.vertex_constraints)
    edge_intersect = intersect(cbase.edge_constraints, cother.edge_constraints)

    return ~(isempty(vertex_intersect)) || ~(isempty(edge_intersect))

end

# Not setting agent_idx and constraints as they will always be there for functions
@with_kw mutable struct Grid2DEnvironment <: MAPFEnvironment
    dimx::Int64
    dimy::Int64
    obstacles::Set{Grid2DLocation}
    goals::Vector{Grid2DLocation}
    state_graph::SimpleVListGraph{Grid2DState}    = SimpleVListGraph{Grid2DState}()
    state_to_idx::Dict{Grid2DState,Int64}         = Dict{Grid2DState,Int64}()
    last_goal_constraint::Int64             = -1
    agent_idx::Int64                        = 0
    curr_goal_idx::Int64                    = 0
end

function get_mapf_state_from_idx(env::Grid2DEnvironment, idx::Int64)
    return env.state_graph.vertices[idx]
end

function get_mapf_action(env::Grid2DEnvironment, source::Int64, target::Int64)

    state1 = env.state_graph.vertices[source]
    state2 = env.state_graph.vertices[target]

    if equal_except_time(state1, state2)
        return Grid2DAction(Wait::Action)
    else
        if state2.x > state1.x
            return Grid2DAction(Right::Action)
        elseif state2.x < state1.x
            return Grid2DAction(Left::Action)
        elseif state2.y > state1.y
            return Grid2DAction(Up::Action)
        else
            return Grid2DAction(Down::Action)
        end
    end
    return nothing
end

function set_low_level_context!(env::Grid2DEnvironment, agent_idx::Int64, constraints::Grid2DConstraints)

    env.last_goal_constraint = -1
    env.agent_idx = agent_idx

    for vc in constraints.vertex_constraints
        if (vc.x == env.goals[agent_idx].x && vc.y == env.goals[agent_idx].y)
            env.last_goal_constraint = max(env.last_goal_constraint, vc.time)
        end
    end
end

function admissible_heuristic_grid2d(env::Grid2DEnvironment, s::Grid2DState)
    return abs(s.x - env.goals[env.agent_idx].x) + abs(s.y - env.goals[env.agent_idx].y)
end

function is_solution(env::Grid2DEnvironment, s::Grid2DState)
    return s.x == env.goals[env.agent_idx].x && s.y == env.goals[env.agent_idx].y &&
            s.time > env.last_goal_constraint
end

function get_state(agent_idx::Int64, solution::Vector{PlanResult{Grid2DState,Grid2DAction,Int64}}, rel_time::Int64)

    @assert agent_idx <= length(solution)

    # t must be relative time index from start of solution
    if rel_time < length(solution[agent_idx].states)
        return solution[agent_idx].states[rel_time+1][1]
    end

    # Just return the last state after verifying it is not empty
    @assert ~(isempty(solution[agent_idx].states))
    return solution[agent_idx].states[end][1]

end

function state_valid(env::Grid2DEnvironment, constraints::Grid2DConstraints, s::Grid2DState)

    con = constraints.vertex_constraints

    return s.x >=0 && s.x < env.dimx && s.y >=0 && s.y < env.dimy &&
            ~(Grid2DLocation((x=s.x, y=s.y)) in env.obstacles) &&
            ~(VertexConstraint(time=s.time, x=s.x, y=s.y) in con)
end

function transition_valid(env::Grid2DEnvironment, constraints::Grid2DConstraints, s1::Grid2DState, s2::Grid2DState)

    con = constraints.edge_constraints

    return ~(EdgeConstraint(time=s1.time, x1=s1.x, y1=s1.y, x2=s2.x, y2=s2.y) in con)
end

function get_env_state_idx!(env::Grid2DEnvironment, s::Grid2DState)

    idx = get(env.state_to_idx, s, 0)

    if idx == 0
        add_vertex!(env.state_graph, s)
        nbr_idx = num_vertices(env.state_graph)
        env.state_to_idx[s] = nbr_idx
        idx = nbr_idx
    end

    return idx
end

function get_first_conflict(env::Grid2DEnvironment, solution::Vector{PlanResult{Grid2DState,Grid2DAction,Int64}})

    max_rel_time = 0
    for sol in solution
        max_rel_time = max(max_rel_time, length(sol.states) - 1)
    end

    for rel_time = 0:max_rel_time-1
        # Drive-drive vertex collisions
        for i = 1:length(solution)-1
            state1 = get_state(i, solution, rel_time)
            for j = i+1:length(solution)
                state2 = get_state(j, solution, rel_time)

                if equal_except_time(state1, state2)
                    return Grid2DConflict(time=rel_time, agent_1=i, agent_2=j,
                                    type=Vertex::ConflictType, x1=state1.x, y1=state1.y)
                end
            end
        end

        # Drive-drive edge
        for i = 1:length(solution)-1
            state1a = get_state(i, solution, rel_time)
            state1b = get_state(i, solution, rel_time + 1)

            for j = i+1:length(solution)
                state2a = get_state(j, solution, rel_time)
                state2b = get_state(j, solution, rel_time + 1)

                if equal_except_time(state1a, state2b) &&
                    equal_except_time(state1b, state2a)

                    return Grid2DConflict(time=rel_time, agent_1=i, agent_2=j,
                                    type=Edge::ConflictType, x1=state1a.x, y1=state1a.y,
                                    x2=state1b.x, y2=state1b.y)
                end
            end
        end
    end

    # No conflict
    return nothing
end

# Return a Dict{Int64, Grid2DConstraints}
function create_constraints_from_conflict(env::Grid2DEnvironment, conflict::Grid2DConflict)

    res_constraints = Dict{Int64, Grid2DConstraints}()

    if conflict.type == Vertex::ConflictType

        c1 = Grid2DConstraints()
        push!(c1.vertex_constraints, VertexConstraint(time=conflict.time, x=conflict.x1, y=conflict.y1))
        res_constraints[conflict.agent_1] = c1
        res_constraints[conflict.agent_2] = c1

    elseif conflict.type == Edge::ConflictType

        c1 = Grid2DConstraints()
        push!(c1.edge_constraints, EdgeConstraint(time=conflict.time, x1=conflict.x1, y1=conflict.y1,
                                                  x2=conflict.x2, y2=conflict.y2))
        res_constraints[conflict.agent_1] = c1

        c2 = Grid2DConstraints()
        push!(c2.edge_constraints, EdgeConstraint(time=conflict.time, x1=conflict.x2, y1=conflict.y2,
                                                  x2=conflict.x1, y2=conflict.y1))
        res_constraints[conflict.agent_2] = c2
    end

    return [res_constraints]
end

function low_level_search!(solver::CBSSolver, agent_idx::Int64, s::Grid2DState, constraints::Grid2DConstraints)

    env = solver.env

    # Reset env index
    env.curr_goal_idx = 0

    # Retrieve vertex index for state
    idx = get_env_state_idx!(env, s)

    # Set all edge weights to one
    edge_wt_fn(u, v) = 1

    # Set the heuristic
    heuristic(v) = admissible_heuristic_grid2d(env, v)

    # Run the search
    # @info "RUNNING SEARCH!"
    vis = CBSGoalVisitorImplicit(env, constraints)
    a_star_states =  a_star_implicit_shortest_path!(env.state_graph, edge_wt_fn,
                                                    idx, vis,
                                                    heuristic, Int64)

    plan_result = get_plan_result_from_astar(env, a_star_states.dists, a_star_states.parent_indices, idx, env.curr_goal_idx, 0)

    # Return empty solution
    if plan_result == nothing
        return PlanResult{Grid2DState,Grid2DAction,Int64}()
    end

    return plan_result
end

mutable struct CBSGoalVisitorImplicit <: AbstractDijkstraVisitor
    env::Grid2DEnvironment
    constraints::Grid2DConstraints
end

function Graphs.include_vertex!(vis::CBSGoalVisitorImplicit, u::Grid2DState, v::Grid2DState,
                                d::N, nbrs::Vector{Int64}) where {N <: Number}

    if is_solution(vis.env, v)
        # @info "Found low-level solution!"
        vis.env.curr_goal_idx = get_env_state_idx!(vis.env, v)
        return false
    end

    # Need to generate neighbors of v using Actions
    empty!(nbrs)

    new_time = v.time + 1

    # Action = Wait
    temp_state = Grid2DState(time=new_time, x=v.x, y=v.y)
    if state_valid(vis.env, vis.constraints, temp_state) &&
        transition_valid(vis.env, vis.constraints, v, temp_state)
        nbr_idx = get_env_state_idx!(vis.env, temp_state)
        push!(nbrs, nbr_idx)
    end

    # Action = Left
    temp_state = Grid2DState(time=new_time, x=v.x-1, y=v.y)
    if state_valid(vis.env, vis.constraints, temp_state) &&
        transition_valid(vis.env, vis.constraints, v, temp_state)
        nbr_idx = get_env_state_idx!(vis.env, temp_state)
        push!(nbrs, nbr_idx)
    end

    # Action = Right
    temp_state = Grid2DState(time=new_time, x=v.x+1, y=v.y)
    if state_valid(vis.env, vis.constraints, temp_state) &&
        transition_valid(vis.env, vis.constraints, v, temp_state)
        nbr_idx = get_env_state_idx!(vis.env, temp_state)
        push!(nbrs, nbr_idx)
    end

    # Action = Up
    temp_state = Grid2DState(time=new_time, x=v.x, y=v.y+1)
    if state_valid(vis.env, vis.constraints, temp_state) &&
        transition_valid(vis.env, vis.constraints, v, temp_state)
        nbr_idx = get_env_state_idx!(vis.env, temp_state)
        push!(nbrs, nbr_idx)
    end

    # Action = Down
    temp_state = Grid2DState(time=new_time, x=v.x, y=v.y-1)
    if state_valid(vis.env, vis.constraints, temp_state) &&
        transition_valid(vis.env, vis.constraints, v, temp_state)
        nbr_idx = get_env_state_idx!(vis.env, temp_state)
        push!(nbrs, nbr_idx)
    end

    # @show nbrs
    return true
end


text:
"""
Table Schema generic data structure
https://github.com/frictionlessdata/TableSchema.jl#table
"""
mutable struct Table
    source
    headers::Array{String}
    schema::Schema
    errors::Array{ConstraintError}

    function Table(csvfilename::String, schema::Schema=Schema() ; encoding::Encoding=enc"UTF-8")
        if match(r"^https?://", csvfilename) !== nothing
            source = read_remote_csv(csvfilename)
        else
			if encoding === enc"UTF-8"
				source = readdlm(csvfilename, DELIMITER)
			else
				io = StringEncodings.open(csvfilename, encoding)
				source = readdlm(io, DELIMITER)
				close(io)
			end
        end
        headers, source = get_headers(source)
        new(source, headers, schema, [])
    end
    function Table(csvdata::Base.GenericIOBuffer, schema::Schema=Schema())
        source = readdlm(csvdata, DELIMITER)
        headers, source = get_headers(source)
        new(source, headers, schema, [])
    end
    function Table(source, headers::Array{String}, schema::Schema=Schema())
        new(source, headers, schema, [])
    end
    function Table()
        new(Nothing, [], Schema(), [])
    end
end

function read_remote_csv(url::String)
    req = request("GET", url)
    readdlm(req.body, DELIMITER)
end

function get_headers(source::Array)
    headers = [ String(s) for s in source[1,:] ]
    headers, source[2:end,:] # clear the headers
end

function read(t::Table ; data=nothing, keyed=false, extended=false, cast=true, relations=false, limit=nothing)
    (keyed == false && extended == false && relations == false && limit == nothing) || throw(ErrorException("Not implemented"))
    if data != nothing
        if typeof(data) == Array{Any, 2}
            t.headers, t.source = get_headers(data)
        else
            throw(ErrorException("Data must be a 2-dimensional array"))
        end
    end
    if cast
        if !is_valid(t.schema)
            throw(ErrorException("Schema must be valid to cast Table"))
        end
        if t.source == nothing
            throw(ErrorException("Data must be available to cast Table"))
        end
        newtable = Nothing
		# Iterate table rows
        for row in t
			# Apply cast to the row's values
            newrow = cast_row(t.schema, row, false, false)
			# Reshape back into a non-elemental array
			newrow = reshape(newrow, 1, length(newrow))
            if newtable == Nothing
                newtable = newrow
            else
                newtable = vcat(newtable, newrow)
            end
        end
        t.source = newtable
    end
    t.source
end

function infer(t::Table ; limit::Int64=-1)
	# what
	limit !== -1 && throw(ErrorException("limit parameter not implemented"))
	tr = read(t, cast=false)
	t.schema = Schema()
	infer(t.schema, tr, t.headers)
end

function save(t::Table, target::String)
	# TODO: should we check the schema too?
	# !valid(t.schema) &&
    #     throw(TableValidationException("Schema not valid"))
	@debug "Building table data"
	headers = reshape(t.headers, 1, length(t.headers))
	tabledata = vcat(headers, t.source)
	@debug "Saving table data to $target"
	delim = ','
	open(target, "w") do io
		writedlm(io, tabledata, delim)
	end
end

function validate(t::Table)
    is_empty(t.schema) &&
        throw(TableValidationException("No schema available"))
    # TODO: should we check the schema too?
    # !valid(t.schema) &&
    #     throw(TableValidationException("Schema not valid"))
    tr = t.source
    for fld in t.schema.fields
        ix = findall(in([fld.name]), t.headers)
        if length(ix) != 1
            # TODO: shouldn't this just cause a ConstraintError?
            throw(TableValidationException("Missing field defined in Schema: $(fld.name)"))
        end
        try
            column = tr[:,ix]
            for r = 1:size(tr, 2)
                row = tr[r,:]
                checkrow(fld, row[ix[1]], column)
            end
        catch ex
            if isa(ex, ConstraintError)
                push!(t.errors, ex)
            end
        end
    end
    return length(t.errors) == 0
end

Base.eltype(it::Table) = Table
Base.length(it::Table) = size(it.source, 1)
function Base.iterate(it::Table, (el, i)=(nothing, 1))
    if i > length(it); return nothing; end
    return (it.source[i,:], (nothing, i + 1))
end

# Base.start(t::Table) = 1
# Base.done(t::Table, i) = i > size(t.source, 1)
# Base.next(t::Table, i) = t.source[i,:], i+1

max_stars_repo_name: loleg/tableschema-jl

max_issues_repo_name: loleg/tableschema-jl

max_forks_repo_name: loleg/tableschema-jl


text:
module Validation

using SearchLight

export ValidationResult, ValidationError, ValidationRule, ModelValidator

const valid = true
const invalid = false
export valid, invalid

export validator, validate, haserrors, haserrorsfor, errorsfor, errors_messages_for, errors_to_string

abstract type ValidationAbstractType end

struct ValidationResult <: ValidationAbstractType
  validation_status::Bool
  error_type::Symbol
  error_message::String

  ValidationResult(validation_status, error_type, error_message) = new(validation_status, error_type, error_message)
end
ValidationResult(; validation_status = false, error_type = :validation_error, validation_message = "Model's data is invalid") = ValidationResult(validation_status, error_type, validation_message)
ValidationResult(validation_status::Bool) = validation_status ? ValidationResult(true, :no_error, "") : ValidationResult(validation_status = false)

struct ValidationError <: ValidationAbstractType
  field::Symbol
  error_type::Symbol
  error_message::String
end

struct ValidationRule <: ValidationAbstractType
  field::Symbol
  validator_function::Function
  validator_arguments::Tuple

  ValidationRule(field, validator_function, validator_arguments = ()) = new(field, validator_function, validator_arguments)
end

"""
The object that defines the rules and stores the validation errors associated with the fields of a `model`.
"""
struct ModelValidator
  rules::Vector{ValidationRule}
  errors::Vector{ValidationError}
end
ModelValidator(rules::Vector{ValidationRule}) = ModelValidator(rules, ValidationError[])

# overwrite!
function validator(m::Type{T})::ModelValidator where {T<:SearchLight.AbstractModel}
  ModelValidator(ValidationRule[])
end

#
# validation API
#

"""
    validate(m::T)::Bool where {T<:AbstractModel}

Validates `m`'s data. A `bool` is return and existing errors are pushed to the validator's error stack.
"""
function validate(m::T)::ModelValidator where {T<:SearchLight.AbstractModel}
  mv = validator(typeof(m))

  for r in mv.rules
    vr::ValidationResult = r.validator_function(r.field, m, r.validator_arguments...)
    vr.validation_status || push!(mv.errors, ValidationError(r.field, vr.error_type, vr.error_message))
  end

  mv
end

function haserrors(mv::ModelValidator) :: Bool
  ! isempty(mv.errors)
end

function haserrorsfor(mv::ModelValidator, field::Symbol) :: Bool
  ! isempty(errorsfor(mv, field))
end

const has_errors_for = haserrorsfor # for Genie compatibility

function errorsfor(mv::ModelValidator, field::Union{Symbol,Nothing} = nothing) :: Vector{ValidationError}
  result = ValidationError[]

  for err in mv.errors
    (field === nothing || err.field == field) && push!(result, err)
  end

  result
end

function errorsmessagesfor(mv::ModelValidator, field::Union{Symbol,Nothing} = nothing) :: Vector{String}
  result = String[]

  for err in errorsfor(mv, field)
    push!(result, err.error_message)
  end

  result
end

function errors_to_string(mv::ModelValidator, field::Union{Symbol,Nothing} = nothing;
                          separator::String = "\n", prepend_field::Bool = true,
                          uppercase_first::Bool = prepend_field) :: String
  errors = String[]

  for err in mv.errors
    (field === nothing || err.field == field) && push!(errors, (prepend_field ? "$(err.field) " : "") * err.error_message)
  end

  join( map(x ->  (uppercase_first ? uppercasefirst(x) : x), errors), separator)
end

end


text:
# TODO: interface should be changed so that we can control deformula aprameters
function phfit(f::Any, cf1::CF1{Tv}, bounds::Tuple{Tv,Tv} = (Tv(0), Tv(Inf)), ::Type{MatT} = SparseMatrixCSC;
    initialize = true, eps::Tv = Tv(1.0e-8), ufact::Tv = Tv(1.01), cf1sort=true, verbose = [false, false],
    steps = 50,
    shapes = Tv[1, 4, 16, 64, 256, 1024], scales = Tv[0.5, 1.0, 2.0], init_maxiter = 5,
    maxiter = 5000, reltol = Tv(1.0e-8)) where {Tv,MatT}

    data = WeightedSample(f, bounds)
    phfit(cf1, data, MatT, initialize=initialize, eps=eps, ufact=ufact, cf1sort=cf1sort, verbose=verbose, shapes=shapes,
        scales=scales, init_maxiter=init_maxiter, steps=steps, maxiter=maxiter, reltol=reltol)
end

function phfit(cf1::CF1{Tv}, data::AbstractPHSample, ::Type{MatT} = SparseMatrixCSC;
    initialize = true, eps::Tv = Tv(1.0e-8), ufact::Tv = Tv(1.01), cf1sort=true, verbose = [false, false],
    steps = 50,
    shapes = Tv[1, 4, 16, 64, 256, 1024], scales = Tv[0.5, 1.0, 2.0], init_maxiter = 5,
    maxiter = 5000, reltol = Tv(1.0e-8)) where {Tv,MatT}

    local newcf1::CF1{Tv}
    if initialize
        newcf1 = initializePH(cf1, data, MatT, shapes=shapes, scales=scales,
            maxiter=init_maxiter, verbose = verbose[1], eps=eps, ufact=ufact)
    else
        newcf1 = copy(cf1)
    end
    llf, conv, iter, rerror = phfit!(newcf1, data, MatT, eps=eps, ufact=ufact,
        cf1sort=cf1sort, verbose = verbose[2], steps=steps, maxiter=maxiter, reltol=reltol)
    return newcf1, llf, conv, iter, rerror
end

function phfit!(cf1::CF1{Tv}, data::AbstractPHSample, ::Type{MatT} = SparseMatrixCSC;
    eps::Tv = Tv(1.0e-8), ufact::Tv = Tv(1.01), cf1sort = true, verbose = false,
    steps = 50, maxiter = 5000, reltol = Tv(1.0e-8)) where {Tv,MatT}
    iter = 0
    ph = GPH(cf1, MatT)
    eres = Estep(ph)
    conv = false
    local rerror::Tv

    llf = estep!(cf1, data, eres, eps=eps, ufact=ufact)
    mstep!(cf1, eres)
    while true
        prevllf = llf
        for k = 1:steps
            llf = estep!(cf1, data, eres, eps=eps, ufact=ufact)
            mstep!(cf1, eres, cf1sort=cf1sort)
        end
        iter += steps
        rerror = abs((llf - prevllf) / prevllf)
        verbose && println("iter=$(iter) llf=$(llf) rerror=$(rerror)")

        if llf < prevllf
            println("Warning: llf does not increase at iter=$(iter)")
        end

        if rerror < reltol
            conv = true
            break
        end

        if iter >= maxiter
            break
        end
    end
    return llf, conv, iter, rerror
end


text:
using DecisionTree
using ScikitLearn
using PyPlot

# Create a random dataset
srand(100)
X = sort(5 * rand(80))
XX = reshape(X, 80, 1)
y = sin(X)
y[1:5:end] += 3 * (0.5 - rand(16))

# Fit regression model
regr_1 = DecisionTreeRegressor()
regr_2 = DecisionTreeRegressor(pruning_purity_threshold=0.05)

fit!(regr_1, XX, y)

fit!(regr_2, XX, y)

# Predict
X_test = 0:0.01:5.0
y_1 = predict(regr_1, hcat(X_test))
y_2 = predict(regr_2, hcat(X_test))

# Plot the results
scatter(X, y, c="k", label="data")
plot(X_test, y_1, c="g", label="no pruning", linewidth=2)
plot(X_test, y_2, c="r", label="pruning_purity_threshold=0.05", linewidth=2)

xlabel("data")
ylabel("target")
title("Decision Tree Regression")
legend(prop=Dict("size"=>10))

max_stars_repo_name: PacktPublishing/Julia-for-data-science

max_issues_repo_name: PacktPublishing/Julia-for-data-science

max_forks_repo_name: PacktPublishing/Julia-for-data-science


text:
"""
```
augment_states{T<:AbstractFloat}(m::AbstractDSGEModel, TTT::Matrix{T}, RRR::Matrix{T}, CCC::Matrix{T})
```

### Arguments

-`m`: the m object
-`TTT`, `RRR`, and `CCC`: matrices of the state transition equation

### Return values

- `TTT_aug`, `RRR_aug`, and `CCC_aug`: extend the corresponding input matrices to include
  jobservables which are growth rates.

### Description

Some observables in the model are growth rates, which are calculated as a linear combination
of a present and lagged state (which is not yet accounted for in the `TTT`, `RRR`,and `CCC`
matrices). To improve the performance of `gensys`, these additional states are added after
the model is solved. `augment_states` assigns an index to each lagged state, and extends the
input `TTT`, `RRR`, and `CCC` matrices to accommodate the additional states and capture the
lagged state value in the current state vector. `RRR` and `CCC` are mostly augmented with
zeros.

The diagram below shows how `TTT` is extended to `TTT_aug`.

                TTT_aug
     (m.endogenous_states_additional
                  x
     m.endogenous_states_additional)
     _________________________________
    |                     |           |
    |          TTT        | endog_    |
    | (endogenous_states  | states_   |
    |          x          | augmented |
    |  endogenous_states) |           |
    |_____________________|           |
    |                                 |
    |    endogenous_states_augmented  |
    |_________________________________|

"""
function augment_states(m::Model805{T}, TTT::Matrix{T}, RRR::Matrix{T}, CCC::Vector{T}) where {T<:AbstractFloat}
    endo = m.endogenous_states
    endo_addl = m.endogenous_states_augmented
    exo = m.exogenous_shocks

    n_endo = n_states(m)
    n_exo = n_shocks_exogenous(m)
    @assert (n_endo, n_endo) == size(TTT)
    @assert (n_endo, n_exo) == size(RRR)
    @assert n_endo == length(CCC)

    # Initialize augmented matrices
    n_addl_states = length(m.endogenous_states_augmented)
    n_addl_eqs = n_addl_states
    TTT_aug = zeros(n_endo + n_addl_eqs, n_endo + n_addl_states)
    TTT_aug[1:n_endo, 1:n_endo] = TTT
    RRR_aug = [RRR; zeros(n_addl_eqs, n_exo)]
    CCC_aug = [CCC; zeros(n_addl_eqs)]

    ### TTT modifications

    # Track Lags
    TTT_aug[endo_addl[:y_t1], endo[:y_t]] = 1.0
    TTT_aug[endo_addl[:c_t1], endo[:c_t]] = 1.0
    TTT_aug[endo_addl[:i_t1], endo[:i_t]] = 1.0
    TTT_aug[endo_addl[:w_t1], endo[:w_t]] = 1.0
    TTT_aug[endo_addl[:π_t1], endo[:π_t]] = 1.0
    TTT_aug[endo_addl[:L_t1], endo[:L_t]]  = 1.0

    ## We construct state for expected inflation using the fact
    ##
    ##   E_t[s_{t+1}] = CCC + TTT*s_{t} = CCC + TTT * (CCC + TTT*s_{t-1} + RRR)
    ##                = (CCC+TTT*CCC) + (TTT^2)s_{t-1} + TTT*RRR
    ##
    ## So to construct the state for E_t[p_{t+1}], we need to make all these
    ## objects, and just index out the row relevant to pi_t

    T2 = TTT^2
    TR = TTT*RRR
    CTC = CCC+TTT*CCC

    TTT_aug[endo_addl[:Et_π_t],:] = [T2[endo[:π_t],:]' zeros(n_addl_states)'];

    RRR_aug[endo_addl[:Et_π_t],:] = TR[endo[:π_t],:]
    CCC_aug[endo_addl[:Et_π_t],:] = CTC[endo[:π_t],:]

    return TTT_aug, RRR_aug, CCC_aug
end


text:
"""
    ContinuousDEQSolver(alg=VCABM4(); mode::Symbol=:rel_deq_default, abstol=1e-8, reltol=1e-8, tspan=Inf)

Solver for Continuous DEQ Problem ([pal2022mixing](@cite)). Similar to `DynamicSS` but provides more flexibility needed
for solving DEQ problems.

## Arguments

* `alg`: Algorithm to solve the ODEProblem. (Default: `VCABM4()`)
* `mode`: Termination Mode of the solver. See below for a description of the various termination conditions (Default: `:rel_deq_default`)
* `abstol`: Absolute tolerance for termination. (Default: `1e-8`)
* `reltol`: Relative tolerance for termination. (Default: `1e-8`)
* `tspan`: Time span. Users should not change this value, instead control termination through `maxiters` in `solve` (Default: `Inf`)

## Termination Modes

#### Termination on Absolute Tolerance

* `:abs`: Terminates if ``all \\left( | \\frac{\\partial u}{\\partial t} | \\leq abstol \\right)``
* `:abs_norm`: Terminates if ``\\| \\frac{\\partial u}{\\partial t} \\| \\leq abstol``
* `:abs_deq_default`: Essentially `abs_norm` + terminate if there has been no improvement for the last 30 steps + terminate if the solution blows up (diverges)
* `:abs_deq_best`: Same as `:abs_deq_default` but uses the best solution found so far, i.e. deviates only if the solution has not converged

#### Termination on Relative Tolerance

* `:rel`: Terminates if ``all \\left(| \\frac{\\partial u}{\\partial t} | \\leq reltol \\times | u | \\right)``
* `:rel_norm`: Terminates if ``\\| \\frac{\\partial u}{\\partial t} \\| \\leq reltol \\times \\| \\frac{\\partial u}{\\partial t} + u \\|``
* `:rel_deq_default`: Essentially `rel_norm` + terminate if there has been no improvement for the last 30 steps + terminate if the solution blows up (diverges)
* `:rel_deq_best`: Same as `:rel_deq_default` but uses the best solution found so far, i.e. deviates only if the solution has not converged

#### Termination using both Absolute and Relative Tolerances

* `:norm`: Terminates if ``\\| \\frac{\\partial u}{\\partial t} \\| \\leq reltol \\times \\| \\frac{\\partial u}{\\partial t} + u \\|`` &
           ``\\| \\frac{\\partial u}{\\partial t} \\| \\leq abstol``
* `fallback`: Check if all values of the derivative is close to zero wrt both relative and absolute tolerance. This is usable for small problems
              but doesn't scale well for neural networks, and should be avoided unless absolutely necessary

See also: [`DiscreteDEQSolver`](@ref)

!!! note
    This  will be upstreamed to DiffEqSensitivity in the later releases of the package
"""
struct ContinuousDEQSolver{M,A,AT,RT,TS} <: SteadyStateDiffEq.SteadyStateDiffEqAlgorithm
    alg::A
    abstol::AT
    reltol::RT
    tspan::TS
end

function ContinuousDEQSolver(alg=VCABM4(); mode::Symbol=:rel_deq_default, abstol=1e-8, reltol=1e-8, tspan=Inf)
    return ContinuousDEQSolver{Val(mode),typeof(alg),typeof(abstol),typeof(reltol),typeof(tspan)}(alg, abstol, reltol, tspan)
end

function terminate_condition_reltol(integrator, abstol, reltol, min_t)
    return all(abs.(DiffEqBase.get_du(integrator)) .<= reltol .* abs.(integrator.u))
end

function terminate_condition_reltol_norm(integrator, abstol, reltol, min_t)
    du = DiffEqBase.get_du(integrator)
    return norm(du) <= reltol * norm(du .+ integrator.u)
end

function terminate_condition_abstol(integrator, abstol, reltol, min_t)
    return all(abs.(DiffEqBase.get_du(integrator)) .<= abstol)
end

function terminate_condition_abstol_norm(integrator, abstol, reltol, min_t)
    return norm(DiffEqBase.get_du(integrator)) <= abstol
end

function terminate_condition(integrator, abstol, reltol, min_t)
    return all((abs.(DiffEqBase.get_du(integrator)) .<= reltol .* abs.(integrator.u)) .&
               (abs.(DiffEqBase.get_du(integrator)) .<= abstol))
end

function terminate_condition_norm(integrator, abstol, reltol, min_t)
    du = DiffEqBase.get_du(integrator)
    du_norm = norm(du)
    return (du_norm <= reltol * norm(du .+ integrator.u)) && (du_norm <= abstol)
end

get_terminate_condition(::ContinuousDEQSolver{Val(:abs)}, args...; kwargs...) = terminate_condition_abstol
get_terminate_condition(::ContinuousDEQSolver{Val(:abs_norm)}, args...; kwargs...) = terminate_condition_abstol_norm
get_terminate_condition(::ContinuousDEQSolver{Val(:rel)}, args...; kwargs...) = terminate_condition_reltol
get_terminate_condition(::ContinuousDEQSolver{Val(:rel_norm)}, args...; kwargs...) = terminate_condition_reltol_norm
get_terminate_condition(::ContinuousDEQSolver{Val(:norm)}, args...; kwargs...) = terminate_condition_norm
get_terminate_condition(::ContinuousDEQSolver, args...; kwargs...) = terminate_condition

# Termination conditions used in the original DEQ Paper
function get_terminate_condition(::ContinuousDEQSolver{Val(:abs_deq_default),A,T}, args...; kwargs...) where {A,T}
    nstep = 0
    protective_threshold = T(1e6)
    objective_values = T[]
    function terminate_condition_closure(integrator, abstol, reltol, min_t)
        du = DiffEqBase.get_du(integrator)
        objective = norm(du)
        # Main termination condition
        objective <= abstol && return true

        # Terminate if there has been no improvement for the last 30 steps
        nstep += 1
        push!(objective_values, objective)

        objective <= 3 * abstol &&
            nstep >= 30 &&
            maximum(objective_values[(end - nstep):end]) < 1.3 * minimum(objective_values[(end - nstep):end]) &&
            return true

        # Protective break
        objective >= objective_values[1] * protective_threshold * length(du) && return true

        return false
    end
    return terminate_condition_closure
end

function get_terminate_condition(::ContinuousDEQSolver{Val(:rel_deq_default),A,T}, args...; kwargs...) where {A,T}
    nstep = 0
    protective_threshold = T(1e3)
    objective_values = T[]
    function terminate_condition_closure(integrator, abstol, reltol, min_t)
        du = DiffEqBase.get_du(integrator)
        u = integrator.u
        objective = norm(du) / (norm(du .+ u) + eps(T))
        # Main termination condition
        objective <= reltol && return true

        # Terminate if there has been no improvement for the last 30 steps
        nstep += 1
        push!(objective_values, objective)

        objective <= 3 * reltol &&
            nstep >= 30 &&
            maximum(objective_values[(end - nstep + 1):end]) < 1.3 * minimum(objective_values[(end - nstep + 1):end]) &&
            return true

        # Protective break
        objective >= objective_values[1] * protective_threshold * length(du) && return true

        return false
    end
    return terminate_condition_closure
end

function get_terminate_condition(::ContinuousDEQSolver{Val(:rel_deq_best),A,T}, terminate_stats::Dict, args...;
                                 kwargs...) where {A,T}
    nstep = 0
    protective_threshold = T(1e3)
    objective_values = T[]

    terminate_stats[:best_objective_value] = T(Inf)
    terminate_stats[:best_objective_value_iteration] = 0

    function terminate_condition_closure(integrator, abstol, reltol, min_t)
        du = DiffEqBase.get_du(integrator)
        u = integrator.u
        objective = norm(du) / (norm(du .+ u) + eps(T))

        if objective < terminate_stats[:best_objective_value]
            terminate_stats[:best_objective_value] = objective
            terminate_stats[:best_objective_value_iteration] = nstep + 1
        end

        # Main termination condition
        objective <= reltol && return true

        # Terminate if there has been no improvement for the last 30 steps
        nstep += 1
        push!(objective_values, objective)

        objective <= 3 * reltol &&
            nstep >= 30 &&
            maximum(objective_values[(end - nstep + 1):end]) < 1.3 * minimum(objective_values[(end - nstep + 1):end]) &&
            return true

        # Protective break
        objective >= objective_values[1] * protective_threshold * length(du) && return true

        return false
    end

    return terminate_condition_closure
end

function get_terminate_condition(::ContinuousDEQSolver{Val(:abs_deq_best),A,T}, terminate_stats::Dict, args...;
                                 kwargs...) where {A,T}
    nstep = 0
    protective_threshold = T(1e3)
    objective_values = T[]

    terminate_stats[:best_objective_value] = T(Inf)
    terminate_stats[:best_objective_value_iteration] = 0

    function terminate_condition_closure(integrator, abstol, reltol, min_t)
        du = DiffEqBase.get_du(integrator)
        objective = norm(du)

        if objective < terminate_stats[:best_objective_value]
            terminate_stats[:best_objective_value] = objective
            terminate_stats[:best_objective_value_iteration] = nstep + 1
        end

        # Main termination condition
        objective <= reltol && return true

        # Terminate if there has been no improvement for the last 30 steps
        nstep += 1
        push!(objective_values, objective)

        objective <= 3 * reltol &&
            nstep >= 30 &&
            maximum(objective_values[(end - nstep + 1):end]) < 1.3 * minimum(objective_values[(end - nstep + 1):end]) &&
            return true

        # Protective break
        objective >= objective_values[1] * protective_threshold * length(du) && return true

        return false
    end

    return terminate_condition_closure
end

has_converged(du, u, alg::ContinuousDEQSolver) = all(abs.(du) .<= alg.abstol .& abs.(du) .<= alg.reltol .* abs.(u))
has_converged(du, u, alg::ContinuousDEQSolver{Val(:norm)}) = norm(du) <= alg.abstol && norm(du) <= alg.reltol * norm(du .+ u)
has_converged(du, u, alg::ContinuousDEQSolver{Val(:rel)}) = all(abs.(du) .<= alg.reltol .* abs.(u))
has_converged(du, u, alg::ContinuousDEQSolver{Val(:rel_norm)}) = norm(du) <= alg.reltol * norm(du .+ u)
has_converged(du, u, alg::ContinuousDEQSolver{Val(:rel_deq_default)}) = norm(du) <= alg.reltol * norm(du .+ u)
has_converged(du, u, alg::ContinuousDEQSolver{Val(:rel_deq_best)}) = norm(du) <= alg.reltol * norm(du .+ u)
has_converged(du, u, alg::ContinuousDEQSolver{Val(:abs)}) = all(abs.(du) .<= alg.abstol)
has_converged(du, u, alg::ContinuousDEQSolver{Val(:abs_norm)}) = norm(du) <= alg.abstol
has_converged(du, u, alg::ContinuousDEQSolver{Val(:abs_deq_default)}) = norm(du) <= alg.abstol
has_converged(du, u, alg::ContinuousDEQSolver{Val(:abs_deq_best)}) = norm(du) <= alg.abstol

struct EquilibriumSolution{T,N,uType,R,P,A,TEnd} <: SciMLBase.AbstractNonlinearSolution{T,N}
    u::uType
    resid::R
    prob::P
    alg::A
    retcode::Symbol
    t::TEnd
    λₜ::T
end

function transform_solution(soln::EquilibriumSolution)
    # Creates a NonlinearSolution/SteadyStateSolution
    return DiffEqBase.build_solution(soln.prob, soln.alg, soln.u, soln.resid; retcode=soln.retcode)
end

function DiffEqBase.__solve(prob::DiffEqBase.AbstractSteadyStateProblem, alg::ContinuousDEQSolver, args...;
                            regularize_endpoint=false, kwargs...)
    tspan = alg.tspan isa Tuple ? alg.tspan : convert.(real(eltype(prob.u0)), (zero(alg.tspan), alg.tspan))
    _prob = ODEProblem(prob.f, prob.u0, tspan, prob.p)

    terminate_stats = Dict{Symbol,Any}(:best_objective_value => real(eltype(prob.u0))(Inf),
                                       :best_objective_value_iteration => nothing)

    sol = solve(_prob, alg.alg, args...; kwargs...,
                callback=TerminateSteadyState(alg.abstol, alg.reltol, get_terminate_condition(alg, terminate_stats)))

    u, t = terminate_stats[:best_objective_value_iteration] === nothing ? (sol.u[end], sol.t[end]) :
           (sol.u[terminate_stats[:best_objective_value_iteration] + 1],
            sol.t[terminate_stats[:best_objective_value_iteration] + 1])

    du = prob.f(u, prob.p, t)

    T = eltype(eltype(u))
    N = ndims(u)
    retcode = (sol.retcode == :Terminated && has_converged(du, u, alg) ? :Success : :Failure)
    _t = regularize_endpoint isa Bool ? (regularize_endpoint ? t : nothing) : t
    regularize_endpoint = regularize_endpoint isa Bool ? (regularize_endpoint ? T(1e-5) : T(0)) : T(regularize_endpoint)

    return EquilibriumSolution{T,N,typeof(u),typeof(du),typeof(prob),typeof(alg),typeof(_t)}(u, du, prob, alg, retcode,
                                                                                _t, regularize_endpoint)
end

function clear_zero(x::T) where T
    ϵ = eps(T)
    if -ϵ <= x < 0
        return -ϵ
    elseif 0 <= x < ϵ
        return ϵ
    end
    return x
end

@noinline function DiffEqSensitivity.SteadyStateAdjointProblem(sol::EquilibriumSolution,
                                                               sensealg::DiffEqSensitivity.SteadyStateAdjoint, g, dg;
                                                               save_idxs=nothing)
    @unpack f, p, u0 = sol.prob

    discrete = false

    p === DiffEqBase.NullParameters() &&
        error("Your model does not have parameters, and thus it is impossible to calculate the derivative of the solution with respect to the parameters. Your model must have parameters to use parameter sensitivity calculations!")

    sense = DiffEqSensitivity.SteadyStateAdjointSensitivityFunction(g, sensealg, discrete, sol, dg, f.colorvec, false)
    @unpack diffcache, y, sol, λ, vjp, linsolve = sense

    _save_idxs = save_idxs === nothing ? Colon() : save_idxs
    if dg !== nothing
        if g !== nothing
            dg(vec(diffcache.dg_val), y, p, nothing, nothing)
        else
            if typeof(_save_idxs) <: Number
                diffcache.dg_val[_save_idxs] = dg[_save_idxs]
            elseif typeof(dg) <: Number
                @. diffcache.dg_val[_save_idxs] = dg
            else
                @. diffcache.dg_val[_save_idxs] = dg[_save_idxs]
            end
        end
    else
        if g !== nothing
            DiffEqSensitivity.gradient!(vec(diffcache.dg_val), diffcache.g, y, sensealg, diffcache.g_grad_config)
        end
    end

    _val, back = Zygote.pullback(x -> f(x, p, nothing), y)
    s_val = size(_val)
    op = DiffEqSensitivity.ZygotePullbackMultiplyOperator{eltype(y),typeof(back),typeof(s_val)}(back, s_val)

    b = vec(diffcache.dg_val)
    # println("Original Mean: $(mean(b)) & Residual Mean: $(mean(sol.resid)) Norm: $(norm(sol.resid))")
    if sol.t !== nothing
        @. b = (b + clamp(sol.λₜ ./ norm(sol.resid), -Inf, mean(b))) / 2
    end
    # println("Updated mean: $(mean(b))")
    linear_problem = LinearProblem(op, b)

    copyto!(vec(λ), solve(linear_problem, linsolve).u)
    _, back = Zygote.pullback(p -> vec(f(y, p, nothing)), p)
    vjp .= -vec(back(λ)[1])

    if g !== nothing
        # compute del g/del p
        dg_dp_val = zero(p)
        dg_dp = DiffEqSensitivity.ParamGradientWrapper(g, nothing, y)
        dg_dp_config = DiffEqSensitivity.build_grad_config(sensealg, dg_dp, p, p)
        DiffEqSensitivity.gradient!(dg_dp_val, dg_dp, p, sensealg, dg_dp_config)

        @. dg_dp_val = dg_dp_val + vjp
        return dg_dp_val
    else
        return vjp
    end
end


text:
# This file is a part of RunStatistics.jl, licensed under the MIT License (MIT).#

using RunStatistics
using Test

@testset "partitions" begin

    p = RunStatistics.init_partition(6, 3)

    for i = 0:2

        sum = 0

        for j = 2:(p.h + 1)

            sum += p.c[j] * p.y[j]
        end

        @test sum == 6

        RunStatistics.next_partition!(p)
    end

    @test RunStatistics.next_partition!(p)

end

# NOTE:
#=

the number of partitions of an integer n into k parts satisfies the recurrence:
p_k(n) = p_k(n - k) + p_k - 1(n - 1)
with p_0(0) = 1 and p_k(n) = 0 uf n <= 0 or k <= 0 and n,k are not both 0

The set of partitions of 6 into 3 parts should therefore contain:

2*1 + 1*4
1*1 + 1*2 + 1*3
3*2

as the partitions

=#


text:
push!(LOAD_PATH, "src")
using DalitzPlotAnalysis
using Cuba
using Plots

const mπ=0.14; const mπ2 = mπ^2;

for v in [("I","-"), ("II","+")]
    @eval function $(Symbol("fρ_"*v[1]))(σ)
        mρ=0.7755
        Γρ=0.156
        return 1.0/(mρ^2-σ + $(Symbol(v[2]))(1.0im*mρ*Γρ))
    end
end
fρ_I(1.1)

function qtb(s)
  function integrand(x,f)
    s1min = 4.0 * mπ2
    s1max = (sqrt(s) - mπ)^2

    s1 =  s1min + x[1] * (s1max - s1min)
    f[1] = fρ_I(s1)*fρ_II(s1) * sqrt(λ(s,s1,mπ2)*λ(s1,mπ2,mπ2)) / (s*s1)*(s1max-s1min)
  end
  return cuhre(integrand, 1, 1)[1][1]
end

let ev = linspace(0.5,2.5,100)
    cal = [qtb(e^2) for e in ev]
    plot(ev, cal)
end

function angular_pow2(s, s1, z, m1sq, m2sq, m3sq)
    θ23 = acos(z)
    θ3 = acos(cross_basis_cosθ3(s1, z, m1sq, m2sq, m3sq, s))
    θ12 = acos(cross_basis_cosθ12(s1, z, m1sq, m2sq, m3sq, s))
    res = Wignerd(1,0,0,θ3+θ12-θ23);
end

function interference(s)
  function integrand(x,f)
    s1min = 4.0 * mπ2
    s1max = (sqrt(s) - mπ)^2

    s1 =  s1min + x[1] * (s1max - s1min)
    z = -1.0 + 2.0*x[2]
    s3 = cross_basis_s3(s1, z, mπ2, mπ2, mπ2, s)
    ampl2 = fρ_I(s1)*fρ_II(s3)*angular_pow2(s, s1, z, mπ2, mπ2, mπ2)
    f[1],f[2] = reim(ampl2 * sqrt(λ(s,s1,mπ2)*λ(s1,mπ2,mπ2)) / (s*s1)*(s1max-s1min))
  end
  return complex(cuhre(integrand, 2, 2)[1]...)
end

@time interference(1.1)

@profile interference(1.1)

let ev = linspace(0.5,2.5,100)
    cal = [interference(e^2) for e in ev]
    plot!(ev, real.(cal))
end

Juno.Profile.clear()
Juno.profiler()


text:
module DAQnidaqmx

using AbstractDAQs
using Dates
import NIDAQ

export NIDev, daqaddinput, daqstart, daqstop, daqread, daqacquire, daqconfig, daqconfigdev
export numchannels, daqchannels, isdaqfinished, isreading, samplesread

export NIException
export TermConfig, NIDefault, NIRSE, NINRSE, NIDiff, NIPseudoDiff

"Terminal configuration"
@enum TermConfig::Int32 NIDefault=NIDAQ.DAQmx_Val_Cfg_Default NIRSE=NIDAQ.DAQmx_Val_RSE NINRSE=NIDAQ.DAQmx_Val_NRSE NIDiff=NIDAQ.DAQmx_Val_Diff NIPseudoDiff=NIDAQ.DAQmx_Val_PseudoDiff

mutable struct NIDev <: AbstractDAQs.AbstractDAQ
    devname::String
    handle::NIDAQ.TaskHandle
    rate::Float64
    nsamples::Int64
    time::DateTime
    conf::DAQConfig
    chanidx::Dict{String,Int}
    NIDev(devname::String, handle::NIDAQ.TaskHandle) = new(devname, handle, -1.0, -1, now(),
                                                           DAQConfig(devname=devname),
                                                           Dict{String,Int}())
end

"""
`error_message(errnum, buflen)`

Return an error/warning message from its code.
"""
function error_message(errnum, buflen=1000)

    buf = zeros(UInt8, buflen)

    r = NIDAQ.DAQmxGetErrorString(errnum, buf, buflen)
    if r != 0
        return "$errnum: unknown error"
    end

    idx = findfirst(iszero, buf)
    idx1 =  isnothing(idx) ? buflen : idx-1
    return String(buf[1:idx1])
end

struct NIException <: Exception
    code::Int
    msg::String
end

"""
`NIException(code)`

Returns an `Exception` object specific for NIDAQmx. From a numeric code it
returns a description of the exception.
"""
function NIException(code)
    msg = error_message(code)
    return NIException(code, msg)
end

"""
`NIDev(devname, loadmaxtask=false)`

Creates a device that handles NIDAQmx data acquisition.
In the parlance used by NIDAQmx, it creates a task.

NI has tools such as NI MAX that can handle the creation and configuration of tasks.
With argument `loadmaxtask=true` these tasks can be used directly.

The device name (`devname`) is the task name to be used by NIDAQmx.

"""
function NIDev(devname::String, loadmaxtask::Bool=false)

    th = Ref(NIDAQ.TaskHandle(0))
    if !loadmaxtask
        r = NIDAQ.DAQmxCreateTask(devname, th)
        r != 0 && throw(NIException(r))
        handle = th[]
    else
        r = NIDAQ.DAQmxLoadTask(devname, th)
        r != 0 && throw(NIException(r))
        handle = th[]
    end

    return NIDev(devname, handle)
end

"""
`stoptask(dev)`

Stops the execution of a task. It calls NIDAQmx function `DAQmxStopTask`.
"""
function stoptask(dev::NIDev)

    r = NIDAQ.DAQmxStopTask(dev.handle)

    r != 0 && throw(NIException(r))

end

"""
`cleartask(dev)`

Clears a task. It calls NIDAQmx function `DAQmxClearTask`. After calling this function,
the task must be reconfigured.
"""
function cleartask(dev::NIDev)

    r = NIDAQ.DAQmxClearTask(dev.handle)

    r != 0 && throw(NIException(r))

end

"""
`daqstop(dev)`

Just a wrapper around [`stoptask`](@ref) implementing the [`AbstractDAQs` ](@ref) interface.
"""
function AbstractDAQs.daqstop(dev::NIDev)
    stoptask(dev)
end

"""
`setchanidx!(dev::NIDev)`

Reads the channel names and sets up `chanidx` Dict field that
maps from channel names (`String`) to index (Integer)

"""
function setchanidx!(dev::NIDev)

    chn = daqchannels(dev)
    chd = Dict{String,Int}()
    for (i,c) in enumerate(chn)
        chd[c] = i
    end
    dev.chanidx = chd
    return

end

"""
`daqaddinput(dev::NIDev, chans, names; termconf, minval, maxmval, units, customscale)`

Adds input to `NIDev` object. For now it handles only AI Voltage channels. In the future,
more specific input types will be implemented.

## Arguments
 * `chans`: Physical channels usually something like "Dev1/ai0". See NIDAQmx documentation.
 * `names`: Channel names, if "" the physical names are used.
 * `termconf`: Terminal configuration (RSE, NRSE, Differential, Pseudodifferential). See [`TermConf`] and NIDAQmx docs for further information
 * `minval` and `maxval`:

"""
function AbstractDAQs.daqaddinput(dev::NIDev, chans::AbstractString; names::AbstractString="",
                     termconf=NIDefault, minval=0.0, maxval=5.0,
                     units=NIDAQ.DAQmx_Val_Volts, customscalename="")
    r = NIDAQ.DAQmxCreateAIVoltageChan(dev.handle, chans, names,
                                       termconf, minval, maxval,
                                       units, customscalename)
    r != 0 && throw(NIException(r))

    dev.conf.fpars["minval"] = minval
    dev.conf.fpars["minval"] = maxval
    dev.conf.ipars["units"] = units
    dev.conf.ipars["termconf"] = Int(termconf)
    dev.conf.spars["customscalename"] = customscalename
    dev.conf.spars["nichans"] = chans

    setchanidx!(dev)

    return
end

function AbstractDAQs.daqaddinput(dev::NIDev, devname::AbstractString,
                     chans::AbstractVector{<:Integer}; names="",
                     termconf=NIDiff, s = "ai",
                     minval=0.0, maxval=5.0,
                     units=NIDAQ.DAQmx_Val_Volts, customscalename="")

    # Create a string appropriate for

    if isa(chans, UnitRange)
        ch = "$devname/$s$(first(chans)):$(last(chans))"
    else
        ch = join(("$devname/$s$i" for i in chans), ",")
    end
    if names==""
        chnames = ""
    elseif isa(names, AbstractString)
        chnames = join((names[1]*string(c) for c in chans), ",")
    elseif isa(names, AbstractVector)
        #if length(names) == length(chans)
        chnames = join(names, ",")
    else
        throw("names: $names. Illegal format!")
    end
    daqaddinput(dev, ch; names=chnames, termconf=termconf, minval=minval, maxval=maxval,
                units=units, customscalename=customscalename)

    setchanidx!(dev)

    return

end

function AbstractDAQs.numchannels(dev::NIDev)
    num = zeros(UInt32, 1)
    r = NIDAQ.DAQmxGetTaskNumChans(dev.handle, num)

    r != 0 && throw(NIException(r))

    return Int(num[1])

end

function AbstractDAQs.daqchannels(dev::NIDev)
    nch = numchannels(dev)

    nperch = 30 # Number of characters per channel.
    buflen = nperch * nch
    buf = zeros(UInt8, buflen)
    r = NIDAQ.DAQmxGetTaskChannels(dev.handle, buf, buflen)

    r != 0 && throw(NIException(r))

    idx = findfirst(iszero, buf)
    idx1 =  isnothing(idx) ? buflen : idx-1
    s = String(buf[1:idx1])
    return [strip(ch) for ch in split(s, ',')]
end

function AbstractDAQs.daqconfigdev(dev::NIDev; source="",
                       samplemode=NIDAQ.DAQmx_Val_FiniteSamps,
                       activeedge=NIDAQ.DAQmx_Val_Rising, kw...)

    if haskey(kw, :rate) && haskey(kw, :dt)
        error("Parameters `rate` and `dt` can not be specified simultaneously!")
    elseif haskey(kw, :rate) || haskey(kw, :dt)
        if haskey(kw, :rate)
            rate = kw[:rate]
        else
            dt = kw[:dt]
            rate = 1/dt
        end
    else
        rate = dev.conf.fpars["rate"]
    end

    if haskey(kw, :nsamples) && haskey(kw, :time)
        error("Parameters `nsamples` and `time` can not be specified simultaneously!")
    elseif haskey(kw, :nsamples) || haskey(kw, :time)
        if haskey(kw, :nsamples)
            nsamples = UInt64(kw[:nsamples])
        else
            tt = kw[:time]
            dt = 1/rate
            nsamples = round(UInt64, tt / dt)
        end
    else
        nsamples = UInt64(dev.conf.ipars("nsamples"))
    end

    r = NIDAQ.DAQmxCfgSampClkTiming(dev.handle, source, rate,
                                    activeedge, samplemode, nsamples)
    r != 0 && throw(NIException(r))

    dev.rate = rate
    dev.nsamples = nsamples
    dev.conf.fpars["rate"] = rate
    dev.conf.ipars["nsamples"] = nsamples
    dev.conf.ipars["samplemode"] = samplemode
    dev.conf.ipars["activeedge"] = activeedge

    return

end

function AbstractDAQs.daqconfig(dev::NIDev; source="",
                                samplemode=NIDAQ.DAQmx_Val_FiniteSamps,
                                activeedge=NIDAQ.DAQmx_Val_Rising, kw...)
    daqconfigdev(dev; source=source, samplemode=samplemode, activeedge=activeedge, kw...)
end

function AbstractDAQs.isdaqfinished(dev::NIDev)

    data = zeros(Int32,1)
    #data = reinterpret(NIDAQ.Bool32, zeros(Int32, 1))

    r = NIDAQ.DAQmxIsTaskDone(dev.handle, data)
    r != 0 && throw(NIException(r))
    return data[1] != 0
end

AbstractDAQs.isreading(dev::NIDev) = !isdaqfinished(dev)

function AbstractDAQs.daqstart(dev::NIDev, usethread=true)
    r = NIDAQ.DAQmxStartTask(dev.handle)
    r != 0 && throw(NIException(r))
    dev.time = now()
    return
end

function AbstractDAQs.daqread(dev::NIDev)

    nch = numchannels(dev)
    nsamples = dev.nsamples
    buffer = zeros(Float64, nch, nsamples)
    fm1 =  Int32[NIDAQ.DAQmx_Val_GroupByScanNumber]
    fm =  reinterpret(NIDAQ.Bool32, fm1)[1]
    nsr = Int32[1]

    # Wait for reading to end:
    while !isdaqfinished(dev)
        sleep(0.05)
    end

    r = NIDAQ.DAQmxReadAnalogF64(dev.handle, -1, -1.0, fm, buffer, nsamples*nch, nsr,
                                 Ptr{NIDAQ.Bool32}(0))
    r != 0 && throw(NIException(r))

    stoptask(dev)

    return MeasData{Matrix{Float64},Int}(devname(dev), devtype(dev), dev.time,
                                         dev.rate, buffer, 0, dev.chanidx)

end

function AbstractDAQs.daqacquire(dev::NIDev)
    daqstart(dev)
    return daqread(dev)
end

function AbstractDAQs.samplesread(dev::NIDev)

    data = zeros(UInt64,1)
    r = NIDAQ.DAQmxGetReadTotalSampPerChanAcquired(dev.handle, data)
    r != 0 && throw(NIException(r))

    return Int64(data[1])
end

end # module DAQnidaqmx


text:
# This file was generated, do not modify it. # hide
tuning = Grid(resolution=3)
resampling = CV(nfolds=6)

tm = TunedModel(model=krb, tuning=tuning, resampling=resampling,
                ranges=ranges, measure=rmsl)

max_stars_repo_name: tlienart/MLJTutorials

max_issues_repo_name: tlienart/MLJTutorials

max_forks_repo_name: tlienart/MLJTutorials


text:
"""
    pythagorean_triplets(n)

Find all positive integer triplets `(a, b, c)` s.t. `a + b + c = n` and `a < b < c` and `a^2 + b^2 == c^2`.
"""
# cmcaine's answer, with thanks to akshu3398.
function pythagorean_triplets(n)
    triplets = NTuple{3, Int}[]
    # Lower bound because the smallest triple is 3, 4, 5.
    # Upper bound implied by a < b < c && a + b + c == n.
    for a in 3:cld(n, 3) - 1
        # Derived by eliminating c from these simultaneous
        # equations and solving for b:
        #     a^2 + b^2 = c^2
        #     a + b + c = n
        b = (n^2 - 2n * a) / (2 * (n - a))
        if a < b && isinteger(b)
            c = n - a - b
            # Proof that b < c:
            # Let b = c. Then a^2 + b^2 = c^2 ≡ a^2 = 0 but we know that a ≠ 0.
            # Let b > c. a^2 + b^2 = c^2 ≡ a^2 = c^2 - b^2.
            #    If b > c, a^2 < 0, but that is impossible for real numbers.
            push!(triplets, (a, b, c))
        end
    end
    return sort!(triplets)
end


text:
my_funct(x, y) = 2*x + y
der(x,y) = ForwardDiff.derivative(x->my_funct(x,y), x)


text:
# set initial values of rhoVec
function set_rho_vec!(ws::COSMO.Workspace)
	m = ws.p.model_size[1]
	# nEQ = p.K.f
	# nINEQ = p.m - p.K.f
	ws.ρ = ws.settings.rho
	ws.ρvec = ws.ρ * ones(m)

	# scale ρ values that belong to equality constraints with a factor of 1e3
	set_ind = findall(x -> typeof(x) == COSMO.ZeroSet{Float64}, ws.p.C.sets)
	if length(set_ind) > 0
		row_ind = COSMO.get_set_indices(ws.p.C.sets)
		for (i, rows) in enumerate(row_ind[set_ind])
			ws.ρvec[rows] *= 1e3
		end
	end
	push!(ws.Info.rho_updates, ws.ρ)
	return nothing
end

# adapt rhoVec based on residual ratio
function adapt_rho_vec!(ws::COSMO.Workspace)
	settings = ws.settings
	# compute normalized residuals based on the working variables (dont unscale)
	ignore_scaling = true
	r_prim::Float64, r_dual::Float64 = calculate_residuals!(ws, ignore_scaling)
	max_norm_prim::Float64, max_norm_dual::Float64  = max_res_component_norm(ws, ignore_scaling)

	r_prim = r_prim / (max_norm_prim + 1e-10)
	r_dual = r_dual / (max_norm_dual + 1e-10)

	new_rho = ws.ρ * sqrt(r_prim / (r_dual + 1e-10))
	new_rho = min(max(new_rho, settings.RHO_MIN), settings.RHO_MAX)
	# only update rho if significantly different than current rho
	if (new_rho > settings.adaptive_rho_tolerance * ws.ρ) || (new_rho < (1 ./ settings.adaptive_rho_tolerance) * ws.ρ)
		update_rho_vec!(new_rho, ws)
	end
	return nothing
end

function update_rho_vec!(new_rho::Float64, ws::COSMO.Workspace)

	ws.ρ     = new_rho
	ws.ρvec .= new_rho

	# scale ρ values that belong to equality constraints with a factor of 1e3
	set_ind = findall(x -> typeof(x) == COSMO.ZeroSet{Float64}, ws.p.C.sets)
	if length(set_ind) > 0
		row_ind = COSMO.get_set_indices(ws.p.C.sets)
		for (i, rows) in enumerate(row_ind[set_ind])
			ws.ρvec[rows] *= 1e3
		end
	end

	# log rho updates to info variable
	push!(ws.Info.rho_updates, new_rho)

	if ws.settings.verbose_timing
		ws.times.factor_time += @elapsed update_rho!(ws.kkt_solver,ws.ρvec)
	else
		update_rho!(ws.kkt_solver,ws.ρvec)
	end

	return nothing
end