text:
module LazyApply

using Compat
using Zeros

import Base:
  *,
  ^,
  +,
  -,
  /,
  exp,
  adjoint,
  reverse,
  show,
  ==,
  convert,
  getindex,
  length,
  size,
  iterate,
  lastindex

export coefficient, expand, Sum, Prod, coefficient

struct Applied{F,Args}
  f::F
  args::Args
  function Applied{F,Args}(f, args::Tuple) where {F,Args}
    return new{F,Args}(f, args)
  end
end
Applied(f, args::Tuple) = Applied{typeof(f),typeof(args)}(f, args)
Applied(f, args...) = Applied(f, args)

# TODO: This makes shorthands like `Add(1, 2)` work, but probably
# it is bad to use `F.instance` to get the function from the type.
Applied{F,Args}(args::Tuple) where {F,Args} = Applied{F,Args}(F.instance, args)
Applied{F,Args}(args...) where {F,Args} = Applied{F,Args}(args)
Applied{F}(args::Tuple) where {F} = Applied{F,typeof(args)}(args)
Applied{F}(args...) where {F} = Applied{F}(args)

# if VERSION < v"1.6"
const AppliedTupleAB{F} = Applied{F,Tuple{A,B}} where {A,B}
const AppliedTupleA{F,B} = Applied{F,Tuple{A,B}} where {A}
const AppliedTupleVector{F} = Applied{F,Tuple{Vector{T}}} where {T}

# For `Scaled(3.2, "X")`
## if VERSION > v"1.5"
##   (Applied{F,Tuple{A,B}} where {A,B})(args::Tuple) where {F} = Applied{F}(args)
##   (Applied{F,Tuple{A,B}} where {A,B})(args...) where {F} = Applied{F}(args)
## else
AppliedTupleAB{F}(args::Tuple) where {F} = Applied{F}(args)
AppliedTupleAB{F}(args...) where {F} = Applied{F}(args)

# For `Scaled{ComplexF64}(3.2, "X")`
## if VERSION > v"1.5"
##   function (Applied{F,Tuple{A,B}} where {A})(args::Tuple{Arg1,Arg2}) where {F,B,Arg1,Arg2}
##     return Applied{F,Tuple{Arg1,B}}(args)
##   end
##   function (Applied{F,Tuple{A,B}} where {A})(args...) where {F,B}
##     return (Applied{F,Tuple{A,B}} where {A})(args)
##   end
## end
function AppliedTupleA{F,B}(args::Tuple{Arg1,Arg2}) where {F,B,Arg1,Arg2}
  return Applied{F,Tuple{Arg1,B}}(args)
end
function AppliedTupleA{F,B}(args...) where {F,B}
  return (Applied{F,Tuple{A,B}} where {A})(args)
end

# For `Sum([1, 2, 3])` and `Prod([1, 2, 3])`
## (Applied{F,Tuple{Vector{T}}} where {T})(args::Vector) where {F} = Applied{F}((args,))
AppliedTupleVector{F}(args::Vector) where {F} = Applied{F}((args,))

_empty(::Type{T}) where {T} = error("_empty not implemented for type $T.")
_empty(::Type{Tuple{T}}) where {T} = (_empty(T),)
_empty(::Type{Vector{T}}) where {T} = Vector{T}()

function Applied{F,Args}() where {F,Args}
  return Applied{F,Args}(_empty(Args))
end

_initialvalue_type(::Type{typeof(+)}) = Zero
_initialvalue_type(::Type{typeof(sum)}) = Zero
_initialvalue_type(::Type{typeof(*)}) = One
_initialvalue_type(::Type{typeof(prod)}) = One

# For `Sum() == Sum{Zero}()` and `Prod() == Prod{One}()`
## function (Applied{F,Tuple{Vector{T}}} where {T})() where {F}
##   return Applied{F,Tuple{Vector{_initialvalue_type(F)}}}()
## end
function AppliedTupleVector{F}() where {F}
  return Applied{F,Tuple{Vector{_initialvalue_type(F)}}}()
end

function (arg1::Applied == arg2::Applied)
  return (arg1.f == arg2.f) && (arg1.args == arg2.args)
end

# Shorthands
const Add{T} = Applied{typeof(+),T}
const Mul{T} = Applied{typeof(*),T}
# By default, `A` is the scalar, but the type constraint isn't
# working for some reason.
# A scaled operator scaled by a `Float64` of type `Scaled{Op,Float64}`.
# A scaled operator with an unspecified scalar type is of type `Scaled{Op}`.
const Scaled{B,A} = Applied{typeof(*),Tuple{A,B}}
const Sum{T} = Applied{typeof(sum),Tuple{Vector{T}}}
const Prod{T} = Applied{typeof(prod),Tuple{Vector{T}}}
const ∑ = Sum
const ∏ = Prod
const α = Scaled

const Exp{T} = Applied{typeof(exp),Tuple{T}}

coefficient(arg::Applied) = 𝟏

length(arg::Union{Sum,Prod}) = length(only(arg.args))
size(arg::Union{Sum,Prod}, args...) = size(only(arg.args), args...)
lastindex(arg::Union{Sum,Prod}) = length(arg)
getindex(arg::Union{Sum,Prod}, n) = getindex(arg.args..., n)
iterate(arg::Union{Sum,Prod}, args...) = iterate(arg.args..., args...)

length(arg::Scaled) = length(arg.args[2])
getindex(arg::Scaled, n) = getindex(arg.args[2], n)
coefficient(arg::Scaled) = arg.args[1]
iterate(arg::Scaled, args...) = iterate(arg.args[2], args...)

Base.convert(::Type{Applied{F,Args}}, arg::Applied{F,Args}) where {F,Args} = arg

# For some reasons this conversion isn't being done automatically.
function convert(::Type{Applied{F,Args1}}, arg2::Applied{F,Args2}) where {F,Args1,Args2}
  return Applied{F,Args1}(arg2.f, convert(Args1, arg2.args))
end

# Just like `Base.promote`, but this doesn't error if
# a conversion doesn't happen.
function try_promote(x::T, y::S) where {T,S}
  R = promote_type(T, S)
  return (convert(R, x), convert(R, y))
end

# Conversion
Sum(arg::Add) = Sum(collect(arg.args))
Sum(arg::Sum) = arg

# Scalar multiplication (general rules)
_mul(arg1::Number, arg2) = Mul(arg1, arg2)
(arg1::Number * arg2::Applied) = _mul(arg1, arg2)
(arg1::Number * arg2::Prod) = _mul(arg1, arg2)

# Scalar division
(arg1::Applied / arg2::Number) = inv(arg2) * arg1

# Put the scalar value first by convention
_mul(arg1, arg2::Number) = Mul(arg2, arg1)
(arg1::Applied * arg2::Number) = _mul(arg1, arg2)
(arg1::Prod * arg2::Number) = _mul(arg1, arg2)

# Scalar multiplication (specialized rules)
(arg1::Number * arg2::Scaled) = Mul(arg1 * arg2.args[1], arg2.args[2])
(arg1::Scaled * arg2::Scaled) =
  Mul(arg1.args[1] * arg2.args[1], arg1.args[2] * arg2.args[2])
(arg1::Scaled * arg2) = Mul(arg1.args[1], arg1.args[2] * arg2)
(arg1 * arg2::Scaled) = Mul(arg2.args[1], arg1 * arg2.args[2])
# Scalars are treated special for the sake of multiplication
Mul(arg1::Number, arg2::Number) = arg1 * arg2
Mul(arg1::Number, arg2::Scaled) = arg1 * arg2
Mul(arg1::Scaled, arg2::Number) = arg1 * arg2
(arg1::Number * arg2::Sum) = Sum(map(a -> Mul(arg1, a), arg2))
(arg1::Number * arg2::Add) = Add(map(a -> Mul(arg1, a), arg2.args))

# Types should implement `__sum`.
_sum(arg1, arg2) = __sum(try_promote(arg1, arg2)...)

# Addition (general rules)
__sum(arg1, arg2) = Sum(vcat(arg1, arg2))
(arg1::Applied + arg2::Applied) = _sum(arg1, arg2)
(arg1::Applied + arg2) = _sum(arg1, arg2)
(arg1 + arg2::Applied) = _sum(arg1, arg2)

# Subtraction (general rules)
_subtract(arg1, arg2) = _sum(arg1, Mul(-𝟏, arg2))
(arg1::Applied - arg2::Applied) = _subtract(arg1, arg2)
(arg1::Applied - arg2) = _subtract(arg1, arg2)
(arg1 - arg2::Applied) = _subtract(arg1, arg2)

# Addition (specialized rules)
__sum(arg1::Sum, arg2::Sum) = Sum(vcat(arg1.args..., arg2.args...))
(arg1::Sum + arg2::Sum) = _sum(arg1, arg2)
(arg1::Add + arg2::Add) = Add(arg1.args..., arg2.args...)

__sum(arg1::Sum, arg2) = Sum(vcat(arg1.args..., arg2))
__sum(arg1, arg2::Sum) = Sum(vcat(arg1, arg2.args...))
(arg1::Sum + arg2) = _sum(arg1, arg2)
(arg1::Sum + arg2::Applied) = _sum(arg1, arg2)

(arg1 + arg2::Sum) = _sum(arg1, arg2)
(arg1::Add + arg2) = Add(arg1.args..., arg2)
(arg1 + arg2::Add) = Add(arg1, arg2.args...)

# Multiplication (general rules)
(arg1::Applied * arg2::Applied) = Prod([arg1, arg2])

# Multiplication (specialized rules)
(arg1::Prod * arg2::Prod) = Prod(vcat(arg1.args..., arg2.args...))
(arg1::Sum * arg2::Sum) = Prod([arg1, arg2])

_prod(arg1::Prod{One}, arg2) = Prod(vcat(arg2))
_prod(arg1::Prod{One}, arg2::Vector) = Prod(arg2)
_prod(arg1::Prod, arg2) = Prod(vcat(arg1.args..., arg2))
(arg1::Prod * arg2) = _prod(arg1, arg2)
(arg1::Prod * arg2::Applied) = _prod(arg1, arg2)

_prod(arg1, arg2::Prod) = Prod(vcat(arg1, arg2...))
(arg1 * arg2::Prod) = _prod(arg1, arg2)
(arg1::Applied * arg2::Prod) = _prod(arg1, arg2)

# Generically make products
(arg1::Applied * arg2) = Prod(vcat(arg1, arg2))
(arg1 * arg2::Applied) = Prod(vcat(arg1, arg2))

function (arg1::Applied^arg2::Integer)
  res = ∏()
  for n in 1:arg2
    res *= arg1
  end
  return res
end

# Other lazy operations
exp(arg::Applied) = Applied(exp, arg)

# adjoint
adjoint(arg::Applied) = Applied(adjoint, arg)
adjoint(arg::Applied{typeof(adjoint)}) = only(arg.args)
adjoint(arg::Prod) = ∏(reverse(adjoint.(arg)))

# reverse
reverse(arg::Prod) = Prod(reverse(arg.args...))

# Materialize
materialize(a::Number) = a
materialize(a::AbstractString) = a
materialize(a::Vector) = materialize.(a)
materialize(a::Applied) = a.f(materialize.(a.args)...)

function _expand(a1::Sum, a2::Sum)
  return ∑(vec([a1[i] * a2[j] for i in 1:length(a1), j in 1:length(a2)]))
end

# Expression manipulation
function expand(a::Prod{<:Sum})
  if length(a) == 1
    return a[1]
  elseif length(a) ≥ 2
    a12 = _expand(a[1], a[2])
    return expand(∏(vcat(a12, a[3:end])))
  end
end

_print(io::IO, args...) = print(io, args...)
function _print(io::IO, a::AbstractVector, args...)
  print(io, "[")
  for n in 1:length(a)
    _print(io, a[n], args...)
    if n < length(a)
      print(io, ",\n")
    end
  end
  return print(io, "]")
end

function show(io::IO, m::MIME"text/plain", a::Applied)
  print(io, a.f, "(\n")
  for n in 1:length(a.args)
    _print(io, a.args[n])
    if n < length(a.args)
      print(io, ", ")
    end
  end
  return print(io, "\n)")
end
show(io::IO, a::Applied) = show(io, MIME("text/plain"), a)
end


text:
using Gridap
using ForwardDiff
using LinearAlgebra
using Test
using GridapODEs.ODETools
using GridapODEs.TransientFETools
using Gridap.FESpaces: get_algebraic_operator

import Gridap: ∇
import GridapODEs.TransientFETools: ∂t
∂
θ = 0.2

u(x,t) = (1.0-x[1])*x[1]*(1.0-x[2])*x[2]*t
u(t::Real) = x -> u(x,t)
f(t) = x -> ∂t(u)(x,t)-Δ(u(t))(x)

domain = (0,1,0,1)
partition = (2,2)
model = CartesianDiscreteModel(domain,partition)

order = 2

V0 = FESpace(
  reffe=:Lagrangian, order=order, valuetype=Float64,
  conformity=:H1, model=model, dirichlet_tags="boundary")
U = TransientTrialFESpace(V0,u)

trian = Triangulation(model)
degree = 2*order
quad = CellQuadrature(trian,degree)

#
a(u,v) = ∇(v)⋅∇(u)
b(v,t) = v*f(t)

res(t,u,ut,v) = a(u,v) + ut*v - b(v,t)
jac(t,u,ut,du,v) = a(du,v)
jac_t(t,u,ut,dut,v) = dut*v

t_Ω = FETerm(res,jac,jac_t,trian,quad)
op = TransientFEOperator(U,V0,t_Ω)

t0 = 0.0
tF = 1.0
dt = 0.1

U0 = U(0.0)
uh0 = interpolate_everywhere(U0,u(0.0))

ls = LUSolver()
using Gridap.Algebra: NewtonRaphsonSolver
nls = NLSolver(ls;show_trace=true,method=:newton) #linesearch=BackTracking())
odes = ThetaMethod(ls,dt,θ)
solver = TransientFESolver(odes)

sol_t = solve(solver,op,uh0,t0,tF)

# Juno.@enter Base.iterate(sol_t)

l2(w) = w*w

tol = 1.0e-6
_t_n = t0

for (uh_tn, tn) in sol_t
  global _t_n
  _t_n += dt
  e = u(tn) - uh_tn

  el2 = sqrt(sum( integrate(l2(e),trian,quad) ))

  @test el2 < tol
end

ext: jl

max_stars_repo_name: stevengogogo/ProgrammableAging

max_issues_repo_name: stevengogogo/ProgrammableAging

max_forks_repo_name: stevengogogo/ProgrammableAging


text:
mutable struct ContactMethods
	J
	dJ
	ϕ
	d
	dλ1
	cf
	dcf
	vs
	vsq2
	vsq1h
	mdvs
	mdψη
	rc
	rcz
	rcθ
end

function ContactMethods()
	function f()
		error("Not Implemented: use instantiate_contact_methods!")
		return nothing
	end
	return ContactMethods(fill(f, 15)...)
end

function contact_forces(model::ContactModel, env::Environment{<:World,LinearizedCone}, γ1, b1, q2, k)
	nc = model.dim.c
	nb = nc * friction_dim(env)
	nf = Int(nb / nc)
	ne = dim(env)
	λ1 = vcat([transpose(rotation(env, k[(i-1) * (ne - 1) .+ (1:ne)])) * [friction_mapping(env) * b1[(i-1) * nf .+ (1:nf)]; γ1[i]] for i = 1:nc]...) # TODO: make efficient
end

function contact_forces(model::ContactModel, env::Environment{<:World,NonlinearCone}, γ1, b1, q2, k)
	nc = model.dim.c
	nb = nc * friction_dim(env)
	ne = dim(env)
	λ1 = vcat([transpose(rotation(env, k[(i-1) * (ne - 1) .+ (1:ne)])) * [b1[(i-1) * (ne-1) .+ (1:(ne-1))]; γ1[i]] for i = 1:nc]...) # TODO: make efficient
end

function velocity_stack(model::ContactModel, env::Environment{<:World,LinearizedCone}, q1, q2, k, h)
	nc = model.dim.c
	ne = dim(env)
	v = J_func(model, env, q2) * (q2 - q1) / h[1]
	v_surf = [rotation(env, k[(i-1) * (ne - 1) .+ (1:ne)]) * v[(i-1) * ne .+ (1:ne)] for i = 1:nc]
	vT_stack = vcat([[v_surf[i][1:ne-1]; -v_surf[i][1:ne-1]] for i = 1:nc]...)
end

function velocity_stack(model::ContactModel, env::Environment{<:World,NonlinearCone}, q1, q2, k, h)
	nc = model.dim.c
	ne = dim(env)
	v = J_func(model, env, q2) * (q2 - q1) / h[1]
	v_surf = [rotation(env, k[(i-1) * ne .+ (1:ne)]) * v[(i-1) * ne .+ (1:ne)] for i = 1:nc]
	vT_stack = vcat([v_surf[i][1:ne-1] for i = 1:nc]...)
end


text:
"""
    import_bibtex(input)
Import a BibTeX file or parse a BibTeX string and convert it to the internal bibliography format.
"""
function import_bibtex(input)
    return isfile(input) ? BibParser.parse_file(input) : BibParser.parse_entry(input)
end

"""
    int_to_spaces(n)

Make a string of `n` spaces.
"""
int_to_spaces(n) = repeat(" ", n)

const spaces = Dict{String,String}(map(
    s -> (string(s) => int_to_spaces(BibInternal.space(s))),
    BibInternal.fields)
)

"""
    field_to_bibtex(key, value)

Convert an entry field to BibTeX format.
"""
function field_to_bibtex(key, value)
    space = get(spaces, key, int_to_spaces(BibInternal.space(Symbol(key))))
    swp = length(key) > 3 && key[1:3] == "swp"
    o,f = isnothing(match(r"@", value)) ? ('{','}') : ('"','"')
    return value == "" || swp ? "" : " $key$space = $o$value$f,\n"
end

"""
    name_to_string(name::BibInternal.Name)

Convert a name in an `Entry` to a string.
"""
function name_to_string(name)
    str = "$(name.particle)"
    if str != "" != name.last
        str *= " "
    end
    str *= name.last
    str *= name.junior == "" ? "" : ", $(name.junior)"
    if name.first != ""
        str *= ", $(name.first)"
    end
    if name.middle != ""
        str *= " $(name.middle)"
    end
    return str
end

"""
    names_to_strings(names)

Convert a collection of names to a BibTeX string.
"""
function names_to_strings(names)
    return join(map(name_to_string, names), " and ")
end

"""
    access_to_bibtex!(fields, a)

Transform the how-to-`access` field to a BibTeX string.
"""
function access_to_bibtex!(fields, a)
    fields["doi"] = a.doi
    fields["howpublished"] = a.howpublished
    fields["url"] = a.url
end

"""
    date_to_bibtex!(fields, date)

Convert a date to a BibTeX string.
"""
function date_to_bibtex!(fields, date)
    fields["day"] = date.day
    fields["month"] = date.month
    fields["year"] = date.year
end

"""
    eprint_to_bibtex!(fields, eprint)

Convert eprint information to a BibTeX string.
"""
function eprint_to_bibtex!(fields, ep)
    fields["archivePrefix"] = ep.archive_prefix
    fields["eprint"] = ep.eprint
    fields["primaryClass"] = ep.primary_class
end

"""
    in_to_bibtex!(fields::BibInternal.Fields, i::BibInternal.In)

Convert the "published `in`" information to a BibTeX string.
"""
function in_to_bibtex!(fields, in_)
    fields["address"] = in_.address
    fields["chapter"] = in_.chapter
    fields["edition"] = in_.edition
    fields["institution"] = in_.institution
    fields["journal"] = in_.journal
    fields["number"] = in_.number
    fields["organization"] = in_.organization
    fields["pages"] = in_.pages
    fields["publisher"] = in_.publisher
    fields["school"] = in_.school
    fields["series"] = in_.series
    fields["volume"] = in_.volume
end

"""
    export_bibtex(e::Entry)

Export an `Entry` to a BibTeX string.
"""
function export_bibtex(e::Entry)
    access_to_bibtex!(e.fields, e.access)
    e.fields["author"] = names_to_strings(e.authors)
    e.fields["booktitle"] = e.booktitle
    date_to_bibtex!(e.fields, e.date)
    e.fields["editor"] = names_to_strings(e.editors)
    eprint_to_bibtex!(e.fields, e.eprint)
    in_to_bibtex!(e.fields, e.in)
    e.fields["title"] = e.title

    str = "@$(e.type == "eprint" ? "misc" : e.type){$(e.id),\n"
    for (name, value) in collect(e.fields)
        m = match(r"swp-",name)
        if m === nothing || m.offset > 1
            str *= value == "" ? "" : field_to_bibtex(name, value)
        end
    end
    return str[1:end - 2] * "\n}"
end

"""
    export_bibtex(bibliography)

Export a bibliography to a BibTeX string.
"""
function export_bibtex(bibliography)
    str = ""
    for e in values(bibliography)
        str *= export_bibtex(e) * "\n\n"
    end
    return str[1:end - 1]
end


text:
# Elementary vector spaces
#--------------------------
abstract ElementarySpace{F} <: VectorSpace
typealias IndexSpace ElementarySpace # index spaces for tensors
# Elementary finite-dimensional vector space over a field F that can be used as the index space corresponding to the indices
# of a tensor. Every elementary vector space should respond to the methods conj and dual, returning the complex conjugate space
# and the dual space respectively. The complex conjugate of the dual space is obtained as dual(conj(V)) === conj(dual(V)). These
# different spaces should be of the same type, so that a tensor can be defined as an element of a homogeneous tensor product
# of these spaces

Base.eltype{F}(V::ElementarySpace{F}) = F
Base.eltype{F}(::Type{ElementarySpace{F}}) = F
Base.eltype{S<:ElementarySpace}(::Type{S}) = eltype(super(S))

Base.isfinite(::ElementarySpace) = true

# The complex conjugate vector space of a real vector space is equal to itself
Base.conj{F<:Real}(V::ElementarySpace{F}) = V

abstract InnerProductSpace{F} <: ElementarySpace{F}
abstract HilbertSpace{F} <: InnerProductSpace{F}
# An inner product space, has the possibility to raise or lower indices of tensors using the metric

const ℝ=Real
const ℂ=Complex{Real}

abstract EuclideanSpace{F<:Union{ℝ,ℂ}} <: HilbertSpace{F}
# Elementary finite-dimensional space R^d or C^d with standard (Euclidean) inner product (i.e. orthonormal basis)

Base.conj(V::EuclideanSpace) = dual(V)
# The complex conjugate space of a Euclidean space is naturally isomorphic to its dual space

dual(V::EuclideanSpace{ℝ}) = V
# For a real Euclidean space, the dual space is naturally isomorphic to the vector space

# Functionality for extracting and iterating over elementary space
Base.length(V::ElementarySpace) = 1
Base.endof(V::ElementarySpace) = 1
Base.getindex(V::ElementarySpace, n::Integer) = (n == 1 ? V : throw(BoundsError()))

Base.start(V::ElementarySpace) = false
Base.next(V::ElementarySpace, state) = (V,true)
Base.done(V::ElementarySpace, state) = state

# CartesianSpace: standard real space R^N with euclidean structure (canonical orthonormal/cartesian basis)
include("spaces/cartesianspace.jl")

# ComplexSpace: standard complex space C^N with euclidean structure (canonical orthonormal basis)
include("spaces/complexspace.jl")

# GeneralSpace: a general elementary vector space (without inner product structure)
include("spaces/generalspace.jl")

# Graded elementary vector spaces
#---------------------------------
# Vector spaces which have a natural decomposition as a direct sum of several sectors

# Sector: hierarchy of types for labelling different sectors, e.g. irreps, quantum numbers, anyon types, ...
abstract Sector
# Abelian: sectors that have an abelian fusion structure
abstract Abelian <: Sector
*(s::Abelian)=s

# Explicit realizations
include("spaces/sectors.jl")

# UnitaryRepresentationSpace: Family of euclidean spaces that are graded corresponding to unitary representations
abstract UnitaryRepresentationSpace{G<:Sector} <: EuclideanSpace{ℂ}

sectortype{G}(V::UnitaryRepresentationSpace{G}) = G
sectortype{G}(::Type{UnitaryRepresentationSpace{G}}) = G
sectortype{S<:UnitaryRepresentationSpace}(::Type{S}) = sectortype(super(S))

#This is a useful shorthand for handling, e.g. the output of sectors()
Base.conj{N,S<:Sector}(nsec::NTuple{N,S}) = map(conj, nsec)

# AbelianSpace: general space that is graded by abelian sectors (i.e. one-dimensional representations)
include("spaces/abelianspace.jl")

# Auxiliary functionality:
#This returns the dimension of the degeneracy space specified by the irrep label assignment s.
_dim{G<:Sector,N}(spaces::NTuple{N,UnitaryRepresentationSpace{G}},s::NTuple{N,G})=(d=1;for n=1:N;d*=dim(spaces[n],s[n]);end;return d)

#This function returns all possible irrep label assignments for a space of indices
@generated function _sectors{G<:Sector,N}(spaces::NTuple{N,UnitaryRepresentationSpace{G}})
    quote
        numsectors=1
        @nexprs $N i->(s_i=collect(sectors(spaces[i]));numsectors*=length(s_i))
        sectorlist=Array(NTuple{N,G},numsectors)
        counter=0
        @nloops $N i d->1:length(s_d) begin
            counter+=1
            sectorlist[counter]=@ntuple $N k->s_k[i_k]
        end
        return sectorlist
    end
end

#This function finds the combinations of irrep labels that can habe nonzero value in an invariant tensor
@generated function _invariantsectors{G<:Abelian,N}(spaces::NTuple{N,AbelianSpace{G}})
    quote
        @nexprs $N i->(s_i=collect(sectors(spaces[i]))) #generates s_n = collect(sectors(spaces[i])) for n in 1:N
        sectorlist=Array(NTuple{N,G},0) #a sector is specified by N irrep labels (instances of G)
        @nloops $N i d->1:length(s_d) begin #generates N nested loops: for i_n in 1:length(s_n)
            sector=@ntuple $N k->s_k[i_k]   #creates an N-tuple: (s_1[i_1], s_2[i_2], ..., s_N[i_N])
            c=one(G)
            for n=1:N
                c*=sector[n] #fuses all the irrep labels in sector
            end
            if c==one(G) #if they fuse to the trivial irrep, we've found an allowed block: add to the list!
                push!(sectorlist,sector)
            end
        end
        return sectorlist
    end
end


text:
module RadauBig

using ForwardDiff
using StaticArrays
using LinearAlgebra
using DelimitedFiles
using Polynomials: Poly, polyval, polyder, coeffs
using PolynomialRoots: roots
using GenericLinearAlgebra: eigvals
using GenericSVD: svd

include("big.jl")
include("generate_butcher_table.jl")
include("save_radau_tables.jl")

export
    # big.jl
    big_eigen,

    # generate_butcher_table.jl
    find_real_eigenvalue,
    radau_butcher_table_plus,

    # save_radau_tables.jl
    save_radau_table_to_file

end


text:
function default_tls_debug(level, filename, number, msg)
    @debug(level, filename, number, msg)
end

function default_tls_rng()
    entropy = MbedTLS.Entropy()
    rng = MbedTLS.CtrDrbg()
    MbedTLS.seed!(rng, entropy)
    rng
end

"""
    amqps_configure(;
        cacerts = nothing,
        verify = MbedTLS.MBEDTLS_SSL_VERIFY_NONE,
        client_cert = nothing,
        client_key = nothing
    )

Creates and returns a configuration for making AMQPS connections.
- cacerts: A CA certificate file (or it's contents) to use for certificate verification.
- verify: Whether to verify server certificate. Default is false if cacerts is not provided and true if it is.
- client_cert and client_key: The client certificate and corresponding private key to use. Default is nothing (no client certificate). Values can either be the file name or certificate/key contents.
"""
function amqps_configure(;
    rng = default_tls_rng(),
    cacerts::Union{String,Nothing} = nothing,
    verify::Int64 = (cacerts === nothing) ? MbedTLS.MBEDTLS_SSL_VERIFY_NONE : MbedTLS.MBEDTLS_SSL_VERIFY_REQUIRED,
    client_cert::Union{String,Nothing} = nothing,
    client_key::Union{String,Nothing} = nothing,
    debug::Union{Function,Nothing} = nothing)

    conf = MbedTLS.SSLConfig()
    MbedTLS.config_defaults!(conf)
    MbedTLS.rng!(conf, rng)
    (debug === nothing) || MbedTLS.dbg!(conf, debug)

    if cacerts !== nothing
        if isfile(cacerts)
            # if it is a file name instead of certificate contents, read the contents
            cacerts = read(cacerts, String)
        end
        MbedTLS.ca_chain!(conf, MbedTLS.crt_parse(cacerts))
    end
    MbedTLS.authmode!(conf, verify)

    if (client_cert !== nothing) && (client_key !== nothing)
        if isfile(client_cert)
            # if it is a file name instead of certificate contents, read the contents
            client_cert = read(client_cert, String)
        end
        if isfile(client_key)
            client_key = read(client_key, String)
        end
        key = MbedTLS.PKContext()
        MbedTLS.parse_key!(key, client_key)
        MbedTLS.own_cert!(conf, MbedTLS.crt_parse(client_cert), key)
    end

    conf
end

function setup_tls(sock::TCPSocket, hostname::String, ssl_options::MbedTLS.SSLConfig)
    @debug("setting up TLS")

    ctx = MbedTLS.SSLContext()
    MbedTLS.setup!(ctx, ssl_options)
    MbedTLS.set_bio!(ctx, sock)
    MbedTLS.hostname!(ctx, hostname)
    MbedTLS.handshake(ctx)
    @debug("TLS setup done")

    BufferedTLSSocket(ctx)
end


text:
using Distributions, PDMats
using Test, LinearAlgebra

# RealInterval
r = RealInterval(1.5, 4.0)
@test minimum(r) == 1.5
@test maximum(r) == 4.0
@test extrema(r) == (1.5, 4.0)

@test partype(Gamma(1, 2)) == Float64
@test partype(Gamma(1.1, 2)) == Float64
@test partype(Normal(1//1, 2//1)) == Rational{Int}
@test partype(MvNormal(rand(Float32, 5), Matrix{Float32}(I, 5, 5))) == Float32

# special cases
@test partype(Kolmogorov()) == Float64
@test partype(Hypergeometric(2, 2, 2)) == Float64
@test partype(DiscreteUniform(0, 4)) == Float64

A = rand(1:10, 5, 5)
B = rand(Float32, 4)
C = 1//2
L = rand(Float32, 4, 4)
D = PDMats.PDMat(L * L')

# Ensure that utilities functions works with abstract arrays

@test isprobvec(GenericArray([1, 1, 1])) == false
@test isprobvec(GenericArray([1/3, 1/3, 1/3]))

# Positive definite matrix
M = GenericArray([1.0 0.0; 0.0 1.0])
# Non-invertible matrix
N = GenericArray([1.0 0.0; 1.0 0.0])

@test Distributions.isApproxSymmmetric(N) == false
@test Distributions.isApproxSymmmetric(M)


text:
#!/bin/env julia

# __precompile__()

using CSV
using Statistics
using DataFrames
using Printf

if length(ARGS) < 2
    println("give me input file name")
    exit(0)
end

input = ARGS

function process(input)
    df = CSV.read(input, DataFrame)
    p = df."Power(mAh)"
    println("power: $p")
    m = Statistics.mean(p)
    d = Statistics.stdm(p, m)
    @printf("power: %.2f ± %.2f\n", m, d)
end

@time process(input[1])
@time process(input[2])

max_stars_repo_name: alsam/scriptics-samples

max_issues_repo_name: alsam/scriptics-samples

max_forks_repo_name: alsam/scriptics-samples


text:
module VAN

include("sampler.jl")
include("arm.jl")
include("psa.jl")
include("exact.jl")
include("loss.jl")
include("utils.jl")
include("zygote_patch.jl")

end # module


text:
###############################################################################
#
# Initialize CanopyOpticals
#
###############################################################################
"""
    create_canopy_opticals( FT, rt_dim::RTDimensions)

Create a [`CanopyOpticals`](@ref) struct, given
- `FT` Floating number type
- `rt_dim` [`RTDimensions`](@ref) type struct
"""
function create_canopy_opticals(FT, rt_dim::RTDimensions)
    @unpack nAzi, nIncl, nLayer, nWL = rt_dim;

    return CanopyOpticals{FT}(nAzi   = nAzi  ,
                              nIncl  = nIncl ,
                              nLayer = nLayer,
                              nWL    = nWL   )
end

max_stars_repo_name: bbuman/Land

max_issues_repo_name: bbuman/Land

max_forks_repo_name: bbuman/Land


text:
### A Pluto.jl notebook ###
# v0.16.1

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 94f8e29e-ef91-11eb-1ae9-29bc46fa505a
begin
	using Pkg
    Pkg.activate(joinpath(Pkg.devdir(), "MLCourse"))
	using Plots, DataFrames, Random, CSV, MLJ, MLJLinearModels
end

# ╔═╡ 12d5824c-0873-49a8-a5d8-f93c73b633ae
using PlutoUI; PlutoUI.TableOfContents()

# ╔═╡ 20c5c7bc-664f-4c04-8215-8f3a9a2095c9
begin
    using MLCourse
    MLCourse.list_notebooks(@__FILE__)
end

# ╔═╡ 8217895b-b120-4b08-b18f-d921dfdddf10
md"# Linear Regression

## Wind speed prediction with one predictor
"

# ╔═╡ 9f84bcc5-e5ab-4935-9076-c19e9bd668e8
weather = CSV.read(joinpath(@__DIR__, "..", "data", "weather2015-2018.csv"), DataFrame);

# ╔═╡ 34e527f2-ef80-4cb6-be3a-bee055eca125
begin
    training_set1 = (X = (LUZ_pressure = weather.LUZ_pressure[1:end-5],),
                     y = weather.LUZ_wind_peak[6:end])
    m1 = machine(LinearRegressor(), training_set1.X, training_set1.y) |> fit!
end;

# ╔═╡ 006fc1eb-50d5-4206-8c87-53e873f158f4
begin
    scatter(training_set1.X.LUZ_pressure[1:10:end],
            training_set1.y[1:10:end], label = "data")
    plot!(fitted_linear_func(m1), label = "linear fit", w = 3)
end

# ╔═╡ e4712ebe-f395-418b-abcc-e10ada4b05c2
md"Let us inspect the results.
First we look at the fitted parameters.
We find that there is a negative correlation between the pressure in Luzern
and the wind speed 5 hours later in Luzern.
"

# ╔═╡ 8c9ae8f7-81b2-4d60-a8eb-ded5364fe0cc
fitted_params(m1)

# ╔═╡ f4f890b6-0ad4-4155-9321-15d673e15489
md"Next we predict the mean for different pressure values the distribution of
    wind speeds. In the probabilistic interpretation of supervised learning,
    standard linear regression finds conditional normal distributions with
    input-dependent mean ``\mu = \hat y = \theta_0 + \theta_1 x`` and
    constant standard deviation σ.
"

# ╔═╡ 9b62c374-c26e-4990-8ffc-790928e62e88
predict(m1, (LUZ_pressure = [930., 960., 990.],))

# ╔═╡ 7923a0a8-3033-4dde-91e8-22bf540c9866
md"We use the root-mean-squared-error (`rmse`) = ``\sqrt{\frac1n\sum_{i=1}^n(y_i - \hat y_i)^2}`` to evaluate our training error. To compute the test error, we use hourly data from 2019 to 2020.
"

# ╔═╡ 57f352dc-55ee-4e14-b68d-698938a97d92
rmse(predict(m1, training_set1.X), training_set1.y)

# ╔═╡ b0de002f-3f39-4378-8d68-5c4606e488b7
begin
    weather_test = CSV.read(joinpath(@__DIR__, "..", "data", "weather2019-2020.csv"), DataFrame);
    test_set1 = (X = (LUZ_pressure = weather_test.LUZ_pressure[1:end-5],),
                     y = weather_test.LUZ_wind_peak[6:end])
end

# ╔═╡ ce536f60-68b3-4901-bd6a-c96378054b12
rmse(predict(m1, test_set1.X), test_set1.y)

# ╔═╡ c65b81bd-395f-4461-a73b-3535903cb2d7
md"## Multiple Linear Regression

In multiple linear regression there are multiple predictors ``x_1, x_2, \ldots, x_p``
and the response is ``\hat y = \beta_0 + \beta_1 x_1 + \beta_2 x_2 + \cdots + \beta_p x_p
`` (we use here ``\beta`` instead of ``\theta`` for the parameters).
With ``p = 2`` predictors we can visualize linear regression as the plane
that is closest to the data. Use the sliders below to get a feeling for the
parameters. You can also change the viewing angle by clicking and dragging in
the figure."

# ╔═╡ 51c9ea74-3110-4536-a4af-7cc73b45a4a6
md"β₀ = $(@bind β₀ Slider(-1:.02:1, default = .4, show_value = true))

β₁ = $(@bind β₁ Slider(-1:.02:1, default = .2, show_value = true))

β₂ = $(@bind β₂ Slider(-1:.02:1, default = -.6, show_value = true))
"

# ╔═╡ 0f544053-1b7a-48d6-b18b-72092b124305
begin
    Random.seed!(3)
	X = DataFrame(X1 = randn(20), X2 = randn(20))
    f0(X1, X2, β₀, β₁, β₂) = β₀ + β₁ * X1 + β₂ * X2
    f0(β₀, β₁, β₂) = (X1, X2) -> f0(X1, X2, β₀, β₁, β₂)
	data_generator(X1, X2; β₀, β₁, β₂, σ = 0.8) = f0(X1, X2, β₀, β₁, β₂) + σ * randn()
	y = data_generator.(X.X1, X.X2, β₀ = .4, β₁ = .5, β₂ = -.6)
end;

# ╔═╡ d541a8cd-5aa4-4c2d-bfdf-5e6297bb65a8
begin
    plotly()
    p1 = scatter3d(X.X1, X.X2, y, markersize = 1,
                   xlims = (-3, 3), xlabel = "X1",
                   ylims = (-3, 3), ylabel = "X2",
                   zlims = (-4, 4), zlabel = "y", label = "data")
    wireframe!(-3:.1:3, -3:.1:3, f0(β₀, β₁, β₂),
               label = "function", title = "data & function", color = :green)
    plot_residuals!(X.X1, X.X2, y, f0(β₀, β₁, β₂))
    p2 = contour(-1:.1:1, -1:.1:1, (β₀, β₁) -> mean((β₀ .+ β₁ .* X.X1 .+ β₂ .* X.X2 .- y).^2), levels = 100, ylabel = "β₁", cbar = false, title = "loss")
    scatter!([β₀], [β₁], label = nothing)
    p3 = contour(-1:.1:1, -1:.1:1, (β₀, β₂) -> mean((β₀ .+ β₁ .* X.X1 .+ β₂ .* X.X2 .- y).^2), levels = 100, xlabel = "β₀", ylabel = "β₂")
    scatter!([β₀], [β₂], label = "current loss")
    plot(p1, plot(p2, p3, layout = (2, 1)), layout = (1, 2),
         size = (700, 400), legend = false)
end

# ╔═╡ da6462d8-3343-41d8-82dd-48770176d4ba
md"## Wind speed prediction with multiple predictors"

# ╔═╡ 753ec309-1363-485d-a2bd-b9fa100d9058
m2 = machine(LinearRegressor(), select(weather[1:end-5,:], Not([:LUZ_wind_peak, :time])),
             weather.LUZ_wind_peak[6:end]) |> fit!;

# ╔═╡ 618ef3c7-0fda-4970-88e8-1dac195545de
sort!(DataFrame(predictor = names(select(weather, Not([:LUZ_wind_peak, :time]))),
                value = fitted_params(m2).coefs), :value)

# ╔═╡ 2d25fbb6-dc9b-40ad-bdce-4c952cdad077
rmse(predict(m2, select(weather[1:end-5,:], Not([:LUZ_wind_peak, :time]))),
     weather.LUZ_wind_peak[6:end])

# ╔═╡ c9f10ace-3299-45fb-b98d-023a35dd405a
rmse(predict(m2, select(weather_test[1:end-5,:], Not([:LUZ_wind_peak, :time]))),
     weather_test.LUZ_wind_peak[6:end])

# ╔═╡ 99a371b2-5158-4c42-8f50-329352b6c1f2
md"# Error Decomposition

"

# ╔═╡ f10b7cad-eda3-4ec9-99ee-d43ed013a057
begin
    f(x) = sin(2x) + 2*(x - .5)^3 - .5x
    conditional_generator(x; n = 50) = f.(x) .+ .2*randn(n)
end;

# ╔═╡ 05354df5-a803-422f-87a3-1c56a34e8a48
f̂(x) = 0.1 + x

# ╔═╡ 9e61b4c3-1a9f-41a7-9882-25ed797a7b8d
expected_error(f, x) = mean((conditional_generator(x, n = 10^6) .- f(x)).^2);

# ╔═╡ c6a59b85-d031-4ad4-9e24-691494d08cde
expected_error(f̂, .1)

# ╔═╡ e50b8196-e804-473a-b3b5-e22fdb9d2f45
(f(.1) - f̂(.1))^2

# ╔═╡ f413ea94-36ca-4afc-8ca8-9a7e88101980
expected_error(f, .1)

# ╔═╡ dbf7fc72-bfd0-4c57-a1a9-fb5881e16e7e
let x = rand(100), grid = 0:.05:1
    gr()
    p1 = scatter(x, vcat(conditional_generator.(x, n = 1)...), label = "samples")
    plot!(f, label = "f")
    plot!(f̂, label = "f̂")
    p2 = plot(grid, expected_error.(f̂, grid), label = "expected error f̂", w = 3)
    plot!(grid, (f.(grid) .- f̂.(grid)).^2, label = "reducible error", w = 3)
    hline!([.2^2], label = "irreducible error", ylims = (0, .15), w = 3, xlabel = "x")
    plot(p1, p2, layout = (2, 1), legend = :right, ylabel = "y")
end

# ╔═╡ ad5b293d-c0f4-4693-84f4-88308639a501
md"# Logistic Regression

## Preparing the spam data

The text in our spam data set is already preprocessed. But we do not yet have
a format similar to our weather prediction data set with a fixed number ``p`` of
predictors for each email. In this section we create a very simple feature
representation of our emails:
1. We create a lexicon of words that are neither very frequent nor very rare.
2. For each email we count how often every word in this lexicon appears.
3. Our feature matrix will consist of ``n`` rows (one for each email) and ``p``
    predictors (one for each word in the lexicon) with ``x_{ij}`` measuring how
    often word ``j`` appears in document ``i``, normalized by the number of
    lexicon words in each email (such that the elements in every row sum to 1).
"

# ╔═╡ 210b977d-7136-407f-a1c9-eeea869d0312
begin
    spamdata = CSV.read(joinpath(@__DIR__, "..", "data", "spam.csv"), DataFrame)
    dropmissing!(spamdata) # remove entries without any text (missing values).
end

# ╔═╡ 72969aca-b203-4d83-8923-74e523aa1c01
import TextAnalysis: Corpus, StringDocument, DocumentTermMatrix, lexicon,
                     update_lexicon!, tf

# ╔═╡ 4cbb3057-01f4-4e80-9029-4e80d6c9e5e6
md"In the next cell we create the full lexicon of words appearing in the first
2000 emails. Each lexicon entry is of the form `\"word\" => count`."

# ╔═╡ c50c529f-d393-4854-b5ed-91e90d557d12
begin
    crps = Corpus(StringDocument.(spamdata.text[1:2000]))
    update_lexicon!(crps)
    lexicon(crps)
end

# ╔═╡ 72b50cee-d436-42ce-add9-07b0c012cb31
md"Now we select only those words of the full lexicon that appear at least 100
times and at most 10^3 times. These numbers are pulled out of thin air (like
all the design choses of this very crude feature engineering).
"

# ╔═╡ bf4110a9-31a4-48a3-bd6d-85c404d0e72d
begin
    small_lex = Dict(k => lexicon(crps)[k]
                     for k in findall(x -> 100 <= last(x) <= 10^3, lexicon(crps)))
    m = DocumentTermMatrix(crps, small_lex)
end

# ╔═╡ 534681d5-71d8-402a-b455-f491cfbb353e
begin
    spam_or_ham = coerce(String.(spamdata.label[1:2000]), Binary)
    normalized_word_counts = float.(DataFrame(tf(m), :auto))
end

# ╔═╡ ec1c2ea5-29ce-4371-be49-08798305ff50
Markdown.parse("Here we go: now we have a matrix of size
 $(join(size(normalized_word_counts), " x ")) as input and a vector of binary label as
 output. We will be able to use this as input in multiple logistic regression.
 For future usage we save this preprocessed representation of the spam data to
 a file.")

# ╔═╡ 681cb7b9-f041-4aea-907e-4d85135c005a
CSV.write(joinpath(dirname(pathof(MLCourse)), "..", "data", "spam_preprocessed.csv"),
          [normalized_word_counts DataFrame(spam_or_ham = spam_or_ham)])

# ╔═╡ f7117513-283f-4e32-a2a1-3594c794c94d
md"## Multiple Logistic Regression

In the top row of the figure below we see in two different ways (once as a 3D
plot and once as a contour plot) the probability of class A for the selected
parameter values. The bottom row shows samples (large points, red = class A)
obtained with this probability distribution and predictions (small points)
at decision threshold 0.5.
Play with the parameters to get a feeling for how they affect the probability
and the samples.

θ₀ = $(@bind θ₀ Slider(-3:3, default = 0, show_value = true))

θ₁ = $(@bind θ₁ Slider(-8:8, default = 3, show_value = true))

θ₂ = $(@bind θ₂ Slider(-8:8, default = 0, show_value = true))
"

# ╔═╡ fd4165dc-c3e3-4c4c-9605-167b5b4416da
md"## Confusion Matrix, ROC and AUC"

# ╔═╡ 7738c156-8e1b-4723-9818-fba364822171
md"s = $(@bind s Slider(-4:.1:4, default = 0, show_value = true))

seed = $(@bind seed Slider(1:100, show_value = true))

threshold = $(@bind threshold Slider(.01:.01:.99, default = 0.5, show_value = true))
"

# ╔═╡ 0fcfd7d2-6ea3-4c75-bad3-7d0fdd6fde11
begin
    logistic(x) = 1/(1 + exp(-x))
    logodds(p) = log(p/(1-p))
    function error_rates(x, y, t)
        P = sum(y)
        N = length(y) - P
        pos_pred = y[(logistic.(x) .> t)]
        TP = sum(pos_pred)
        FP = sum(1 .- pos_pred)
        FP/N, TP/P
    end
end;

# ╔═╡ 4f89ceab-297f-4c2c-9029-8d2d7fad084f
let f(x1, x2) = logistic(θ₀ + θ₁ * x1 + θ₂ * x2)
    p1 = wireframe(-3:.1:3, -3:.1:3, f, zlims = (0, 1))
    p2 = contour(-3:.1:3, -3:.1:3, f, contour_labels = true, levels = 20, cbar = false)
	plotly()
    samples = (X1 = 6 * rand(200) .- 3, X2 = 6 * rand(200) .- 3)
    labels = f.(samples.X1, samples.X2) .> rand(200)
    xgrid = MLCourse.grid(-3:.2:3, -3:.2:3, names = (:X1, :X2))
    scatter(xgrid.X1, xgrid.X2, color = (f.(xgrid.X1, xgrid.X2) .> .5) .+ 1,
            markersize = 2, markerstrokewidth = 0, label = nothing)
    p3 = scatter!(samples.X1, samples.X2, color = labels .+ 1, xlabel = "X1")
    plot(p1, p2, plot(), p3, layout = (2, 2), size = (700, 600),
         ylabel = "X2", legend = false)
end

# ╔═╡ 285c6bfc-5f29-46e0-a2c1-8abbec74501b
begin
    Random.seed!(seed)
    auc_samples_x = 2 * randn(200)
end;

# ╔═╡ c98524b5-d6b3-469c-82a1-7d231cc792d6
begin
    auc_samples_y = logistic.(2.0^s * auc_samples_x) .> rand(200)
    auc = [error_rates(2.0^s * auc_samples_x, auc_samples_y, t)
           for t in .01:.01:.99]
    push!(auc, (0., 0.))
    prepend!(auc, [(1., 1.)])
end;

# ╔═╡ 3336ab15-9e9b-44af-a7d5-1d6472241e62
let
    gr()
    p1 = scatter(auc_samples_x, auc_samples_y, markershape = :vline, label = nothing, color = :black)
    plot!(x -> logistic(2.0^s * x), color = :blue, label = nothing, xlims = (-8, 8))
    vline!([1/(2.0^s) * logodds(threshold)], w = 3, color = :red,
           label = nothing, xlabel = "x", ylabel = "y")
    p2 = plot(first.(auc), last.(auc), title = "ROC", label = nothing)
    fp, tp = auc[floor(Int, threshold * 100)]
    scatter!([fp], [tp], color = :red, xlims = (-.01, 1.01), ylims = (-.01, 1.01),
            labels = nothing, ylabel = "true positive rate",
            xlabel = "false positive rate")
    plot(p1, p2, size = (700, 400))
end

# ╔═╡ 62ad57e5-1366-4635-859b-ccdab2efd3b8
md"## Multiple Logistic Regression on the spam data"

# ╔═╡ 29e1d9ff-4375-455a-a69b-8dd0c2cac57d
m3 = fit!(machine(LogisticClassifier(penalty = :none),
                  normalized_word_counts,
                  spam_or_ham));

# ╔═╡ 1d1a24c6-c166-49a2-aa21-7acf50b55a66
predict(m3, normalized_word_counts)

# ╔═╡ 32bafa9e-a35e-4f54-9857-d269b47f95c3
confusion_matrix(predict_mode(m3, normalized_word_counts), spam_or_ham)

# ╔═╡ 4e4f4adf-364f-49b9-9391-5050a4c1286a
md"With our simple features, logistic regression can classify the training data
almost always correctly. Let us see how well this works for test data.
"

# ╔═╡ 50c035e6-b892-4157-a52f-824578366977
begin
    test_crps = Corpus(StringDocument.(spamdata.text[2001:4000]))
    test_input = float.(DataFrame(tf(DocumentTermMatrix(test_crps, small_lex)), :auto))
    test_labels = coerce(String.(spamdata.label[2001:4000]), Binary)
    confusion_matrix(predict_mode(m3, test_input), test_labels)
end

# ╔═╡ 21b66582-3fda-401c-9421-73ae2f455a75
predict_mode(m3, normalized_word_counts)

# ╔═╡ ba4b5683-5932-415e-8772-8b3eef5eb63d
md"We save also the test data for future usage."

# ╔═╡ a30578dd-aecb-46eb-b947-f009282cf2fc
md"Let us evaluate the fit in terms of commonly used losses for binary classification."

# ╔═╡ 8ed39cdc-e99e-48ff-9973-66df41aa0f78
function losses(machine, input, response)
    (loglikelihood = -sum(log_loss(predict(machine, input), response)),
     misclassification_rate = mean(predict_mode(machine, input) .!= response),
     accuracy = accuracy(predict_mode(machine, input), response),
     auc = MLJ.auc(predict(machine, input), response)
	)
end;

# ╔═╡ dd463687-b73d-4e70-b2cf-97a56a0ad409
losses(m3, normalized_word_counts, spam_or_ham)

# ╔═╡ 57dcadc0-2da2-4521-aeaf-6fd01f4bd82b
spam_or_ham

# ╔═╡ 935adbcd-48ab-4a6f-907c-b04137ca3abe
losses(m3, test_input, test_labels)

# ╔═╡ 8b0451bf-59b0-4e71-be84-549e23b5bfe7
md"""# Exercises

## Conceptual

1. Suppose we have a data set with three predictors, ``X_1`` = Final Grade, ``X_2`` = IQ, ``X_3`` = Level (1 for College and 0 for High School).  The response is starting salary after graduation (in thousands of dollars). Suppose we use least squares to fit the model, and get ``\hat\beta_0 = 25, \hat\beta_1 = 2, \hat\beta_2 = 0.07, \hat\beta_3 = 15``.
   - Which answer is correct, and why?
      - For a fixed value of IQ and Final Grade, high school graduates earn more, on average, than college graduates.
      - For a fixed value of IQ and Final Grade, college graduates earn more, on average, than high school graduates.
   - Predict the salary of a college graduate with IQ of 110 and a Final Grade of 4.0.
2. Suppose we collect data for a group of students in a machine learning class with variables ``X_1 =`` hours studied, ``X_2 =`` grade in statistics class, and ``Y =`` receive a 6 in the machine learning class. We fit a logistic regression and produce estimated coefficients, ``\hat{\beta}_0 = -6``, ``\hat{\beta}_1 = 0.025``, ``\hat{\beta}_2 = 1``.
   - Estimate the probability that a student who studies for 75 hours and had a 4 in the statistics class gets a 6 in the machine learning class.
   - How many hours would the above student need to study to have a 50% chance of getting an 6 in the machine learning class?
3. In this exercise we will derive the loss function implicitly defined by maximum likelihood estimation of the parameters in a classification setting with multiple classes. Remember that the input ``f(x)`` of the softmax function ``s`` is a vector-valued function. Here we assume a linear function ``f`` and write the ``i``th component of this function as ``f_i(x) = \theta_{i0} + \theta_{i1}x_1 + \cdots + \theta_{ip}x_p``. Note that each component ``i`` has now its own parameters ``\theta_{i0}`` to ``\theta_{ip}``. Using matrix multiplication we can also write ``f(x) = \theta x`` where ``\theta = \left(\begin{array}{ccc}\theta_{10} & \cdots & \theta_{1p}\\\vdots & \ddots & \cdots\\\theta_{K0} & \cdots & \theta_{Kp}\end{array}\right)`` is a ``K\times(p+1)`` dimensional matrix and ``x = (1, x_1, x_2, \ldots, x_p)`` is a column vector of length ``p+1``.
    - Write the log-likelihood function for a classification problem with ``K`` classes. *Hint*: to simplify the notation we can use the convention ``s_y(f(x)) = P(y|x)`` to write the conditional probability of class ``y`` given input ``x``. This convention makes sense when the classes are identified by the integers ``1, 2, \ldots, K``; in this case ``s_y(f(x))`` is the ``y``th component of ``s(f(x))``. Otherwise we would could specify a mapping from classes ``C_1, C_2, \ldots, C_K`` to the integers ``1, 2, \ldots, K`` for this convention to make sense.
    - Assume now ``K = 3`` and ``p = 2``. Explicitly write the loss function for the training set ``\mathcal D = ((x_1 = (0, 0), y_1 = C), (x_2 = (3, 0), y_2 = A), (x_3 = (0, 2), y_3 = B))``.
    - Assume ``K = 2`` and ``p = 1`` and set ``\theta_{20} = 0`` and ``\theta_{21} = 0``. Show that we recover standard logistic regression in this case. *Hint*: show that ``s_1(f(x)) = \sigma(f_1(x))`` and ``s_2(f(x)) = 1 - \sigma(f_1(x))``, where ``s`` is the softmax function and ``\sigma(x) = 1/(1 + e^{-x})`` is the logistic function.
    - Show that one can always set ``\theta_{K0}, \theta_{K1}, \ldots, \theta_{Kp}`` to zero. *Hint* Show that the softmax function with the transformed parameters ``\tilde\theta_{ij}=\theta_{ij} - \theta_{Kj}`` has the same value as the softmax function in the original parameters.

## Applied
1. In the multiple linear regression of the weather data set above we used all
   available predictors. We do not know if all of them are relevant. In this exercise our aim is to find models with fewer predictors and quantify the loss in prediction accuracy.
    - Systematically search for the model with at most 2 predictors that has the lowest test rmse. *Hint* write a function `train_and_evaluate` that takes the training and the test data as input as well as an array of two predictors; remember that `data[:, [\"A\", \"B\"]]` returns a sub-dataframe with columns \"A\" and \"B\". This function should fit a `LinearRegressor` on the training set with those two predictors and return the test rmse for the two predictors. To get a list of all pairs of predictors you can use something like `predictors = setdiff(names(train), ["time", "LUZ_wind_peak"]); predictor_pairs = [[p1, p2] for p1 in predictors, p2 in predictors if p1 != p2 && p1 > p2]`
    - How much higher is the test error compared to the fit with all available predictors?
    - How many models did you have to fit to find your result above?
    - How many models would you have to fit to find the best model with at most 5 predictors? *Hint* the function `binomial` may be useful.
2. In this exercise we perform linear classification of the MNIST handwritten digits
   dataset.
   - Load the MNIST data set with `using OpenML; mnist = OpenML.load(554) |> DataFrame; dropmissing!(mnist);`
   - Usually the first 60'000 images are taken as training set, but for this exercise I recommend to use fewer rows, e.g. the first 5000.
   - Scale the input values to the interval [0, 1) with `mnist[:, 1:784] ./= 255`
   - Fit a `MLJLinearModels.MultinomialClassifier(penalty = :none)` to the data. Be patient! This can take a few minutes.
   - Compute the misclassification rate and the confusion matrix on the training set.
   - Use as test data rows 60001 to 70000 and compute the misclassification rate
     and the confusion matrix on this test set.
   - Plot some of the wrongly classified training and test images.
     Are they also difficult for you to classify?
"""

# ╔═╡ 7f08fcaa-000d-422d-80b4-e58a2f489d74
MLCourse.footer()

# ╔═╡ Cell order:
# ╟─12d5824c-0873-49a8-a5d8-f93c73b633ae
# ╠═94f8e29e-ef91-11eb-1ae9-29bc46fa505a
# ╟─8217895b-b120-4b08-b18f-d921dfdddf10
# ╠═9f84bcc5-e5ab-4935-9076-c19e9bd668e8
# ╠═34e527f2-ef80-4cb6-be3a-bee055eca125
# ╠═006fc1eb-50d5-4206-8c87-53e873f158f4
# ╟─e4712ebe-f395-418b-abcc-e10ada4b05c2
# ╠═8c9ae8f7-81b2-4d60-a8eb-ded5364fe0cc
# ╟─f4f890b6-0ad4-4155-9321-15d673e15489
# ╠═9b62c374-c26e-4990-8ffc-790928e62e88
# ╟─7923a0a8-3033-4dde-91e8-22bf540c9866
# ╠═57f352dc-55ee-4e14-b68d-698938a97d92
# ╠═b0de002f-3f39-4378-8d68-5c4606e488b7
# ╠═ce536f60-68b3-4901-bd6a-c96378054b12
# ╟─c65b81bd-395f-4461-a73b-3535903cb2d7
# ╟─51c9ea74-3110-4536-a4af-7cc73b45a4a6
# ╟─d541a8cd-5aa4-4c2d-bfdf-5e6297bb65a8
# ╟─0f544053-1b7a-48d6-b18b-72092b124305
# ╟─da6462d8-3343-41d8-82dd-48770176d4ba
# ╠═753ec309-1363-485d-a2bd-b9fa100d9058
# ╠═618ef3c7-0fda-4970-88e8-1dac195545de
# ╠═2d25fbb6-dc9b-40ad-bdce-4c952cdad077
# ╠═c9f10ace-3299-45fb-b98d-023a35dd405a
# ╟─99a371b2-5158-4c42-8f50-329352b6c1f2
# ╠═f10b7cad-eda3-4ec9-99ee-d43ed013a057
# ╠═05354df5-a803-422f-87a3-1c56a34e8a48
# ╠═9e61b4c3-1a9f-41a7-9882-25ed797a7b8d
# ╠═c6a59b85-d031-4ad4-9e24-691494d08cde
# ╠═e50b8196-e804-473a-b3b5-e22fdb9d2f45
# ╠═f413ea94-36ca-4afc-8ca8-9a7e88101980
# ╟─dbf7fc72-bfd0-4c57-a1a9-fb5881e16e7e
# ╟─ad5b293d-c0f4-4693-84f4-88308639a501
# ╠═210b977d-7136-407f-a1c9-eeea869d0312
# ╠═72969aca-b203-4d83-8923-74e523aa1c01
# ╟─4cbb3057-01f4-4e80-9029-4e80d6c9e5e6
# ╠═c50c529f-d393-4854-b5ed-91e90d557d12
# ╟─72b50cee-d436-42ce-add9-07b0c012cb31
# ╠═bf4110a9-31a4-48a3-bd6d-85c404d0e72d
# ╠═534681d5-71d8-402a-b455-f491cfbb353e
# ╟─ec1c2ea5-29ce-4371-be49-08798305ff50
# ╠═681cb7b9-f041-4aea-907e-4d85135c005a
# ╟─f7117513-283f-4e32-a2a1-3594c794c94d
# ╟─4f89ceab-297f-4c2c-9029-8d2d7fad084f
# ╟─fd4165dc-c3e3-4c4c-9605-167b5b4416da
# ╟─7738c156-8e1b-4723-9818-fba364822171
# ╠═0fcfd7d2-6ea3-4c75-bad3-7d0fdd6fde11
# ╟─3336ab15-9e9b-44af-a7d5-1d6472241e62
# ╟─285c6bfc-5f29-46e0-a2c1-8abbec74501b
# ╟─c98524b5-d6b3-469c-82a1-7d231cc792d6
# ╟─62ad57e5-1366-4635-859b-ccdab2efd3b8
# ╠═29e1d9ff-4375-455a-a69b-8dd0c2cac57d
# ╠═1d1a24c6-c166-49a2-aa21-7acf50b55a66
# ╠═32bafa9e-a35e-4f54-9857-d269b47f95c3
# ╟─4e4f4adf-364f-49b9-9391-5050a4c1286a
# ╠═50c035e6-b892-4157-a52f-824578366977
# ╠═21b66582-3fda-401c-9421-73ae2f455a75
# ╟─ba4b5683-5932-415e-8772-8b3eef5eb63d
# ╟─a30578dd-aecb-46eb-b947-f009282cf2fc
# ╠═8ed39cdc-e99e-48ff-9973-66df41aa0f78
# ╠═dd463687-b73d-4e70-b2cf-97a56a0ad409
# ╠═57dcadc0-2da2-4521-aeaf-6fd01f4bd82b
# ╠═935adbcd-48ab-4a6f-907c-b04137ca3abe
# ╟─8b0451bf-59b0-4e71-be84-549e23b5bfe7
# ╟─20c5c7bc-664f-4c04-8215-8f3a9a2095c9
# ╟─7f08fcaa-000d-422d-80b4-e58a2f489d74

max_stars_repo_name: jbrea/MLCourse

max_issues_repo_name: davibarreira/MLCourse

max_forks_repo_name: davibarreira/MLCourse


text:
"""
`Inventory` represent the inventory of an inv in the BOM of an `InventorySystem`.
During the top-bottom activation, an `Inventory` receives pull orders from `Assembly`s or `Market`s.
During the bottom-top activation, an `Inventory` resolves the pull orders by pushing quantities from its on-hand inventory to pull-order issuers.
No prioritization is considered with respect to external demand's cost penalties. Two alternatives could be considered later:
1) adding a prioritization "action" to the action space of the `inv`;
2) creating a _recourse_ framework where the agent must decide order quantities before the demand realizes but can allocate after.
Consequently the inventory models cannot fully model environments where an inv has both external and internal demands or more than one external demand.
"""
mutable struct Inventory{Dl<:Distribution, F}
    holding_cost::F
    onhand::Float64
    capacity::Float64
    onhand_reset::Dl
    pull_orders::IdDict{Any, Float64}
    on_hand_log::Vector{Float64}
    stockout_log::Vector{Float64}
    fillrate_log::Vector{Float64}
    cost_log::Vector{Float64}
    name::String
end

function Inventory(holding_cost, onhand::NumDist; capacity::Number = Inf, name = "inventory")
    @assert hasmethod(holding_cost, Tuple{Inventory}) "holding cost must have a method with `(::Inventory)` arguments"
    dl = parametrify(onhand)
    Inventory{typeof(dl), typeof(holding_cost)}(holding_cost, min(capacity, Float64(rand(dl))), Float64(capacity), dl,  IdDict(), zeros(0), zeros(0), zeros(0), zeros(0), name)
end

function Inventory(holding_cost::Number, onhand::NumDist; capacity::Number = Inf, name = "inventory")
    Inventory(LinearHoldingCost(holding_cost), onhand, capacity = capacity, name = name)
end

RLBase.state(inv::Inventory) = [inv.onhand]
state_size(::Inventory) = 1
print_state(inv::Inventory) = [inv.name*" on hand" => inv.onhand]

function pull!(inv::Inventory, quantity::Number, issuer)
    push!(inv.pull_orders, issuer => quantity)
    return nothing
end

function Base.push!(inv::Inventory, quantity, source)
    inv.onhand += min(quantity, inv.capacity-inv.onhand)
    nothing
end

function dispatch!(inv::Inventory)
    sumorders = sum(values(inv.pull_orders))
    if sumorders != 0
        proportion = min(one(sumorders), inv.onhand/sumorders)
        for (issuer, quantity) in inv.pull_orders
            push!(issuer, quantity*proportion, inv)
        end
        push!(inv.stockout_log, (1-proportion)*sumorders)
        push!(inv.fillrate_log, proportion)
        inv.onhand -= proportion*sumorders
    else
        push!(inv.stockout_log, 0.0)
        push!(inv.fillrate_log, 1.0)
        for (issuer, quantity) in inv.pull_orders
            push!(issuer, 0.0, inv)
        end
    end
    nothing
end

"""
    reward(inv::Inventory)

1. Resolve the pull orders by pushing quantities from its on-hand inventory to pull-order issuers with proportional distribution in case of shortage.
2. Erase pull_orders.
3. Returns a reward with respect to minus its `cost` function.
"""
function reward!(inv::Inventory)
    push!(inv.on_hand_log, inv.onhand)
    empty!(inv.pull_orders)
    cost = inv.holding_cost(inv)
    push!(inv.cost_log, cost)
    return -cost
end

"""
    reset!(inv::Inventory)

Randomizes inv's on-hand inventory with respect to its initial distribution"
"""
function RLBase.reset!(inv::Inventory)
    empty!(inv.cost_log)
    empty!(inv.on_hand_log)
    empty!(inv.fillrate_log)
    empty!(inv.stockout_log)
    inv.onhand = min(inv.capacity, rand(inv.onhand_reset))
    return nothing
end

inventory_position(inv::Inventory) = inv.onhand

"""
    LinearHoldingCost(h)

A simple linear holding cost function with h as the cost per unit of `onhand` inventory.
"""
mutable struct LinearHoldingCost{T}
    h::T
end
(f::LinearHoldingCost)(inv::Inventory) = f.h*inv.onhand

Base.show(io::IO, inv::Inventory{Dl, F}) where {Dl,F} = print(io, "Inventory($(inv.name), $F)")

max_stars_repo_path: src/elementary blocks/inventory.jl

max_issues_repo_path: src/elementary blocks/inventory.jl

max_forks_repo_path: src/elementary blocks/inventory.jl


text:
### A Pluto.jl notebook ###
# v0.14.6

using Markdown
using InteractiveUtils

# ╔═╡ 2f761b20-bde2-11eb-190f-077bd9e0c638
begin
    import Pkg
    Pkg.activate(mktempdir())
    Pkg.add([
        Pkg.PackageSpec(name="Plots", version="1"),
        Pkg.PackageSpec(name="PlutoUI", version="0.7"),
	    Pkg.PackageSpec(name="DifferentialEquations", version="6"),
	    Pkg.PackageSpec(name="Parameters", version="0.12"),
        Pkg.PackageSpec(name="LabelledArrays", version="1"),
        Pkg.PackageSpec(name="Setfield", version="0.7"),
    ])
    using Plots, PlutoUI, LinearAlgebra, DifferentialEquations, Parameters, LabelledArrays, Setfield

	Plots.gr(fmt=:png, lw=2)

	PlutoUI.TableOfContents()
end

# ╔═╡ 11f28aab-79d2-4fa3-a9f8-34df3e581c0c
md"""
# Problem 3.7.5 Reduced Michaelis-Menten
"""

# ╔═╡ 708d922b-5c70-42b1-b065-4e47c3fba282
function model!(du, u, p, t)
	@unpack V0, VM1, VM2, VM3, KM1, KM2, KM3, _mm = p
	@unpack s1, s2, s3 = u
	v1 = VM1 * _mm(s1, KM1)
	v2 = VM2 * _mm(s2, KM2)
	v3 = VM3 * _mm(s3, KM3)
	du.s1 = V0 - v1
	du.s2 = v1 - v2
	du.s3 = v2 - v3
	return du
end

# ╔═╡ 14836120-7d01-4e59-8869-3b387fccf461
u0 = LVector(s1=0.3, s2=0.2, s3=0.1)

# ╔═╡ 89ce6bbb-ec25-454a-90f0-c45d4ec8b217
tend = 2.0

# ╔═╡ 0e967887-f222-4aa0-8874-1900b5aa69a4
# second set of ICs
u02 = LVector(s1=6.0, s2=4.0, s3=4.0)

# ╔═╡ bd13dc0b-929e-4270-a010-0da4a6e7b8ea
md"""
## Appendix
Running environment and some auxillary functions.

**You can use [this helper tool](https://fonsp.com/article-test-3/pkghelper.html) to generate these commands!**
"""

# ╔═╡ ebc7e456-ff70-4c51-ac5a-81daff88e287
begin
	# Convenience functions
    hill(x, k) = x / (x + k)
    hill(x, k, n) = hill(x^n, k^n)
    exprel(x) = ifelse(x≈zero(x), one(x), x / expm1(x))
	reducedmm(x, k) = x / k
end

# ╔═╡ b04b6578-c363-4c03-bc26-b8c608676fee
p = (V0 = 2.0,
	VM1 = 9.0,
	VM2 = 12.0,
	VM3 = 15.0,
	KM1 = 1.0,
	KM2 = 0.4,
	KM3 = 3.0,
	_mm = hill)

# ╔═╡ 52f29c2e-45ad-4d1e-820b-0103d80248a9
sol1 = solve(ODEProblem(model!, u0, tend, p));

# ╔═╡ cb3f07ea-a0a4-40f0-83d4-6ffc5cb5ddf6
plot(sol1, ylims=(0.0, 0.8),
     title="Problem 3.7.5 (1)",
     xlabel="Time (arbitrary units)",
     ylabel="Concentration (arbitrary units)")

# ╔═╡ e0081e60-2cd0-48fa-9049-e816666293b6
sol2 = solve(ODEProblem(model!, u02, 4.0, p));

# ╔═╡ b1883077-bbc2-409f-9f10-40d8c7b84868
plot(sol2, ylims=(0.0, 6.0),
     title="Problem 3.7.5 (2)",
     xlabel="Time (arbitrary units)",
     ylabel="Concentration (arbitrary units)")

# ╔═╡ 66444666-0be1-4104-8974-2133aca8e6d0
# Reduced model
param2 = @set p._mm = reducedmm

# ╔═╡ df7fc593-93a6-46fb-9a1a-0e7e488d7caa
sol3 = solve(ODEProblem(model!, u0, tend, param2));

# ╔═╡ 8ff61c0e-b725-48f4-88b1-ad06432105a9
begin
p3 = plot(sol1, ylims=(0.0, 0.8),
     title="Problem 3.7.5 (1) (full vs reduced)",
     xlabel="Time (arbitrary units)",
     ylabel="Concentration (arbitrary units)",
     labels=["S1 " "S2 " "S3 "], ls=:dash)
	plot!(p3, sol3, labels=["S1 (reduced)" "S2 (reduced)" "S3 (reduced)"] )
	p3
end

# ╔═╡ dad7528c-f6d0-4ace-b4e0-eb4c0f9ff6d1
sol4 = solve(ODEProblem(model!, u02, 4.0, param2));

# ╔═╡ b1e5b71d-7150-4cc5-981a-bb0f5910c584
begin
p4 = plot(sol2, ylims=(0.0, 8.0),
     title="Problem 3.7.5 (2) (full vs reduced)",
     xlabel="Time (arbitrary units)",
     ylabel="Concentration (arbitrary units)",
     labels=["S1 " "S2 " "S3 "], ls=:dash)
	plot!(p4, sol4, labels=["S1 (reduced)" "S2 (reduced)" "S3 (reduced)"] )

	p4
end

# ╔═╡ Cell order:
# ╠═11f28aab-79d2-4fa3-a9f8-34df3e581c0c
# ╠═708d922b-5c70-42b1-b065-4e47c3fba282
# ╠═14836120-7d01-4e59-8869-3b387fccf461
# ╠═89ce6bbb-ec25-454a-90f0-c45d4ec8b217
# ╠═b04b6578-c363-4c03-bc26-b8c608676fee
# ╠═52f29c2e-45ad-4d1e-820b-0103d80248a9
# ╠═cb3f07ea-a0a4-40f0-83d4-6ffc5cb5ddf6
# ╠═0e967887-f222-4aa0-8874-1900b5aa69a4
# ╠═e0081e60-2cd0-48fa-9049-e816666293b6
# ╠═b1883077-bbc2-409f-9f10-40d8c7b84868
# ╠═66444666-0be1-4104-8974-2133aca8e6d0
# ╠═df7fc593-93a6-46fb-9a1a-0e7e488d7caa
# ╠═8ff61c0e-b725-48f4-88b1-ad06432105a9
# ╠═dad7528c-f6d0-4ace-b4e0-eb4c0f9ff6d1
# ╠═b1e5b71d-7150-4cc5-981a-bb0f5910c584
# ╠═bd13dc0b-929e-4270-a010-0da4a6e7b8ea
# ╠═ebc7e456-ff70-4c51-ac5a-81daff88e287
# ╠═2f761b20-bde2-11eb-190f-077bd9e0c638

max_stars_repo_name: stevengogogo/pluto-notebooks

max_issues_repo_name: stevengogogo/pluto-notebooks

max_forks_repo_name: stevengogogo/pluto-notebooks


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

using ReTest

FluxTraining.runtests()

module FluxTrainingTests

include("./imports.jl")

include("./metrics.jl")
include("./protected.jl")
include("./training.jl")
include("./callbacks/stoponnanloss.jl")
include("./callbacks/conditional.jl")
include("./callbacks/logging.jl")
include("./callbacks/recorder.jl")
include("./callbacks/scheduler.jl")
include("./callbacks/checkpointer.jl")
include("./callbacks/garbagecollect.jl")
include("./callbacks/sanitycheck.jl")
include("./callbacks/earlystopping.jl")
include("./callbackutils.jl")

end

FluxTrainingTests.runtests()


text:
using LinearAlgebra, Arpack, IterativeSolvers
using Random, Distributions

struct Schur{
    T,
    M<:AbstractMatrix{T},
    N<:AbstractMatrix{T},
    O<:AbstractMatrix{T},
    P<:AbstractMatrix{T}
}
    W::M
    Ku::N
    Kux::O
    Kx::P
    ϵ::T
end

function Matrix(sc::Schur)
    return sc.Kx .- sc.W' * (sc.Kux .- sc.Ku * sc.W) .- sc.Kux' * sc.W + sc.ϵ * I
end
eltype(::Schur{T}) where T = T
size(sc::Schur, args...) = size(sc.Kx, args...)
function mul!(y, sc::Schur, x, α, β)
    p, n = size(sc.W)
    v1 = similar(x, p)
    v2 = similar(x, p)

    y .= y .* β .+ x .* (α * sc.ϵ)

    mul!(y, sc.Kx, x, α, 1)

    mul!(v1, sc.W, x)
    mul!(v2, sc.Ku, v1)

    mul!(y, sc.W', v2, α, 1)
    mul!(y, sc.Kux', v1, -α, 1)

    mul!(v1, sc.Kux, x)

    mul!(y, sc.W', v1, -α, 1)

    return y
end

function lanczos_arpack(A, k, v; maxiter, tol)
    T = eltype(A)
    n = size(A, 1)
    mulA! = (y, x) -> mul!(y, A, x)
    id = x -> x
    # in: (T, mulA!, mulB, solveSI, n, issym, iscmplx, bmat,
    #            nev, ncv, whichstr, tol, maxiter, mode, v0)
    # out: (resid, v, ldv, iparam, ipntr, workd, workl, lworkl, rwork, TOL)
    out = Arpack.aupd_wrapper(T, mulA!, id, id, n, true, false, "I",
                       1, k, "LM", tol, maxiter, 1, v)

    α = out[7][k+1:2*k-1]
    β = out[7][2:k-1]

    return out[2], α, β, out[1]
end

function _lanczos_logdet!(z, acc, sc, k; maxiter, tol, nsamples)
    for i in 1:nsamples
        rand!(Normal(), z)
        z .= sign.(z)
        Q, α, β, resid = lanczos_arpack(sc, k, z; maxiter=maxiter, tol=tol)
        T = SymTridiagonal(α, β)
        Λ = eigen(T)
        wts = Λ.vectors[1, :] .^ 2 .* norm(z) ^ 2
        acc += dot(wts, log.(Λ.values))
    end
    return acc / nsamples
end

function lanczos_logdet(sc, k; maxiter, tol, nsamples)
    z = Vector{eltype(sc)}(undef, size(sc, 1))
    return _lanczos_logdet!(z, 0, sc, k;
                            maxiter=maxiter, tol=tol, nsamples=nsamples)
end

max_stars_repo_name: dbindel/btg

max_issues_repo_name: dbindel/btg

max_forks_repo_name: dbindel/btg


text:
## J e W e L s
function JWL(name::AbstractString)
    Pkg.project(Pkg.Types.Context(),name,tempdir())
    folder = joinpath(tempdir(),name)
    open(joinpath(folder,"include.jl"), "w") do f
        Base.write(f, """

        [compat]
        julia = "1.6"
        """ * string(VERSION))
    end
    open(joinpath(folder,"include.jl"), "w") do f
        Base.write(f, """
        # Use this file to `import` any relevant packages*,
        # and to define your own functions for use in the workbook.
        # It will be `include`d when the jwl workbook is opened.
        # * (don't forget to add packages to the `Project.toml` too)

        """)
    end
    return JWL(name, folder, Workbook())
end

Base.show(io::IO, j::JWL) = Base.print(io, "JWL " * string(j.name) * " with " * string(j.wb) * " in " * dirname(j.folder))


text:
using BayesianMixtures
using Statistics

n_total = 1000  # total number of MCMC sweeps to run

# _____________________________________________________________________________________
# Check samplers on univariate data

n = 50
x = randn(n)
for mode_ in ["Normal","NormalNonoptimized"]
	for model_type_ in ["MFM","DPM"]
		# Run MCMC sampler
		options = BayesianMixtures.options(mode_,model_type_,x,n_total)
		result = BayesianMixtures.run_sampler(options)
	end
end

# _____________________________________________________________________________________
# Check samplers on multivariate data

d = 3
n = 50
x = randn(n,d)
x = (x .- mean(x,dims=1))./std(x,dims=1; corrected=false)
x = [x[i,:] for i=1:n]
for mode_ in ["MVN","MVNaaC","MVNaaN","MVNaaRJ"]
	for model_type_ in ["MFM","DPM"]
		if model_type_=="DPM" && mode_=="MVNaaRJ"; continue; end
		# Run MCMC sampler
		options = BayesianMixtures.options(mode_,model_type_,x,n_total)
		result = BayesianMixtures.run_sampler(options)
	end
end

# _____________________________________________________________________________________
# Check other functions using univariate data

# Settings
mode = "Normal"
model_type = "MFM"
n = 100
B = BayesianMixtures

# Generate data and run sampler
x = randn(n)
options = B.options(mode,model_type,x,n_total)
result = B.run_sampler(options)

# Save and load
B.save_result("temp.jld",result)
result2 = B.load_result("temp.jld")

B.t_running(result)
B.t_posterior(result)
B.k_posterior(result)

B.plot_density_estimate(result)

k,a,w,theta = B.sample_mixture_parameters(result,50)
loglik = B.log_likelihoods(x,w,theta,mode)

B.open_figure(1)
B.traceplot(result.t); P.pause(0.5); P.clf()
B.traceplot_timewise(result,1); P.pause(0.5); P.clf()
B.plot_t_running(result); P.pause(0.5); P.clf()
B.plot_autocorrelation(result,20); P.pause(0.5); P.clf()
B.plot_t_posterior(result); P.pause(0.5); P.clf()
B.plot_t_posterior_average([result]); P.pause(0.5); P.clf()
B.plot_k_posterior(result); P.pause(0.5); P.clf()
B.plot_k_posterior_average([result]); P.pause(0.5); P.clf()
B.plot_similarity_matrix(result); P.pause(0.5); P.clf()
B.plot_clusters(x,result.z[:,end]); P.pause(0.5); P.clf()
B.rug_plot(x); P.pause(0.5); P.clf()
B.plot_density_estimate(result); P.pause(0.5); P.clf()

# _____________________________________________________________________________________
# Check other functions using multivariate data

# Settings
mode = "MVN"
model_type = "MFM"
d = 2
n = 100
B = BayesianMixtures

# Generate data and run sampler
x = randn(n,d)
x = (x .- mean(x,dims=1))./std(x,dims=1; corrected=false)
x = [x[i,:] for i=1:n]
options = B.options(mode,model_type,x,n_total)
result = B.run_sampler(options)

# Save and load
B.save_result("temp.jld",result)
result2 = B.load_result("temp.jld")

B.t_running(result)
B.t_posterior(result)
B.k_posterior(result)

B.plot_density_estimate(result)

k,a,w,theta = B.sample_mixture_parameters(result,50)
loglik = B.log_likelihoods(x,w,theta,mode)

B.open_figure(1)
B.traceplot(result.t); P.pause(0.5); P.clf()
B.traceplot_timewise(result,1); P.pause(0.5); P.clf()
B.plot_t_running(result); P.pause(0.5); P.clf()
B.plot_autocorrelation(result,20); P.pause(0.5); P.clf()
B.plot_t_posterior(result); P.pause(0.5); P.clf()
B.plot_t_posterior_average([result]); P.pause(0.5); P.clf()
B.plot_k_posterior(result); P.pause(0.5); P.clf()
B.plot_k_posterior_average([result]); P.pause(0.5); P.clf()
B.plot_similarity_matrix(result); P.pause(0.5); P.clf()
B.plot_clusters(x,result.z[:,end]); P.pause(0.5); P.clf()
B.plot_density_estimate(result); P.pause(0.5); P.clf()


text:
# Lower-level interface to the zlib library.

if is_windows()
    const zlib = "zlib1"
else
    const zlib = "libz"
end

# Constants
# ---------

# Flush values
const Z_NO_FLUSH      = Cint(0)
const Z_PARTIAL_FLUSH = Cint(1)
const Z_SYNC_FLUSH    = Cint(2)
const Z_FULL_FLUSH    = Cint(3)
const Z_FINISH        = Cint(4)
const Z_BLOCK         = Cint(5)
const Z_TREES         = Cint(6)

# Return codes
const Z_OK            = Cint(0)
const Z_STREAM_END    = Cint(1)
const Z_NEED_DICT     = Cint(2)
const Z_ERRNO         = Cint(-1)
const Z_STREAM_ERROR  = Cint(-2)
const Z_DATA_ERROR    = Cint(-3)
const Z_MEM_ERROR     = Cint(-4)
const Z_BUF_ERROR     = Cint(-5)
const Z_VERSION_ERROR = Cint(-6)

# Compression levels
const Z_NO_COMPRESSION      = Cint(0)
const Z_BEST_SPEED          = Cint(1)
const Z_BEST_COMPRESSION    = Cint(9)
const Z_DEFAULT_COMPRESSION = Cint(-1)

# Compression strategy
const Z_FILTERED         = Cint(1)
const Z_HUFFMAN_ONLY     = Cint(2)
const Z_RLE              = Cint(3)
const Z_FIXED            = Cint(4)
const Z_DEFAULT_STRATEGY = Cint(0)

# Possible values of the data_type field
const Z_BINARY  = Cint(0)
const Z_TEXT    = Cint(1)
const Z_ASCII   = Z_TEXT
const Z_UNKNOWN = Cint(2)

# The deflate compression method
const Z_DEFLATED = Cint(8)

# For initializing zalloc, zfree, opaque
# const Z_NULL = C_NULL

# ZStream
# -------

type ZStream
    next_in::Ptr{UInt8}
    avail_in::Cuint
    total_in::Culong

    next_out::Ptr{UInt8}
    avail_out::Cuint
    total_out::Culong

    msg::Ptr{UInt8}
    state::Ptr{Void}

    zalloc::Ptr{Void}
    zfree::Ptr{Void}
    opaque::Ptr{Void}

    data_type::Cint
    adler::Culong
    reserved::Culong

    function ZStream()
        zstream = new()
        zstream.next_in   = C_NULL
        zstream.avail_in  = 0
        zstream.total_in  = 0
        zstream.next_out  = C_NULL
        zstream.avail_out = 0
        zstream.total_out = 0
        zstream.msg       = C_NULL
        zstream.state     = C_NULL
        zstream.zalloc    = C_NULL
        zstream.zfree     = C_NULL
        zstream.opaque    = C_NULL
        zstream.data_type = 0
        zstream.adler     = 0
        zstream.reserved  = 0
        return zstream
    end
end

# Functions
# ---------

function version()
    return unsafe_string(ccall((:zlibVersion, zlib), Ptr{UInt8}, ()))
end

const zlib_version = version()

function init_inflate!(zstream::ZStream, windowbits::Integer)
    return ccall(
        (:inflateInit2_, zlib),
        Cint,
        (Ref{ZStream}, Cint, Cstring, Cint),
        zstream, windowbits, zlib_version, sizeof(ZStream))
end

function init_inflate!(zstream::ZStream; windowbits::Integer=15 + 32)
    return init_inflate!(zstream, windowbits)
end

function reset_inflate!(zstream::ZStream)
    return ccall((:inflateReset, zlib), Cint, (Ref{ZStream},), zstream)
end

function end_inflate!(zstream::ZStream)
    return ccall((:inflateEnd, zlib), Cint, (Ref{ZStream},), zstream)
end

function inflate!(zstream::ZStream, flush::Integer)
    return ccall((:inflate, zlib), Cint, (Ref{ZStream}, Cint), zstream, flush)
end

function init_deflate!(zstream::ZStream,
                       level::Integer,
                       method::Integer,
                       windowbits::Integer,
                       memlevel::Integer,
                       strategy::Integer)
    return ccall(
        (:deflateInit2_, zlib),
        Cint,
        (Ref{ZStream}, Cint, Cint, Cint, Cint, Cint, Cstring, Cint),
        zstream, level, method, windowbits, memlevel, strategy,
        zlib_version, sizeof(ZStream))
end

function init_deflate!(zstream::ZStream;
                       level::Integer=Z_DEFAULT_COMPRESSION,
                       method::Integer=Z_DEFLATED,
                       windowbits::Integer=15,
                       memlevel::Integer=8,
                       strategy::Integer=Z_DEFAULT_STRATEGY)
    return init_deflate!(zstream, level, method, windowbits, memlevel, strategy)
end

function reset_deflate!(zstream::ZStream)
    return ccall((:deflateReset, zlib), Cint, (Ref{ZStream},), zstream)
end

function end_deflate!(zstream::ZStream)
    return ccall((:deflateEnd, zlib), Cint, (Ref{ZStream},), zstream)
end

function deflate!(zstream::ZStream, flush::Integer)
    return ccall((:deflate, zlib), Cint, (Ref{ZStream}, Cint), zstream, flush)
end

function crc32(crc::UInt32, data::Ptr{UInt8}, n::Int)
    return ccall((:crc32, zlib), Culong, (Culong, Ptr{Cchar}, Cuint), crc, data, n) % UInt32
end

function crc32()
    return crc32(UInt32(0), convert(Ptr{UInt8}, C_NULL), 0)
end

function adler32(adler::UInt32, data::Ptr{UInt8}, n::Int)
    return ccall((:adler32, zlib), Culong, (Culong, Ptr{Cchar}, Cuint), adler, data, n) % UInt32
end

function adler32()
    return adler32(UInt32(0), convert(Ptr{UInt8}, C_NULL), 0)
end

# Utils
# -----

# Convert the return code to a corresponding string.
function code2str(code::Cint)
    if code == Z_OK
        return "Z_OK"
    elseif code == Z_STREAM_END
        return "Z_STREAM_END"
    elseif code == Z_NEED_DICT
        return "Z_NEED_DICT"
    elseif code == Z_ERRNO
        return "Z_ERRNO"
    elseif code == Z_STREAM_ERROR
        return "Z_STREAM_ERROR"
    elseif code == Z_DATA_ERROR
        return "Z_DATA_ERROR"
    elseif code == Z_MEM_ERROR
        return "Z_MEM_ERROR"
    elseif code == Z_BUF_ERROR
        return "Z_BUF_ERROR"
    elseif code == Z_VERSION_ERROR
        return "Z_VERSION_ERROR"
    end
    throw(AssertionError("unknown return code: $(code)"))
end

# throw error exception based on ZStream and the return code
function zerror(zstream::ZStream, code::Cint)
    @assert code < 0
    if zstream.msg == C_NULL
        zerror(code)
    else
        error("zlib error: ", unsafe_string(zstream.msg), " (", code2str(code), ")")
    end
end

# throw error exception based on the return code
function zerror(code::Cint)
    @assert code < 0
    error("zlib error: ", code2str(code))
end

# Reset the error message field (.msg), call a zlib function (ex), and throw an
# exception if the return code is not Z_OK.
macro zcheck(ex)
    @assert ex.head == :call
    zstream = esc(ex.args[2])
    quote
        $(zstream).msg = C_NULL
        ret = $(esc(ex))
        if ret != Z_OK
            zerror($(zstream), ret)
        end
        ret
    end
end

function init_inflate_zstream(gzip::Bool)
    zstream = ZStream()
    @zcheck init_inflate!(zstream, gzip ? 32 + 15 : 15)
    return zstream
end

function init_deflate_zstream(gzip::Bool, level::Integer, mem_level::Integer, strategy)
    if !(0 <= level <= 9 || level == Z_DEFAULT_COMPRESSION)
        throw(ArgumentError("invalid zlib compression level"))
    end

    if !(1 <= mem_level <= 9)
        throw(ArgumentError("invalid zlib memory level"))
    end

    if strategy != Z_DEFAULT_STRATEGY &&
        strategy != Z_FILTERED &&
        strategy != Z_HUFFMAN_ONLY &&
        strategy != Z_RLE &&
        strategy != Z_FIXED
        throw(ArgumentError("invalid zlib strategy"))
    end

    zstream = ZStream()
    window_bits = gzip ? 16 + 15 : 15
    @zcheck init_deflate!(zstream, level, Z_DEFLATED, window_bits, mem_level, strategy)
    return zstream
end


text:
using DiffRules, SpecialFunctions, NaNMath
using Base.FastMath: fast_op, make_fastmath

# TODO use CSE here

for (M, f, arity) in DiffRules.diffrules()
  arity == 1 || continue
  dx = DiffRules.diffrule(M, f, :x)
  @eval begin
    @tangent $M.$f(x::Number) = $M.$f(x), ẋ -> ẋ * $dx
  end
end

for (M, f, arity) in DiffRules.diffrules()
  arity == 2 || continue
  da, db = DiffRules.diffrule(M, f, :a, :b)
  @eval begin
    @tangent $M.$f(a::Number, b::Number) = $M.$f(a, b), (ȧ, ḃ) -> ȧ*$da + ḃ*$db
  end
end

# Some specific overrides
# The DiffRules definitions are suboptimal due to repeated work in the tangent

@tangent function tanh(x)
  y = tanh(x)
  y, ẋ -> ẋ * (1 - y^2)
end

@tangent function exp(x)
  y = exp(x)
  y, ẋ -> ẋ * y
end

for f in [>, <, ==, ===, !=, in]
  @eval @tangent $f(a, b) = $f(a, b), (_, _) -> false
end

@tangent convert(T::Type{<:Real}, x::Real) = convert(T, x), (_, ẋ) -> convert(T, ẋ)

@tangent function Colon()(xs...)
  c = Colon()(xs...)
  c, (_...) -> zerolike(c)
end

zerolike(x::AbstractRange) =
  invoke(zerolike, Tuple{Any}, x)

DiffRules._abs_deriv(x::Complex) = x/abs(x)


text:
###################################################
## network/tools/edges.jl
## functions to study properties of edges/paths
###################################################

"""
	`findNearEdges` : find all edges that share a vertex with input e in network n, as a set of edges
"""
function findNearEdges(n::Network, e::Edge)
	orig = src(e)
	dest = dst(e)
	# find all edges that share a vertex with e
	nearEdges = [Edge(dest, newDest) for newDest in outneighbors(n.graph, dest)]
	append!(nearEdges, [Edge(orig, newDest) for newDest in outneighbors(n.graph, orig)])
	append!(nearEdges, [Edge(newOrig, orig) for newOrig in in_neighbors(n.graph, orig)])
	append!(nearEdges, [Edge(newOrig, dest) for newOrig in in_neighbors(n.graph, dest)])
	# remove duplicates and edge e itself
	nearEdges = Set(nearEdges)
	delete!(nearEdges, e)
	return nearEdges
end

"""
	`findNearEdgesSameType`	: find all edges that can be an out neighbor of input e in network n and are of the same type, as a set of edges
"""
function findNearEdgesSameType(n::Network, e::Edge)
	orig = src(e)
	dest = dst(e)
	# find all edges that share a vertex with e
	nearEdges = Set(out_edges(n.graph, dest))
	# remove duplicates and edge e itself (as well as the reverse of e)
	if reverse(e) in nearEdges
		delete!(nearEdges, reverse(e))
	end
	# remove edges not of same type as e
	for edge in nearEdges
		if n.roads[orig,dest].roadType != n.roads[src(edge),dst(edge)].roadType
			delete!(nearEdges, edge)
		end
	end
	return nearEdges
end

"""
	`clusterEdges` : group edges into rough clusters of neighbors
"""
function clusterEdges(n::Network, nNeighbors::Int=div(nRoads(n), 180))
	# define lists of edges and other inputs and outputs
	edgeList = collect(edges(n.graph))
	edgesAssigned = falses(length(edgeList))
	startingEdges = Set(Edge[])
	clusters = [Edge[] for i=1:nNeighbors]
	nodeMap = Dict{Int, Int}()
	# define starting points of clusters
	Random.seed!(2002)
	for i = 1:nNeighbors
		while true
			idx = rand(eachindex(edgeList))
			edge = edgeList[idx]
			if !(edge in startingEdges) && !(src(edge) in keys(nodeMap)) && !(dst(edge) in keys(nodeMap))
				push!(startingEdges, edge)
				nodeMap[src(edge)] = i
				nodeMap[dst(edge)] = i
				edgesAssigned[idx] = true
				clusters[i] = [edge]
				break
			end
		end
	end
	# as long as edges are not assigned, check if endpoints is already part of a neighborhood
	# if one is, assign edge to that neighborhood (arbitrarily giving priority to src)
	# if they aren't, skip edge and come back to it later
	while !all(edgesAssigned)
		for idx in shuffle(collect(eachindex(edgeList)))
			edge = edgeList[idx]
			if !edgesAssigned[idx]
				if src(edge) in keys(nodeMap)
					push!(clusters[nodeMap[src(edge)]], edge)
					edgesAssigned[idx] = true
					if !(dst(edge) in keys(nodeMap))
						nodeMap[dst(edge)] = nodeMap[src(edge)]
					end
				elseif dst(edge) in keys(nodeMap)
					push!(clusters[nodeMap[dst(edge)]], edge)
					edgesAssigned[idx] = true
					if !(src(edge) in keys(nodeMap))
						nodeMap[src(edge)] = nodeMap[dst(edge)]
					end
				end
			end
		end
	end
	return clusters
end

"""
	`flatten` : recursively concatenate all vectors in a vector of vectors
	Useful helper function for neighborhood continuity
"""
flatten(a::Array) = any(x -> isa(x, Array), a) ? flatten(vcat(map(flatten,a)...)) : a

"""
	`findEdges`	: find all edges that begin and end in the provided list of nodes
"""
function findEdges(n::Network, nodeList::Vector{Int})
	nodeSet = Set(nodeList)
	edges = Set(Edge[])
	for node in nodeList
		for neighbor in outneighbors(n.graph, node)
			if neighbor in nodeSet && n.roads[node, neighbor].roadType > 1
				push!(edges, Edge(node, neighbor))
			end
		end
	end
	return edges
end

"""
	`getFullPathEdges`	: given network timings and origin and destination, get path edges
	Wrapper for RoutingNetworks.getPathEdges that deals with particular case of road projection
	Args:
		t 		: network trip
		timings : current network timings
	Returns:
		Path for network trip as (Edge, Float) dictionary
"""
function getFullPathEdges(t::NetworkTrip, timings::NetworkTimings)
	sp = Dict(edge => 1. for edge in getPathEdges(timings, t.orig[2], t.dest[1]))
    if t.roadProj
        sp[Edge(t.orig[1], t.orig[2])] = t.orig[3]
        sp[Edge(t.dest[1], t.dest[2])] = t.dest[3]
    end
    return sp
end

"""
	`toPathEdgesDict`: given a trip and path (list of nodes), convert to (edge, weight) dictionary
"""
function toPathEdgesDict(t::NetworkTrip, path::Vector{Int})
	pathEdges = Dict(Edge(path[i], path[i+1]) => 1. for i = 1:(length(path)-1))
	if t.roadProj
        pathEdges[Edge(t.orig[1], t.orig[2])] = t.orig[3]
        pathEdges[Edge(t.dest[1], t.dest[2])] = t.dest[3]
    end
	return pathEdges
end


text:
##########################################
# Master file for running all test cases #
##########################################
using Turing; turnprogress(false)
using Test

@debug("[runtests.jl] runtests.jl loaded")

path = dirname(@__FILE__)
cd(path); include("utility.jl")
@debug("[runtests.jl] utility.jl loaded")
@debug("[runtests.jl] testing starts")

if get(ENV, "TRAVIS", "false") == "true"
    # If Travis is testing, separate the tests.
    numerical_tests = [joinpath("hmc.jl", "matrix_support.jl"),
                       joinpath("mh.jl", "mh_cons.jl")]

    if ENV["STAGE"] == "test"
        runtests(exclude = numerical_tests)
    elseif ENV["STAGE"] == "numerical"
        runtests(specific_tests = numerical_tests)
    else
        @warn "Unknown Travis stage, currently set to: $(ENV["STAGE"])"
    end
else
    # Otherwise, test everything.
    runtests()
end

@debug("[runtests.jl] all tests finished")


text:
STATIC_PROPERTIES = [:n, :num_steps, :dt, :seed, :randomize, :nmac]

function initialize(log::Dict{String,Any}, sim::Simulator)
    sim.logging || return

    for p in STATIC_PROPERTIES
        log["$p"] = getproperty(sim, p)
    end

    log["t"] = [sim.step * sim.dt]

    for i in 1:sim.n
    	log["ac_$i"] = [vec(sim.acs[i])]
    	log["cas_$i"] = [vec(sim.cass[i])]
    end

    for (k, v) in sim.metrics.d
        log[k] = [v]
    end
end

function update(log::Dict{String,Any}, sim::Simulator)
    sim.logging || return
    push!(log["t"], sim.step * sim.dt)

    for i in 1:sim.n
    	push!(log["ac_$i"], vec(sim.acs[i]))
    	push!(log["cas_$i"], vec(sim.cass[i]))
    end

    for (k, v) in sim.metrics.d
        push!(log[k], v)
    end
end


text:
module Simple

export simple
include("dict.jl")

simple(a::Any) = string(a)

end # module


text:
using Test, YaoBase, YaoBlocks, YaoArrayRegister

@testset "test chain" begin
    include("chain.jl")
end

@testset "test kron" begin
    include("kron.jl")
end

@testset "test control" begin
    include("control.jl")
end

@testset "test put" begin
    include("put.jl")
end

@testset "test repeated" begin
    include("repeated.jl")
end

@testset "test subroutine" begin
    include("subroutine.jl")
end

@testset "test tag" begin
    include("tag.jl")
    include("cache.jl")
end

@testset "test unitary channel" begin
    include("unitary_channel.jl")
end

@testset "test single block chsubblocks" begin
    @test chsubblocks(chain(X), Y) == chain(Y)
    @test chsubblocks(kron(X), Y) == kron(Y)
    @test chsubblocks(*(X), Y) == *(Y)
    @test chsubblocks(+(X), Y) == +(Y)
end

@testset "extension" begin
    # check extension fallback errors
    struct MockedQFT{D} <: CompositeBlock{D} end
    @test_throws MethodError ishermitian(MockedQFT{2}())
    @test_throws MethodError isunitary(MockedQFT{2}())
    @test Val(im)*X == im*X
end