text:
# Some discussion on this:
# https://groups.google.com/forum/#!msg/julia-users/YP31LM3Qto0/ET-XjN-vQuAJ
#
# All the model parameters
#
# TODO: improve macro hygiene.

__precompile__()

module Parameters
import Base: @__doc__
import DataStructures: OrderedDict
using Compat

export @with_kw, type2dict, reconstruct, @unpack, @pack

## Parser helpers
#################

# To iterate over code blocks dropping the line-number bits:
immutable Lines
    block::Expr
end
Base.start(lns::Lines) = 1
function Base.next(lns::Lines, nr)
    for i=nr:length(lns.block.args)
        if isa(lns.block.args[i], LineNumberNode)
            continue
        end
        if isa(lns.block.args[i], Symbol) || !(lns.block.args[i].head==:line)
            return lns.block.args[i], i+1
        end
    end
    return -1
end
function Base.done(lns::Lines, nr)
    if next(lns::Lines, nr)==-1
        true
    else
        false
    end
end

# Transforms :(a::b) -> :a
decolon2(a::Expr) = (@assert a.head==:(::);  a.args[1])
decolon2(a::Symbol) = a

# Returns the name of the type as Symbol
function typename(typedef::Expr)
    if isa(typedef.args[2], Symbol)
        return typedef.args[2]
    elseif isa(typedef.args[2].args[1], Symbol)
        return typedef.args[2].args[1]
    elseif isa(typedef.args[2].args[1].args[1], Symbol)
        return typedef.args[2].args[1].args[1]
    else
        error("Could not parse type-head from: $typedef")
    end
end

# Transforms:  Expr(:<:, :A, :B) -> :A
stripsubtypes(s::Symbol) = s
function stripsubtypes(e::Expr)
    e.args[1]
end
stripsubtypes(vec::Vector) = [stripsubtypes(v) for v in vec]

function check_inner_constructor(l)
    if length(l.args[1].args)==1
        error("No inner constructors with zero positional arguments allowed!")
    elseif (length(l.args[1].args)==2 #1<length(l.args[1].args)<=3
            && isa(l.args[1].args[2], Expr)
            && l.args[1].args[2].head==:parameters)
        error("No inner constructors with zero positional arguments plus keyword arguments allowed!")
    end
    nothing
end

## Exported helper functions
#####################
"""
Transforms a type-instance into a dictionary.

```
julia> type T
           a
           b
       end

julia> type2dict(T(4,5))
Dict{Symbol,Any} with 2 entries:
  :a => 4
  :b => 5
```
"""
function type2dict(dt)
    di = Dict{Symbol,Any}()
    for n in @compat fieldnames(dt)
        di[n] = getfield(dt, n)
    end
    di
end

"""
Make a new instance of a type with the same values as
the input type except for the fields given in the associative
second argument or as keywords.

```julia
type A; a; b end
a = A(3,4)
b = reconstruct(a, [(:b, 99)]) # ==A(3,99)
```
"""
function reconstruct{T}(pp::T, di)
    di = !isa(di, Associative) ? Dict(di) : di
    ns = @compat fieldnames(pp)
    args = Array(Any, length(ns))
    for (i,n) in enumerate(ns)
        args[i] = get(di, n, getfield(pp, n))
    end
    T(args...)
end
reconstruct{T}(pp::T; kws...) = reconstruct(pp, kws)

###########################
# Keyword constructors with @with_kw
##########################

# A type with fields (r,a) in variable aa becomes
# quote
#     r = aa.r
#     a = aa.a
# end
_unpack(binding, fields) = Expr(:block, [:($f = $binding.$f) for f in fields]...)
# Pack fields back into binding using reconstruct:
function _pack(binding, fields)
    kws = [Expr(:kw, f, f) for f in fields]
    :($binding = Main.Parameters.reconstruct($binding, $(kws...)) )
end

"""
This function is called by the `@with_kw` macro and does the syntax
transformation from:

```julia
@with_kw immutable MM{R}
    r::R = 1000.
    a::R
end
```

into

```julia
immutable MM{R}
    r::R
    a::R
    MM(r,a) = new(r,a)
    MM(;r=1000., a=error("no default for a")) = MM{R}(r,a) # inner kw, type-paras are required when calling
end
MM{R}(r::R,a::R) = MM{R}(r,a) # default outer positional constructor
MM(;r=1000,a=error("no default for a")) =  MM(r,a) # outer kw, so no type-paras are needed when calling
MM(m::MM; kws...) = reconstruct(mm,kws)
MM(m::MM, di::Union{Associative, Tuple{Symbol,Any}}) = reconstruct(mm, di)
macro unpack_MM(varname)
    esc(quote
    r = varname.r
    a = varname.a
    end)
end
macro pack_MM(varname)
    esc(quote
    varname = Main.Parameters.reconstruct(varname,r=r,a=a)
    end)
end
```
"""
function with_kw(typedef)
    if typedef.head!=:type
        error("only works on type-defs")
    end
    const err1str = "Field \'"
    const err2str = "\' has no default, supply it with keyword."

    inner_constructors = Any[]

    # parse a few things
    tn = typename(typedef) # the name of the type
    # Returns M{...} (removes any supertypes)
    if isa(typedef.args[2], Symbol)
        typparas = Any[]
    elseif typedef.args[2].head==:<:
        if isa(typedef.args[2].args[1],Symbol)
            typparas = Any[]
        else
            typparas = typedef.args[2].args[1].args[2:end]
        end
    else
        typparas = typedef.args[2].args[2:end]
    end

    # error on types without fields
    lns = Lines(typedef.args[3])
    if done(lns, start(lns))
        error("@with_kw only supported for types which have at least one field.")
    end
    # default type @deftype
    l, i = next(lns, start(lns))
    if isa(l, Expr) && l.head==:macrocall && l.args[1]==Symbol("@deftype")
        has_deftyp = true
        if length(l.args) != 2
            error("Malformed `@deftype` line")
        end
        deftyp = l.args[2]
        if done(lns, i)
            error("@with_kw only supported for types which have at least one field.")
        end
    else
        has_deftyp = false
    end

    # the vars for the unpack macro
    unpack_vars = Any[]
    # the type def
    fielddefs = quote end # holds r::R etc
    fielddefs.args = Any[] # in julia 0.5 this is [:( # /home/mauro/.julia/v0.5/Parameters/src/Parameters.jl, line 228:)]
    kws = OrderedDict{Any, Any}()
    # assertions in the body
    asserts = Any[]
    for (i,l) in enumerate(lns) # loop over body of typedef
        if i==1 && has_deftyp
            continue
        end
        if isa(l, Symbol)  # no default value and no type annotation
            if has_deftyp
                push!(fielddefs.args, :($l::$deftyp))
            else
                push!(fielddefs.args, l)
            end
            sym = l
            syms = string(sym)
            kws[sym] = :(error($err1str * $syms * $err2str))
            # unwrap-macro
            push!(unpack_vars, sym)
        elseif l.head==:(=)  # default value and with or without type annotation
            if isa(l.args[1], Expr) && l.args[1].head==:call # inner constructors
                check_inner_constructor(l)
                push!(inner_constructors, l)
            else
                fld = l.args[1]
                if isa(fld, Symbol) && has_deftyp # no type annotation
                    fld = :($fld::$deftyp)
                end
                push!(fielddefs.args, fld)
                kws[decolon2(fld)] = l.args[2]
                # unwrap-macro
                push!(unpack_vars, decolon2(fld))
            end
        elseif l.head==:macrocall  && l.args[1]==Symbol("@assert")
            # store all asserts
            push!(asserts, l)
        elseif l.head==:function # inner constructor
            check_inner_constructor(l)
            push!(inner_constructors, l)
        else # no default value but with type annotation
            push!(fielddefs.args, l)
            sym = decolon2(l.args[1])
            syms = string(sym)
            kws[sym] = :(error($err1str *$syms * $err2str))
            # unwrap-macro
            push!(unpack_vars, l.args[1])
        end
    end
    # The type definition without inner constructors:
    typ = Expr(:type, deepcopy(typedef.args[1:2])..., deepcopy(fielddefs))

    # Inner keyword constructor.  Note that this calls the positional
    # constructor under the hood and not `new`.  That way a user can
    # provide a special positional constructor (say enforcing
    # invariants) which also gets used with the keywords.
    args = Any[]
    kwargs = Expr(:parameters)
    for (k,w) in kws
        push!(args, k)
        push!(kwargs.args, Expr(:kw,k,w))
    end
    if length(typparas)>0
        innerc = :($tn($kwargs) = $tn{$(stripsubtypes(typparas)...)}($(args...)) )
    else
        innerc = :($tn($kwargs) = $tn($(args...)) )
    end
    push!(typ.args[3].args, innerc)

    # Inner positional constructor: only make it if no inner
    # constructors are user-defined.  If one or several are defined,
    # assume that one has the standard positional signature.
    if length(inner_constructors)==0
        innerc2 = :($tn($(args...)) = new($(args...)))
        prepend!(innerc2.args[2].args, asserts)
        push!(typ.args[3].args, innerc2)
    else
        if length(asserts)>0
            error("Assertions are only allowed in type-definitions with no inner constructors.")
        end
        append!(typ.args[3].args, inner_constructors)
    end

    # Outer positional constructor which does not need explicit
    # type-parameters when called.  Only make this constructor if
    #  (1) type parameters are used at all
    #  (2) all type parameters are used in the fields (otherwise get a
    #      "method is not callable" warning!)
    #       See also https://github.com/JuliaLang/julia/issues/17186
    if typparas!=Any[] # condition (1)
        outer_positional = :(  $tn{$(typparas...)}($(fielddefs.args...))
                             = $tn{$(stripsubtypes(typparas)...)}($(args...)))
        # Check condition (2)
        used_paras = Any[]
        for f in fielddefs.args
            if !isa(f,Symbol) && f.head==:(::)
                push!(used_paras, f.args[2])
            end
        end
        if !issubset(stripsubtypes(typparas), used_paras)
            outer_positional = :()
        end
    else
        outer_positional = :()
    end

    # Outer keyword constructor, useful to infer the type parameter
    # automatically.  This calls the outer positional constructor.
    # only create if type parameters are used.
    if typparas==Any[]
        outer_kw=:()
    else
        outer_kw = :($tn($kwargs) = $tn($(args...)) )
    end

    ## outer copy constructor
    ###
    outer_copy = quote
        $tn(pp::$tn; kws... ) = reconstruct(pp, kws)
        # $tn(pp::$tn, di::Union(Associative,Vararg{Tuple{Symbol,Any}}) ) = reconstruct(pp, di) # see issue https://github.com/JuliaLang/julia/issues/11537
        # $tn(pp::$tn, di::Union(Associative, Tuple{Vararg{Tuple{Symbol, Any}}}) ) = reconstruct(pp, di) # see issue https://github.com/JuliaLang/julia/issues/11537
        $tn(pp::$tn, di::Associative ) = reconstruct(pp, di)
        $tn(pp::$tn, di::Vararg{Tuple{Symbol,Any}} ) = reconstruct(pp, di)
    end

    # (un)pack macro from https://groups.google.com/d/msg/julia-users/IQS2mT1ITwU/hDtlV7K1elsJ
    unpack_name = Symbol("unpack_"*string(tn))
    pack_name = Symbol("pack_"*string(tn))
    # Finish up
    quote
        Main.Parameters.@__doc__ $typ # use Main.Parameters.@__doc__ for 0.3 compatibility
        $outer_positional
        $outer_kw
        $outer_copy
        function Base.show(io::IO, p::$tn)
            println(io, string(typeof(p)))
            for (i, var) in enumerate($unpack_vars)
                print(io, "  " * string(var) * ": $(getfield(p,var))")
                i == length($unpack_vars) || print(io, "\n")
            end
        end
        macro $unpack_name(ex)
            esc(Main.Parameters._unpack(ex, $unpack_vars))
        end
        macro $pack_name(ex)
            esc(Main.Parameters._pack(ex, $unpack_vars))
        end
        $tn
    end
end

"""
Macro which allows default values for field types and a few other features.

Basic usage:

```julia
@with_kw immutable MM{R}
    r::R = 1000.
    a::Int = 4
end
```

For more details see manual.
"""
macro with_kw(typedef)
    return esc(with_kw(typedef))
end

###########################
# Packing and unpacking @unpack, @pack
##########################
# Below code slightly adapted from Simon Danisch's GLVisualize via PR
# https://github.com/mauro3/Parameters.jl/pull/13

"""
This function is invoked to unpack one entity of some DataType and has
signature:

`unpack(dt::Any, ::Val{field}) -> value of field`

Note that this means the only symbols or immutable field-descriptors
are allowed, as they are used as type parameter in `Val`.

Two definitions are included in the package to unpack a composite type
or a dictionary:
```
@inline unpack{f}(x, ::Val{f}) = getfield(x, f)
@inline unpack{k}(x::Associative{Symbol}, ::Val{k}) = x[k]
```

More methods can be added to allow for specialized unpacking of other datatypes.

See also `pack!`.
"""
function unpack end
@inline unpack{f}(x, ::Val{f}) = getfield(x, f)
@inline unpack{k}(x::Associative{Symbol}, ::Val{k}) = x[k]

"""
This function is invoked to pack one entity into some DataType and has
signature:

`pack!(dt::Any, ::Val{field}, value) -> value`

Note that this means the only symbols or immutable field-descriptors
are allowed, as they are used as type parameter in `Val`.

Two definitions are included in the package to pack into a composite
type or into a dictionary:

```
@inline pack!{f}(x, ::Val{f}, val) = setfield!(x, f, val)
@inline pack!{k}(x::Associative{Symbol}, ::Val{k}, val) = x[k]=val
```

More methods can be added to allow for specialized packing of other
datatypes.

See also `unpack`.
"""
function pack! end
@inline pack!{f}(x, ::Val{f}, val) = setfield!(x, f, val)
@inline pack!{k}(x::Associative{Symbol}, ::Val{k}, val) = x[k]=val

"""
Unpacks fields/keys from a composite type or a `Dict{Symbol}` into variables
```julia_skip
@unpack a, b, c = dict_or_typeinstance
```

Example with dict:
```julia
d = Dict{Symbol,Any}(:a=>5.0,:b=>2,:c=>"Hi!")
@unpack a, c = d
a == 5.0 #true
c == "Hi!" #true
```

Example with type:
```julia
type A; a; b; c; end
d = A(4,7.0,"Hi")
@unpack a, c = d
a == 4 #true
c == "Hi!" #true
```
"""
macro unpack(args)
    args.head!=:(=) && error("Expression needs to be of form `a, b = c`")
    items, suitecase = args.args
    items = isa(items, Symbol) ? [items] : items.args
    suitecase_instance = gensym()
    kd = [:( $key = Parameters.unpack($suitecase_instance, Val{$(Expr(:quote, key))}()) )for key in items]
    kdblock = Expr(:block, kd...)
    expr = quote
        $suitecase_instance = $suitecase # handle if suitecase is not a variable but an expression
        $kdblock
    end
    esc(expr)
end

"""
Packs variables into a composite type or a `Dict{Symbol}`
```julia_skip
@pack dict_or_typeinstance = a, b, c
```

Example with dict:
```julia
a = 5.0
c = "Hi!"
d = Dict{Symbol,Any}()
@pack d = a, c
d # Dict{Symbol,Any}(:a=>5.0,:c=>"Hi!")
```

Example with type:
```julia
a = 99
c = "HaHa"
type A; a; b; c; end
d = A(4,7.0,"Hi")
@pack d = a, c
d.a == 99 #true
d.c == "HaHa" #true
```
"""
macro pack(args)
    args.head!=:(=) && error("Expression needs to be of form `a = b,c`")
    suitecase, items = args.args
    items = isa(items, Symbol) ? [items] : items.args
    suitecase_instance = gensym()
    kd = [:( Parameters.pack!($suitecase_instance, Val{$(Expr(:quote, key))}(), $key) ) for key in items]
    kdblock = Expr(:block, kd...)
    expr = quote
        $suitecase_instance = $suitecase # handle if suitecase is not a variable but an expression
        $kdblock
    end
    esc(expr)
end

# TODO: maybe add @pack_new for packing into a new instance.  Could be
# used with immutables also.

end # module


text:
export StateSpaceModel, LinearStateSpaceModel, TimeVariantLinearStateSpaceModel, GenericStateSpaceModel, ProposalStateSpaceModel, BootstrapStateSpaceModel, initial_distribution, transition_distribution, observation_distribution, proposal_distribution, weight_function,
    initial_mean,
    initial_covariance,
    state_transition_matrix,
    state_covariance,
    state_input_matrix,
    observation_matrix,
    observation_covariance,
    observation_input_matrix,
    state_noise_transformation_matrix

abstract type StateSpaceModel end

# Default definitions of the modeling interface functions
initial_rand(m::StateSpaceModel,θ,N) = rand(initial_distribution(m,θ),N)
initial_logpdf(m::StateSpaceModel,x,θ) = logpdf(initial_distribution(m,θ),x)

initial_gradient!(dg,m::StateSpaceModel,x,θ) = ForwardDiff.gradient!(dg,θ->initial_logpdf(m,x,θ),θ)

transition_rand(m::StateSpaceModel,x,u,t,θ) = rand(transition_distribution(m,x,u,t,θ))
transition_logpdf(m::StateSpaceModel,xn,x,u,t,θ) = logpdf(transition_distribution(m,x,u,t,θ),xn)

observation_rand(m::StateSpaceModel,x,u,t,θ) = rand(observation_distribution(m,x,u,t,θ))
observation_logpdf(m::StateSpaceModel,y,x,u,t,θ) = logpdf(observation_distribution(m,x,u,t,θ),y)

ssm_gradient!(dg,m::StateSpaceModel,xn,x,y,u,t,θ) = ForwardDiff.gradient!(dg,θ->transition_logpdf(m,xn,x,u,t,θ) + observation_logpdf(m,y,xn,u,t,θ),θ)

proposal_rand(m::StateSpaceModel,x,y,u,t,θ) = rand(proposal_distribution(m,x,y,u,t,θ))
proposal_rand!(xn,m::StateSpaceModel,x,y,u,t,θ) = rand!(proposal_distribution(m,x,y,u,t,θ),xn)
proposal_logpdf(m::StateSpaceModel,xn,x,y,u,t,θ) = logpdf(proposal_distribution(m,x,y,u,t,θ),xn)

"""
A linear, time-invariant state space model

A degenerate process noise is allowed by
using Θ ≠ I, though correlated process
and measurement noise are not allowed.

To interface your own state space model
with the filtering/smoothing and parameter
estimation code, provide a function

```
buildFunction(θ)::LinearStateSpaceModel
```

where θ is a vector of parameters.
"""
struct LinearStateSpaceModel <: StateSpaceModel
    F
    Q
    G
    H
    R
    Γ
    Θ
    μ
    Σ0
end

function initial_distribution(m::LinearStateSpaceModel,θ0)
    μ0 = m.μ(θ0)
    Σ0 = m.Σ0(θ0)

    MvNormal(μ0,Σ0)
end

function transition_distribution(m::LinearStateSpaceModel,x,u,t,θ0)
    F = m.F(θ0)
    G = m.G(θ0)
    Q = m.Q(θ0)
    Θ = m.Θ(θ0)

    MvNormal(F*x+G*u,Θ*Q*Θ')
end

function observation_distribution(m::LinearStateSpaceModel,x,u,t,θ0)
    H = m.H(θ0)
    Γ = m.Γ(θ0)
    R = m.R(θ0)

    MvNormal(H*x + Γ*u,R)
end

function proposal_distribution(m::LinearStateSpaceModel,x,y,u,t,θ)
    transition_distribution(m,x,u,t,θ)
end

function weight_function(m::LinearStateSpaceModel,x,y,u,t,θ)
    0.0
end

### Linear model interface

initial_mean(m::LinearStateSpaceModel,θ) = m.μ(θ)
initial_covariance(m::LinearStateSpaceModel,θ) = m.Σ0(θ)
state_transition_matrix(m::LinearStateSpaceModel,θ) = m.F(θ)
state_covariance(m::LinearStateSpaceModel,θ) = m.Q(θ)
state_input_matrix(m::LinearStateSpaceModel,θ) = m.G(θ)
observation_matrix(m::LinearStateSpaceModel,θ) = m.H(θ)
observation_covariance(m::LinearStateSpaceModel,θ) = m.R(θ)
observation_input_matrix(m::LinearStateSpaceModel,θ) = m.Γ(θ)
state_noise_transformation_matrix(m::LinearStateSpaceModel,θ) = m.Θ(θ)

struct LinearStateSpaceGradient
    dF
    dQ
    dH
    dR
    dΘ
    dμ
    dΣ0
end

function simulate(m::LinearStateSpaceModel,θ,N,u=fill([0.0],N))
    μ0 = initial_distribution(m,θ)
    X = fill(rand(μ0),N+1)
    Y = fill(rand(observation_distribution(m,X[1],u[1],0,θ)),N)
    for t in 1:N
        X[t+1] = rand(transition_distribution(m,X[t],u[t],t,θ))
        Y[t] = rand(observation_distribution(m,X[t+1],u[t],t,θ))
    end
    X,Y
end

#=
function simulate(m::LinearStateSpaceModel,N,u=fill([0.0],N))
    d = size(m.μ,1)
    n = size(m.H*m.μ,1)
    X = rand(MvNormal(m.μ,m.Σ0),N+1)
    Y = rand(MvNormal(m.H*m.μ,m.R),N+1)
    for i in 2:N+1
        X[:,i] = m.F*X[:,i-1] + m.G*u[i-1] + rand(MvNormal(zeros(d),m.Q))
        Y[:,i] = m.H*X[:,i] .+ m.Γ*u[i-1] .+ rand(MvNormal(zeros(n),m.R))
    end
    X,Y
end
=#

"""
A linear, time-variant state space model

Each of the model components (with the
exception of the initial state and
covariance estimates) is a vector of
matrices. F[i] gives the matrix F at
the `i`th time step.

A degenerate process noise is allowed by
using Θ ≠ I, though correlated process
and measurement noise are not allowed.

To interface your own state space model
with the filtering/smoothing and parameter
estimation code, provide a function

```
buildFunction(θ)::LinearStateSpaceModel
```

where θ is a vector of parameters.
"""
struct TimeVariantLinearStateSpaceModel <: StateSpaceModel
    F
    Q
    G
    H
    R
    Γ
    Θ
    μ
    Σ0
end

function simulate(m::TimeVariantLinearStateSpaceModel,N,u=fill([0.0],N))
    d = size(m.μ,1)
    n = size(m.H[1]*m.μ,1)
    X = rand(MvNormal(m.μ,m.Σ0),N+1)
    Y = rand(MvNormal(m.H[1]*m.μ,m.R[1]),N+1)
    for i in 2:N+1
        X[:,i] = m.F[i-1]*X[:,i-1] + m.G[i-1]*u[i-1] + rand(MvNormal(zeros(d),m.Q[i-1]))
        Y[:,i] = m.H[i-1]*X[:,i] + m.Γ[i-1]*u[i-1] + rand(MvNormal(zeros(n),m.R[i-1]))
    end
    X,Y
end

"""
A generic state space model for use with particle filters

Three functions are needed:

- `μ0(θ=θ0) :: Distribution`
- `f(x,θ=θ0) :: Distribution`
- `g(x,θ=θ0) :: Distribution`

`μ0(θ)` is the initial distribution over hidden states
i.e. `x = rand(μ0(θ))` is a sample from that initial distribution

`f(x,θ)` is the transition probability distribution to move from
x to a new state. `xn = rand(f(x,θ))` propagates a particle forward.
`pdf(f(x,θ),xn)` evaluates the density of `f` at the new hidden state.

`g(x,θ)` is the observation probability distribution. `pdf(g(x,θ),y)`
evaluates the density of `g` at the observation `y`

If you don't want to make use of the maximum likelihood estimation tools, you can provide one-argument versions of the functions. You would then use the `particle_filter` method without the parameter.
"""
struct GenericStateSpaceModel <: StateSpaceModel
    μ0
    f
    g
end

"""
A state space model with a proposal distribution for auxiliary filtering

`μ0`, `f` and `g` resemble those for a GenericStateSpaceModel

`q(x,y,t,θ)` creates a proposal distribution at current state x,
next observation y, time t and parameter value θ.

`ξ(x,y,t,θ)` gives the value of the weight multiplier so that particle (x,w)
is resampled with probability `w*ξ(x,y,t,θ)`.

"""
struct ProposalStateSpaceModel <: StateSpaceModel
    μ0
    f
    g
    q
    ξ
    dp
    ProposalStateSpaceModel(μ0,f,g,q,ξ,dp=DP(f,g)) = new(μ0,f,g,q,ξ,dp)
end

initial_distribution(m::ProposalStateSpaceModel,θ) = m.μ0(θ)

transition_distribution(m::ProposalStateSpaceModel,x,u,t,θ) = m.f(x,u,t,θ)

observation_distribution(m::ProposalStateSpaceModel,x,u,t,θ) = m.g(x,u,t,θ)

proposal_distribution(m::ProposalStateSpaceModel,x,y,u,t,θ) = m.q(x,y,u,t,θ)

weight_function(m::ProposalStateSpaceModel,x,y,u,t,θ) = m.ξ(x,y,u,t,θ)

function DP(f,g)
    function dp(dh,x,xn,y,t,θp)
        ForwardDiff.gradient!(dh,θ->logpdf(g(xn,t,θ),y) + logpdf(f(x,t,θ),xn),θp)
    end
end

function BootstrapStateSpaceModel(μ0,f,g)
    q(x,y,u,t,θ) = f(x,u,t,θ)
    ξ(x,y,u,t,θ) = 0.0 # We do this in logarithms now
    ProposalStateSpaceModel(μ0,f,g,q,ξ)
end

function simulate(m::ProposalStateSpaceModel,θ,N,u=fill([0.0],N))
    X = fill(rand(m.μ0(θ)),N+1)
    Y = fill(rand(m.g(X[1],u[1],1,θ)),N)
    for t in 1:N
        X[t+1] =  rand(m.f(X[t],u[t],t,θ))
        Y[t] = rand(m.g(X[t+1],u[t],t,θ))
    end
    X,Y
end

(m::ProposalStateSpaceModel)(θ,N) = simulate(m,θ,N)

function simulate(m::StateSpaceModel,θ,N,u=fill([0.0],N))
    X = fill(initial_rand(m,θ),N+1)
    Y = fill(observation_rand(m,X[1],u[1],0,θ),N)
    for t in 1:N
        X[t+1] =  transition_rand(m,X[t],u[t],t,θ)
        Y[t] = observation_rand(m,X[t+1],u[t],t,θ)
    end
    X,Y
end

ext: jl


text:
# analytical
@testset "Chaos" begin
    Random.seed!(SEED)
    @testset "Array" begin
        x = randn(252)
        # helpers
        a, b = Indicators.divide(x)
        @test [a; b] == x
        x = randn(101)
        a, b = Indicators.divide(x)
        @test [a; b] == x
        # workhorses
        h = hurst(x, n=100)
        @test size(h) == size(x)
        rs = rsrange(x)
        @test size(rs) == size(x)
        x = randn(100)
    end
    @testset "Temporal" begin
        # chaos indicators
        x = TS(randn(N))
        tmp = hurst(x)
        @test size(tmp,1) == size(x,1)
        @test size(tmp,2) == 1
        tmp = rsrange(x)
        @test size(tmp,1) == size(x,1)
        @test size(tmp,2) == 1
    end
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

struct TrustedIdProvidersApi <: SwaggerApi
    client::Swagger.Client
end

"""

Creates or updates the specified trusted identity provider. During update, the trusted identity provider with the specified name will be replaced with this new provider
Param: resourceGroupName::String (required)
Param: accountName::String (required)
Param: trustedIdProviderName::String (required)
Param: parameters::TrustedIdProvider (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: TrustedIdProvider
"""
function trustedIdProvidersCreateOrUpdate(_api::TrustedIdProvidersApi, resourceGroupName::String, accountName::String, trustedIdProviderName::String, parameters, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "PUT", TrustedIdProvider, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.DataLakeStore/accounts/{accountName}/trustedIdProviders/{trustedIdProviderName}", ["azure_auth"], parameters)
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "accountName", accountName)  # type String
    Swagger.set_param(_ctx.path, "trustedIdProviderName", trustedIdProviderName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    Swagger.exec(_ctx)
end

"""

Deletes the specified trusted identity provider from the specified Data Lake Store account
Param: resourceGroupName::String (required)
Param: accountName::String (required)
Param: trustedIdProviderName::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: Nothing
"""
function trustedIdProvidersDelete(_api::TrustedIdProvidersApi, resourceGroupName::String, accountName::String, trustedIdProviderName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "DELETE", Nothing, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.DataLakeStore/accounts/{accountName}/trustedIdProviders/{trustedIdProviderName}", ["azure_auth"])
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "accountName", accountName)  # type String
    Swagger.set_param(_ctx.path, "trustedIdProviderName", trustedIdProviderName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    Swagger.exec(_ctx)
end

"""

Gets the specified Data Lake Store trusted identity provider.
Param: resourceGroupName::String (required)
Param: accountName::String (required)
Param: trustedIdProviderName::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: TrustedIdProvider
"""
function trustedIdProvidersGet(_api::TrustedIdProvidersApi, resourceGroupName::String, accountName::String, trustedIdProviderName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", TrustedIdProvider, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.DataLakeStore/accounts/{accountName}/trustedIdProviders/{trustedIdProviderName}", ["azure_auth"])
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "accountName", accountName)  # type String
    Swagger.set_param(_ctx.path, "trustedIdProviderName", trustedIdProviderName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    Swagger.exec(_ctx)
end

"""

Lists the Data Lake Store trusted identity providers within the specified Data Lake Store account.
Param: resourceGroupName::String (required)
Param: accountName::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: DataLakeStoreTrustedIdProviderListResult
"""
function trustedIdProvidersListByAccount(_api::TrustedIdProvidersApi, resourceGroupName::String, accountName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", DataLakeStoreTrustedIdProviderListResult, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.DataLakeStore/accounts/{accountName}/trustedIdProviders", ["azure_auth"])
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "accountName", accountName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    Swagger.exec(_ctx)
end

"""

Updates the specified trusted identity provider.
Param: resourceGroupName::String (required)
Param: accountName::String (required)
Param: trustedIdProviderName::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Param: parameters::UpdateTrustedIdProviderParameters
Return: TrustedIdProvider
"""
function trustedIdProvidersUpdate(_api::TrustedIdProvidersApi, resourceGroupName::String, accountName::String, trustedIdProviderName::String, api_version::String, subscriptionId::String; parameters=nothing, _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "PATCH", TrustedIdProvider, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.DataLakeStore/accounts/{accountName}/trustedIdProviders/{trustedIdProviderName}", ["azure_auth"], parameters)
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "accountName", accountName)  # type String
    Swagger.set_param(_ctx.path, "trustedIdProviderName", trustedIdProviderName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    Swagger.exec(_ctx)
end

export trustedIdProvidersCreateOrUpdate, trustedIdProvidersDelete, trustedIdProvidersGet, trustedIdProvidersListByAccount, trustedIdProvidersUpdate


text:
# # 介绍
#
# *这是一篇平和的关于并行化和 Julia GPU 编程的介绍*
#
# [Julia](https://julialang.org/)拥有对 GPU 编程最高等的支持：
# 你可以用高度抽象或得到细粒度的控制，而不需要离开你最爱的编程语言。
# 这篇指南的目的是帮助 Julia 用户们迈出他们 GPU 计算的第一步。
# 在这片指南中，你将会比较 CPU 和 GPU 实行同一个简单的计算，了解到一些会影响到你得到的性能的因素。
#
# 这篇指南部分地受到了一篇 Mark Harris 的博客的启发 [An Even Easier
# Introduction to CUDA](https://devblogs.nvidia.com/even-easier-introduction-cuda/)。
# 那是一篇使用编程语言 C++ 的更加容易的关于 CUDA 的介绍。
# 你不需要阅读那篇指南，因为这篇会从零开始讲起。

# ## CPU 上的一个简单例子

# 我们来研究一下这个例子，一个 CPU 上的简单计算。

N = 2^20
x = fill(1.0f0, N)  # a vector filled with 1.0 (Float32)
y = fill(2.0f0, N)  # a vector filled with 2.0

y .+= x             #  # 递增y的每个元素和其相应的元素x。

# 检查后发现我们得到了正确的答案
using Test
@test all(y .== 3.0f0)

# 从“Test Passed”显示后我们知道一切已经就绪。
# 我们使用了 Float32 数据类型，为切换到 GPU 计算做准备：
# GPU 在使用 Float32 类型时比用 Float64 更快（有些时候快得非常多）。

# 这个计算的一个显著特征是，y 的每个元素都使用相同的操作进行更新。
# 这意味着我们可以把它并行化。

# ### CPU 上的并行化

# 首先让我们在 CPU 上进行并行化。
# 我们将创建一个“内核函数”(算法的计算核心)，实现两个版本，第一个是顺序版本：

function sequential_add!(y, x)
    for i in eachindex(y, x)
        @inbounds y[i] += x[i]
    end
    return nothing
end

fill!(y, 2)
sequential_add!(y, x)
@test all(y .== 3.0f0)

# 现在是并行的实现版本：

# parallel implementation
function parallel_add!(y, x)
    Threads.@threads for i in eachindex(y, x)
        @inbounds y[i] += x[i]
    end
    return nothing
end

fill!(y, 2)
parallel_add!(y, x)
@test all(y .== 3.0f0)

# 现在，如果我在一台 JULIA_NUM_THREADS=4，意味着至少有4个核的机器上启动 Julia，我将得到以下结果：

@assert Threads.nthreads() == 4     #src

# ```julia
# using BenchmarkTools
# @btime sequential_add!($y, $x)
# ```

# ```
#   487.303 μs (0 allocations: 0 bytes)
# ```

# versus

# ```julia
# @btime parallel_add!($y, $x)
# ```

# ```
#   259.587 μs (13 allocations: 1.48 KiB)
# ```

#可以看到并行化在性能上有提升，但由于启动线程的开销，这一提升并没有提高到原来的4倍。
# 当使用较大的数组时，这种开销会被大量的“实际工作”“稀释”；随着内核数增加，将表现出接近于线性的提升。
# 相反，对于小数组，并行版本可能比串行版本慢。

# ## 你的第一次 GPU 计算

# ### 安装

# 对于本教程的大部分内容，你需要有一台具有兼容 GPU 且安装 [CUDA](https://developer.nvidia.com/cuda-downloads) 的计算机。
# 您还应该使用 Julia [包管理器](https://docs.julialang.org/en/latest/stdlib/Pkg/)安装以下包：

# ```julia
# pkg> add CUDA
# ```

# 如果这是你的第一次使用，通过测试 CUDA.jl 包来测试你的 GPU 是否在工作是个不错的主意。

# ```julia
# pkg> add CUDA
# pkg> test CUDA
# ```

# ### GPU 上的并行化

# 我们将首先演示使用 CuArray 类型的高级 GPU 计算，而不显式地编写内核函数：

using CUDA

x_d = CUDA.fill(1.0f0, N)  # a vector stored on the GPU filled with 1.0 (Float32)
y_d = CUDA.fill(2.0f0, N)  # a vector stored on the GPU filled with 2.0

# 这里 d 的意思是“设备”，与“主机”相对。现在我们来产生一些增量：

y_d .+= x_d
@test all(Array(y_d) .== 3.0f0)

# 语句数组 (y_d) 将 y_d 中的数据移回主机进行测试。
# 如果我们想对它进行基准测试，我们把它放到一个函数中：

function add_broadcast!(y, x)
    CUDA.@sync y .+= x
    return
end

# ```julia
# @btime add_broadcast!($y_d, $x_d)
# ```

# ```
#   67.047 μs (84 allocations: 2.66 KiB)
# ```

# 最有趣的部分是对 `CUDA.@sync` 的调用。
# CPU 可以分配任务给 GPU，然后去做其他的事情（比如分配*更多的*任务给 GPU），而 GPU 完成它的任务。
# 在 `CUDA.@sync` 块中封装执行会使CPU阻塞，直到队列中的 GPU 任务完成，类似 `Base.@sync` 等待 CPU 任务的方式。
# 如果没有这样的同步，你需要考虑的是启动计算所花费的时间，而不是执行计算的时间。
# 但大多数时候，你不需要显式同步：有许多执行方法，如从 GPU 复制内存到 CPU，隐式同步执行。

# 从这个特殊的计算机和 GPU 中，你可以看到 GPU 的计算速度明显快于单线程的 CPU 计算，
# 而 CPU 多线程的使用使得用 CPU 实现的版本具有竞争力。
# 根据硬件的不同，可能会得到不同的结果。

# ### 编写你的第一个 GPU 核函数

# 使用高级的 GPU 数组函数使得在 GPU 上执行这一计算变得很容易。
# 然而，我们还并不了解底层的工作原理，而这正是本教程的主要目标。
# 那么让我们用 GPU 核函数实现相同的功能：

function gpu_add1!(y, x)
    for i = 1:length(y)
        @inbounds y[i] += x[i]
    end
    return nothing
end

fill!(y_d, 2)
@cuda gpu_add1!(y_d, x_d)
@test all(Array(y_d) .== 3.0f0)

# 除了使用 `CuArray`s `x_d` 和 `y_d` 之外，唯一与 GPU 相关的部分是通过 `@cuda` 进行的*核函数启动*。
# 当您第一次发出这个 `@cuda` 语句时，它将编译核函数 (`gpu_add1!`) 用于在GPU上执行。
# 一旦编译，以后的调用会很快。
# 你可以用 ?@cuda` 语句从 Julia 提示中查看 `@cuda` 的扩展。

# 让我们测试一下这个：

function bench_gpu1!(y, x)
    CUDA.@sync begin
        @cuda gpu_add1!(y, x)
    end
end

# ```julia
# @btime bench_gpu1!($y_d, $x_d)
# ```

# ```
#   119.783 ms (47 allocations: 1.23 KiB)
# ```

# 这比上面基于 broadcasting 的版本要*慢得多*。发生了什么？

# ### 分析

# 当您没有获得预期的性能时，通常您的第一步应该是分析代码并查看它在哪里花费了时间。
# 为此，你需要能够运行 NVIDIA的 [`nvprof`tool](https://devblogs.nvidia.com/cuda-pro-tip-nvprof-your-handy-universal-gpu-profiler/) 工具。
# 在 Unix 系统上，这样启动 Julia：
#
# ```sh
# $ nvprof --profile-from-start off /path/to/julia
# ```
#
# 用你的 Julia 二进制文件的路径替换 `/path/to/julia`。
# 注意，我们不会立即启动剖析器，而是调用 CUDA 接口，用 `CUDq.@profile` 手动启动剖析器。
# 这样就排除了编译内核的时间：

bench_gpu1!(y_d, x_d)  # run it once to force compilation
CUDA.@profile bench_gpu1!(y_d, x_d)

# W当我们退出 Julia REPL 时，profiler进程将打印关于执行核函数和接口调用的信息：

# ```
# ==2574== Profiling result:
#             Type  Time(%)      Time     Calls       Avg       Min       Max  Name
#  GPU activities:  100.00%  247.61ms         1  247.61ms  247.61ms  247.61ms  ptxcall_gpu_add1__1
#       API calls:   99.54%  247.83ms         1  247.83ms  247.83ms  247.83ms  cuEventSynchronize
#                     0.46%  1.1343ms         1  1.1343ms  1.1343ms  1.1343ms  cuLaunchKernel
#                     0.00%  4.9490us         1  4.9490us  4.9490us  4.9490us  cuEventRecord
#                     0.00%  4.4190us         1  4.4190us  4.4190us  4.4190us  cuEventCreate
#                     0.00%     960ns         2     480ns     358ns     602ns  cuCtxGetCurrent
# ```

# 您可以看到100%的时间都花在了 `ptxcall_gpu_add1__1` 上。这是当为了这些输入编译 `gpu_add1!` 时，CUDA.jl 分配的那个核函数的名字。
# （你是否创建了多个数据类型的数组，例如 `xu_d = CUDA.fill(0x01, N)`，
# 你也许还见过 `ptxcall_gpu_add1__2` 和其他的。
# 与 Julia 的其他部分一样，您可以定义单个方法，它将在编译时，针对您正在使用的特定数据类型进行专门化。)

# 为了进一步了解情况，使用 `--print-gpu-trace` 运行分析。
# 你也可以用一个包含所有你想要运行的命令的文件路径参数来调用 Julia（包括对 `CUDA.@profile` 的调用）：
#
# ```sh
# $ nvprof --profile-from-start off --print-gpu-trace /path/to/julia /path/to/script.jl
#      Start  Duration   Grid Size   Block Size     Regs*    SSMem*    DSMem*           Device   Context    Stream  Name
#   13.3134s  245.04ms     (1 1 1)      (1 1 1)        20        0B        0B  GeForce GTX TIT         1         7  ptxcall_gpu_add1__1 [34]
# ```

# T这里要注意的关键是“Grid Size”和“Block Size”列中的`(1 1 1)`。
# 这些术语稍后会进行解释，但是现在，只要指出这表明这个计算是按顺序运行的就足够了。
# 值得注意的是，用 GPU 进行顺序处理要比用 CPU 慢得多；GPU 的亮点在于大规模并行。

# ### 编写一个平行的 GPU 核函数

# 为了加速内核，我们希望并行化它，这意味着将不同的任务分配给不同的线程。 To facilitate the assignment of work, each CUDA thread gets access
# 为了方便分配工作，每个CUDA线程都可以访问表示自己唯一身份的变量，就像 [`Threads.threadid()`](https://docs.julialang.org/en/latest/manual/parallel-computing/#Multi-Threading-(Experimental)-1) 对CPU线程一样。
# `threadid` 和 `nthreads` 在 CUDA 中的类似物分别被称为 `threadIdx` 和 `blockDim`；
# 其中一个不同点在于，它们返回一个包含字段 x、y 和 z 的三维结构，以简化最多用于三维数组的笛卡尔索引。
# 因此，我们可以通过以下方式分配独特的工作：

function gpu_add2!(y, x)
    index = threadIdx().x    # this example only requires linear indexing, so just use `x`
    stride = blockDim().x
    for i = index:stride:length(y)
        @inbounds y[i] += x[i]
    end
    return nothing
end

fill!(y_d, 2)
@cuda threads=256 gpu_add2!(y_d, x_d)
@test all(Array(y_d) .== 3.0f0)

# 注意这里的 `threads=256`，它将工作划分为 256 个用线性模式编号的线程。
# （对于二维数组，我们可以使用 `threads=(16, 16)`，然后`x`和`y`都是相关的。）

# 现在让我们尝试对它进行测试：

function bench_gpu2!(y, x)
    CUDA.@sync begin
        @cuda threads=256 gpu_add2!(y, x)
    end
end

# ```julia
# @btime bench_gpu2!($y_d, $x_d)
# ```

# ```
#   1.873 ms (47 allocations: 1.23 KiB)
# ```

# 好多了！

# 但显然，我们还有很长的路要走，以达到最初用 broadcasting 实现的结果。
# 为了做得更好，我们需要更多的并行化。GPU 在*单流多处理器*（SM）上运行的线程数量有限，但与此同时它也有多个 SM。
# 为了充分利用它们，我们需要运行一个具有多个*块*的核函数。
# 我们将这样分配工作:
#
# ![block grid](intro1.png)
#
# 此图[借用了 C/ C++ 库的描述](https://devblogs.nvidia.com/even-easier-introduction-cuda/);
# 在Julia中，线程和块从1而不是0开始编号。
# 在这个图中，由 256 个线程组成的 4096 块（即 1048576 = 2^20 个线程）确保每个线程只递增一个条目；
# 但是，为了确保可以处理任意大小的数组，我们仍然用了一个循环：

function gpu_add3!(y, x)
    index = (blockIdx().x - 1) * blockDim().x + threadIdx().x
    stride = blockDim().x * gridDim().x
    for i = index:stride:length(y)
        @inbounds y[i] += x[i]
    end
    return
end

numblocks = ceil(Int, N/256)

fill!(y_d, 2)
@cuda threads=256 blocks=numblocks gpu_add3!(y_d, x_d)
@test all(Array(y_d) .== 3.0f0)

# The benchmark:

function bench_gpu3!(y, x)
    numblocks = ceil(Int, length(y)/256)
    CUDA.@sync begin
        @cuda threads=256 blocks=numblocks gpu_add3!(y, x)
    end
end

# ```julia
# @btime bench_gpu3!($y_d, $x_d)
# ```

# ```
#   67.268 μs (52 allocations: 1.31 KiB)
# ```

# 最后，我们取得了与 broadcast 版本差不多的性能。
# 让我们再次运行 `nvprof` 来确认这个启动配置：
#
# ```
# ==23972== Profiling result:
#    Start  Duration            Grid Size      Block Size     Regs*    SSMem*    DSMem*           Device   Context    Stream  Name
# 13.3526s  101.22us           (4096 1 1)       (256 1 1)        32        0B        0B  GeForce GTX TIT         1         7  ptxcall_gpu_add3__1 [34]
# ```

# ### 打印

# 在调试时，需要打印一些值是很常见的。 这是通过 `@cuprint` 实现的：

function gpu_add2_print!(y, x)
    index = threadIdx().x    # this example only requires linear indexing, so just use `x`
    stride = blockDim().x
    @cuprintln("thread $index, block $stride")
    for i = index:stride:length(y)
        @inbounds y[i] += x[i]
    end
    return nothing
end

@cuda threads=16 gpu_add2_print!(y_d, x_d)
synchronize()

# 注意，打印输出只有在使用 `synchronize()` 同步整个 GPU 时才生成。
# 这和 `CUDA.@sync` 很相似，是 `cudaDeviceSynchronize` 在 CUDA C++中的对应。

# ### 错误处理

# 本介绍的最后一个主题是关于错误的处理。
# 注意，上面的核函数使用了 `@inbounds`，但是没有检查 `y` 和 `x` 是否有相同的长度。
# 如果你的内核不遵守这些范围，你将会陷入严重的错误：

# ```
# ERROR: CUDA error: an illegal memory access was encountered (code #700, ERROR_ILLEGAL_ADDRESS)
# 堆栈跟踪：
#  [1] ...
# ```

# 如果删除 `@inbounds` 注释，则会得到
#
# ```
# ERROR: a exception was thrown during kernel execution.
#        在调试级别 2 上运行 Julia 以跟踪设备堆栈。
# ```

# 正如错误消息所提到的，更高级别的调试信息将产生更详细的报告。
# 让我们运行与 `-g2` 相同的代码：
#
# ```
# ERROR: a exception was thrown during kernel execution.
# 堆栈跟踪：
#  [1] throw_boundserror at abstractarray.jl:484
#  [2] checkbounds at abstractarray.jl:449
#  [3] setindex! at /home/tbesard/Julia/CUDA/src/device/array.jl:79
#  [4] some_kernel at /tmp/tmpIMYANH:6
# ```

# !!! 警告
#
#     在旧的 GPU 上（计算能力低于 `sm_70`），这些错误是致命的，且会有效地杀死 CUDA 环境。
#     在这样的 GPU 上，使用运行在 CPU 上的代码来执行你的“完整性检查”通常是一个好主意，只有当你认为它安全的时候才将计算交给 GPU。

# ## 总结

# 请记住，CUDA 的高级功能通常意味着您不需要担心如何在底层编写核函数。
# 然而，在许多情况下，可以使用聪明的底层操作来优化计算。
# 希望你现在适应了这种大胆的操作。


text:
# For converting a number of molecules to molars, given the cell volume
# Arguments:  molecules = number of molecules
#             volume = volume in metres cubed
# Returns:    Float64, the molar density

function moleculesToMolars(molecules::Float64, volume::Float64)
  N_av = 6.02214085774e+23 # avagadro's number
  volume_litres = 1000 * volume
  molars = molecules/(N_av * volume_litres)
  return molars
end

function moleculesToMolars(molecules::Array{Float64, 1}, volume::Float64)
  N_av = 6.02214085774e+23 # avagadro's number
  volume_litres = 1000 * volume
  molars = molecules/(N_av * volume_litres)
  return molars
end

function moleculesToMolars(molecules::Array{Float64, 2}, volume::Float64)
  N_av = 6.02214085774e+23 # avagadro's number
  volume_litres = 1000 * volume
  molars = molecules/(N_av * volume_litres)
  return molars
end


text:
# =========================================================================
#    BSEARCHING MECHANISMS
# =========================================================================

function integerCheck(n::node)
    integer = true
    θ = value.(master(n).obj_dict[:θ])

    for val in θ
        if !issinteger(val)
            integer = false
            break
        end
    end

    return integer
end

function leaf(n::node,upperBound::Float64,maxiter::Float64)
    println("leaf search initiated.")
    #INITIALIZE DFS
    terminate = false
    stack = Vector{node}()
    visited = Vector{node}()
    iter = 0

    #PUSH ROOT TO STACK
    push!(stack,n)

    #ITERATION
    while !(terminate || isempty(stack) || iter > maxiter)
        #POP FROM STACK
        u = pop!(stack)

        #CHECK AND PROCESS
        if u.status[end] == "UNVISITED"
            #PROCESS THE NODE
            colGen(u;track=false,maxCG=50.0)

            #NODE STATUSES
            if u.status[end] == "INTEGER"
                obj = objective_value(master(u))

                #CHANGE UPPERBOUND
                if obj < upperBound
                    upperBound = obj
                    push!(u.status,"NEW_BOUND")
                end

                #ALWAYS TERMINATE IF INTEGER NODE IS FOUND
                terminate = true
                println("found integer.")
            elseif u.status[end] == "EVALUATED"
                obj = objective_value(master(u))

                #BRANCH IF NOT ABOVE uB
                if obj <= upperBound
                    append!(stack,createBranch(u))
                end
            end

            #SEND TO VISITED
            push!(visited,u)
            iter += 1
        end
    end

    return (stack=stack,visited=visited,upperBound=upperBound)
end

function createBranch(n::node)
    branches = Vector{node}()
    R = Dict(1:length(n.columns) .=> n.columns)
    θ = value.(master(n).obj_dict[:θ])

    seeds = find_branch(R,θ)
    if !isempty(Q(seeds,R))
        val = s(seeds,R,θ)
        println("branch on $seeds: $val")

        if val - floor(val) < 0.5
            signs = [">=","<="]
        else
            signs = ["<=",">="]
        end

        for br in signs
            push!(branches,
                node(
                    n.self, #parent
                    uuid1(), #self

                    vcat(n.bounds, #bounds
                        bound(
                            seeds,br,
                            if br == "<="
                                floor(s(seeds,R,θ))
                            else
                                ceil(s(seeds,R,θ))
                            end
                        )
                    ),
                    deepcopy(n.columns), #columns
                    initStab(), #stabilizer
                    ["UNVISITED"]
                )
            )
        end
    end

    return branches
end

function traverse(max::Float64)
    #VALUE STORE
    stack = Vector{node}()
    visited = Vector{node}()
    upperBound = Inf

    #FLOW CONTROL
    terminate = false

    #CREATE ROOT
    newroot = root()

    #PUSH ROOT TO STACK
    push!(stack,newroot)

    #ITERATION
    while !(terminate || isempty(stack))
        #POP FIRST ELEMENT OF STACK
        n = pop!(stack)

        #FIND LEAVE (INTEGER NODE) FROM NODE
        results = leaf(n,upperBound,Inf)

        #USE THE RESULTS
        append!(stack,results.stack)
        append!(visited,results.visited)
        if results.upperBound < upperBound
            upperBound = results.upperBound
        end

        #PRINT CURRENT CONDITIONS
        println("size of stack: $(length(stack))")
        println("size of visited: $(length(visited))")
        println("current upperBound: $upperBound")

        #MAX ITER REACHED
        if length(visited) >= max
            terminate = true
        end
    end

    return stack,visited,upperBound
end


text:
module UtilityFunctions890

# Write your package code here.

# Code for utility functions in the Permanent Income Model

# This is the interfaca (API) visible from outside of the module.
export AbstractUtility, UtilityLog, UtilityCRRA
export utility, marg_utility, inv_utility, inv_marg_utility, c_growth, euler_dev

## ----------  Generic
# Here we define and document the interface that is common
# to all concrete types

abstract type AbstractUtility end

"""
    utility(u, c)

Utility level.

This is a generic documentation that applies to all concrete types.
But there is no generic implementation, hence no function body (no methods yet).
"""
function utility end

"""
    marg_utility(u, c)

Marginal utility. Also a generic documentation.
"""
function marg_utility end

"""
	euler_dev(u, cV, betaR)

Euler equation deviation. Generic!
Returns u'(c) / u'(c') - β R

This is a generic function with an implementation that is common to all concrete types.
So it has a function body.
"""
function euler_dev(
    u :: AbstractUtility,
    cV :: AbstractVector{Float64},
    betaR :: Float64
    )
    muV = marg_utility(u, cV);
    T = length(cV);
    devV = muV[1 : (T-1)] ./ muV[2 : T] .- betaR;
    return devV
end

## ----------  Log

struct UtilityLog <: AbstractUtility end

# Note the broadcasting dot. This now works for scalars and arrays.
utility(u :: UtilityLog, c) = log.(c);
marg_utility(u :: UtilityLog, c) = 1.0 ./ c;
inv_utility(u :: UtilityLog, util) = exp.(util);
inv_marg_utility(u :: UtilityLog, mu) = 1.0 ./ mu;
c_growth(u :: UtilityLog, betaR) = betaR;

## -----------  CRRA

struct UtilityCRRA <: AbstractUtility
    sigma :: Float64
end

utility(u :: UtilityCRRA, c) = (c .^ (1.0 - u.sigma)) ./ (1.0 - u.sigma) .- 1.0;
marg_utility(u :: UtilityCRRA, c) = c .^ (-u.sigma);
inv_utility(u :: UtilityCRRA, util) = ((util .+ 1.0) .* (1.0 - u.sigma)) .^ (1.0 / (1.0 - u.sigma));
inv_marg_utility(u :: UtilityCRRA, mu) = mu .^ (-1.0 / u.sigma);
c_growth(u :: UtilityCRRA, betaR) = betaR .^ (1.0 / u.sigma);

end # module


text:
module Automata

# include("util.jl")
# include("history.jl")
# include("render.jl")
# include("elementary.jl")
# include("totalistic.jl")
# include("continuous.jl")

_bounds = Dict(
    "nn" => ((-1, -1), (-1, 1)),
    "nn2t" => ((-2, -1), (-1, 1))
)

export Rule
struct Rule{T, N, F}
    map_fn::F
    bounds::Tuple{Tuple{Int, Int}, Tuple{Int, Int}}
    default_initial_state::Array{T, N}
end

export ElementaryRule
function ElementaryRule(rule::Int)
    rule = UInt8(rule)
    rule_map = zeros(Bool, 2^3)
    for i = 1:length(rule_map)
        rule_map[i] = rule >> (i - 1) & 1
    end

    map_fn = @inline (v) -> rule_map[v[1, 1] << 2 | v[1, 2] << 1 | v[1, 3] + 1]
    return Rule(map_fn, _bounds["nn"], reshape(Bool[ true ], 1, 1))
end

export ElementaryReversibleRule
function ElementaryReversibleRule(rule::Int)
    rule = UInt8(rule)
    rule_map = zeros(Bool, 2^4)
    for p = 0:1
        for i = 1:2^3
            rule_map[p * 8 + i] = rule >> (i - 1) & 1
            if p == 1
                rule_map[p * 8 + i] = !rule_map[p * 8 + i]
            end
        end
    end

    map_fn = @inline (v) -> rule_map[v[1, 2] << 3 | v[2, 1] << 2 | v[2, 2] << 1 | v[2, 3] + 1]
    return Rule(map_fn, _bounds["nn2t"], reshape(Bool[ false true ], 2, 1))
end

export TotalisticRule
function TotalisticRule(code::Int, colors::Int)
    code_count = 3 * (colors - 1) + 1
    code_map = zeros(UInt8, code_count)
    for i in 1:code_count
        code_map[i] = (code ÷ Int64(colors)^(i-1)) % Int64(colors)
    end

    map_fn = @inline (v) -> code_map[sum(v) + 1]
    return Rule(map_fn, _bounds["nn"], reshape(UInt8[ 1 ], 1, 1))
end

export ContinuousRule
function ContinuousRule(map_val_fn)
    map_fn = @inline (v) -> map_val_fn(sum(v))
    return Rule(map_fn, _bounds["nn"], reshape(Float32[ 1 ], 1, 1))
end

export automaton
function automaton(rule::Rule, steps::Int, initial_state::Vector{T}; edges=nothing) where {T}
    return automaton(rule, steps, reshape(initial_state, 1, length(initial_state)); edges=edges)
end
function automaton(rule::Rule{T, 2}, steps::Int; edges=nothing) where {T}
    return automaton(rule, steps, rule.default_initial_state; edges=nothing)
end
function automaton(rule::Rule{T, 2}, steps::Int, initial_state::Array{T, 2}; edges=nothing) where {T}
    rule_fn = rule.map_fn
    bounds = rule.bounds

    if isnothing(edges)
        edges = (steps - 1)
    end
    space = zeros(T, steps, edges * 2 + size(initial_state)[2])

    slice_y, slice_x = bounds
    space[1:size(initial_state)[1], (edges + 1):(edges + size(initial_state)[2])] .= initial_state

    width = size(space)[2]

    left_start = width+slice_x[1]+1
    left_end = 1+slice_x[2]
    right_start = width+slice_x[1]
    right_end = slice_x[2]

    @views for i = (size(initial_state)[1]+1):steps
        vertical_slice = (i+slice_y[1]):(i+slice_y[2])
        for j = 2:(edges*2-1+size(initial_state)[2])
            space[i, j] = rule_fn(space[vertical_slice, (j+slice_x[1]):(j+slice_x[2])])
        end
        space[i, 1] = rule_fn(hcat(space[vertical_slice, left_start:width], space[vertical_slice, 1:left_end]))
        space[i, end] = rule_fn(hcat(space[vertical_slice, right_start:width], space[vertical_slice, 1:right_end]))
    end

    return space
end

export automaton_simple
function automaton_simple(rule::UInt8, steps::Int, initial_state::Union{Array{T, 2}, Nothing}=nothing; edges=nothing) where {T}
    rule_fn = rule.map_fn
    bounds = rule.bounds

    if isnothing(initial_state)
        initial_state = rule.default_initial_state
    end

    if isnothing(edges)
        edges = (steps - 1)
    end
    space = zeros(T, steps, edges * 2 + size(initial_state)[2])

    slice_y, slice_x = bounds
    space[1:size(initial_state)[1], (edges + 1):(edges + size(initial_state)[2])] .= initial_state

    width = size(space)[2]

    left_start = width+slice_x[1]+1
    left_end = 1+slice_x[2]
    right_start = width+slice_x[1]
    right_end = slice_x[2]

    for i = (size(initial_state)[1]+1):steps
        vertical_slice = (i+slice_y[1]):(i+slice_y[2])
        for j = 2:(edges*2-1+size(initial_state)[2])
            v = space[vertical_slice, (j+slice_x[1]):(j+slice_x[2])]
            space[i, j] = rule_fn(v)
        end
        space[i, 1] = rule_fn(hcat(space[vertical_slice, left_start:width], space[vertical_slice, 1:left_end]))
        space[i, end] = rule_fn(hcat(space[vertical_slice, right_start:width], space[vertical_slice, 1:right_end]))
    end

    return space
end

end # module


text:
using VectorizedRNG
using Test

using InteractiveUtils: versioninfo
versioninfo(verbose=true)

using RNGTest, Random, SpecialFunctions, Aqua, Distributions

const α = 1e-4

function smallcrushextrema(res)
    r1 = Base.Cartesian.@ntuple 5 i -> (res[i])::Float64
    r2 = res[6]::Tuple{Float64,Float64}
    r3 = Base.Cartesian.@ntuple 3 i -> (res[i+6])::Float64
    r4 = res[10]::NTuple{5,Float64}
    mi = min(
        minimum(r1), minimum(r2), minimum(r3), minimum(r4)
    )
    ma = max(
        maximum(r1), maximum(r2), maximum(r3), maximum(r4)
    )
    mi, ma
end

const INVSQRT2 = Float64(1/sqrt(big(2)))
# TODO: Get a new SIMD erf implementation
# @inline function normalcdf(v::VectorizedRNG.Vec{W,T}) where {W,T}
#     T(0.5) * ( one(T) + VectorizedRNG.SIMDPirates.verf( v * INVSQRT2 ) )
# end
# function normalcdf!(x::AbstractVector{T}) where {T}
#     _W, Wshift = VectorizedRNG.VectorizationBase.pick_vector_width_shift(T)
#     W = VectorizedRNG.VectorizationBase.pick_vector_width(T)
#     N = length(x)
#     ptrx = pointer(x)
#     i = 0
#     for _ ∈ 1:(N >>> Wshift)
#         ptrxᵢ = VectorizedRNG.VectorizationBase.gep(ptrx, i)
#         v = VectorizedRNG.SIMDPirates.vload(W, ptrxᵢ)
#         VectorizedRNG.SIMDPirates.vstore!(ptrxᵢ, normalcdf(v))
#         i += _W
#     end
#     if i < N
#         ptrxᵢ = VectorizedRNG.VectorizationBase.gep(ptrx, i)
#         mask = VectorizedRNG.VectorizationBase.mask(T, N & (_W - 1))
#         v = VectorizedRNG.SIMDPirates.vload(W, ptrxᵢ, mask)
#         VectorizedRNG.SIMDPirates.vstore!(ptrxᵢ, normalcdf(v), mask)
#     end
#     x
# end
normalcdf(x::T) where {T} = T(0.5) * ( one(T) + erf( x * INVSQRT2 ) )
function normalcdf!(x::AbstractVector{T}) where {T}
    @. x = normalcdf(x)
end

    struct RandNormal01{T<:VectorizedRNG.AbstractVRNG} <: Random.AbstractRNG
        rng::T
    end
    function Random.rand!(r::RandNormal01, x::AbstractArray)
        randn!(r.rng, x)
        normalcdf!(x)
        # x .= cdf.(Normal(0,1), x)
    end
    # rngnorm = RNGTest.wrap(RandNormal01(local_pcg()), Float64);

function test_serial_rng(f)
  res = RNGTest.smallcrushJulia(f)
  mi, ma = smallcrushextrema(res)
  @show mi, ma
  @test mi > α
  @test ma < 1 - α
end

@testset "VectorizedRNG.jl" begin
    Aqua.test_all(VectorizedRNG)#, ambiguities = VERSION < v"1.6-DEV")

    @test isempty(detect_unbound_args(VectorizedRNG))
    @testset "Small Crush" begin
        rngunif64 = RNGTest.wrap(local_rng(), Float64);
        res = RNGTest.smallcrushJulia(rngunif64)
        mi, ma = smallcrushextrema(res)
        @show mi, ma
        @test mi > α
        @test ma < 1 - α

        rngunif32 = RNGTest.wrap(local_rng(), Float32);
        res = RNGTest.smallcrushJulia(rngunif32)
        mi, ma = smallcrushextrema(res)
        @show mi, ma
        @test mi > α
        @test ma < 1 - α

        rngnorm = RNGTest.wrap(RandNormal01(local_rng()), Float64);
        res = RNGTest.smallcrushJulia(rngnorm)
        mi, ma = smallcrushextrema(res)
        @show mi, ma
        @test mi > α
        @test ma < 1 - α

        A = zeros(13, 29);
        randn!(local_rng(), A);
        @test iszero(sum(iszero, A))
        A32 = zeros(Float32, 13, 29);
        randn!(local_rng(), A32);
        @test iszero(sum(iszero, A32))
        # TODO: Support this again
        # rngnorm = RNGTest.wrap(RandNormal01(local_rng()), Float32);
        # res = RNGTest.smallcrushJulia(rngnorm)
        # mi, ma = smallcrushextrema(res)
        # @show mi, ma
        # @test mi > α
        # @test ma < 1 - α

        # scalar mode tests
        test_serial_rng(() -> rand(local_rng()))
        test_serial_rng(() -> normalcdf(randn(local_rng())))
        test_serial_rng(() -> cdf(Gamma(), rand(local_rng(), Gamma())))
        test_serial_rng(() -> VectorizedRNG.floatbitmask((rand(local_rng(), UInt128) >> 64) % UInt64, Float64) - VectorizedRNG.oneopenconst(Float64))

    end
    @testset "Discontiguous in place" begin
        x = zeros(5, 117); xv = view(x, 5, :)
        rand!(local_rng(), xv)
        @test !any(iszero, xv)
        @test all(iszero, view(x, 1:4, :))
    end
    @testset "Scaled sampling" begin
        A = Matrix{Float64}(undef, 89, 100);
        randn!(local_rng(), A, VectorizedRNG.StaticInt(0), 5, 100);
        s,l = extrema(A)
        @test s < -100
        @test l > 100
        randn!(local_rng(), A, VectorizedRNG.StaticInt(0), 100, 10);
        @test 90 < sum(A) / length(A) < 110
        randn!(local_rng(), A, VectorizedRNG.StaticInt(1), 100, 10);
        @test 190 < sum(A) / length(A) < 210
    end
end


text:
# This script tests the internal functions of eqcond.jl
using Test, OrderedCollections, HDF5
include("../../../src/includeall.jl")

verbose = :high # set to :none to avoid printing output
rp = joinpath(dirname(@__FILE__), "../../reference/models/li2020")

@testset "Solve Li2020" begin
    # Try with default method
    m = Li2020()
    vars_for_error = [:p, :Q̂]
    stategrid_default, funcvar_default, derivs_default, endo_default, flag_default = solve(m, verbose = verbose,
                                                                                vars_for_error = vars_for_error)

    # Now try with different error approaches and ensure similar convergence
    m <= Setting(:error_method, :Linf)
    m <= Setting(:tol, 1e-5)
    stategrid_Linf, funcvar_Linf, derivs_Linf, endo_Linf, flag_Linf = solve(m, verbose = verbose,
                                                                            vars_for_error = vars_for_error)

    m <= Setting(:error_method, :L2)
    m <= Setting(:tol, 1e-10)
    stategrid_L2, funcvar_L2, derivs_L2, endo_L2, flag_L2 = solve(m, verbose = verbose,
                                                                  vars_for_error = vars_for_error)

    @test flag_default
    @test flag_Linf
    @test flag_L2
    for k in [:p, :Q̂]
        @test funcvar_default[k] ≈ funcvar_L2[k] atol=1e-5
        @test funcvar_default[k] ≈ funcvar_Linf[k] atol=1e-5
    end

    @test funcvar_default[:p] == h5read(joinpath(rp, "solve.h5"), "p")
    @test funcvar_default[:Q̂] == h5read(joinpath(rp, "solve.h5"), "Qhat")
end

#=
# Uncomment to re-write output
h5open(joinpath(rp, "solve.h5"), "w") do file
    write(file, "p", funcvar_default[:p])
    write(file, "Qhat", funcvar_default[:Q̂])
end
=#

nothing


text:
module ParquetFiles

using Parquet, IteratorInterfaceExtensions, TableTraits, FileIO
import IterableTables, DataValues, TableShowUtils

export load, File, @format_str

struct ParquetFile
    filename::String
end

function Base.show(io::IO, source::ParquetFile)
    TableShowUtils.printtable(io, getiterator(source), "Parquet file")
end

function Base.show(io::IO, ::MIME"text/html", source::ParquetFile)
    TableShowUtils.printHTMLtable(io, getiterator(source))
end
Base.Multimedia.showable(::MIME"text/html", source::ParquetFile) = true

function Base.show(io::IO, ::MIME"application/vnd.dataresource+json", source::ParquetFile)
    TableShowUtils.printdataresource(io, getiterator(source))
end
Base.Multimedia.showable(::MIME"application/vnd.dataresource+json", source::ParquetFile) = true

struct ParquetNamedTupleIterator{T,T_row}
    rc::RecCursor
    nrows::Int
end

function Base.eltype(itr::ParquetNamedTupleIterator{T,T_row}) where {T,T_row}
    return T
end

function Base.length(itr::ParquetNamedTupleIterator)
    return itr.nrows
end

@generated function Base.iterate(itr::ParquetNamedTupleIterator{T,T_row}, state...) where {T,T_row}
    names = fieldnames(T)
    quote
        y = iterate(itr.rc, state...)
        if y === nothing
            return nothing
        else
            v = y[1]
            next_state = y[2]
            return T(($([fieldtype(T, i) <: String ? :(String(copy(v.$(names[i])))) : :(v.$(names[i])) for i = 1:length(names)]...),)), next_state
        end
    end
end

function fileio_load(f::FileIO.File{FileIO.format"Parquet"})
    return ParquetFile(f.filename)
end

IteratorInterfaceExtensions.isiterable(x::ParquetFile) = true
TableTraits.isiterabletable(x::ParquetFile) = true

function IteratorInterfaceExtensions.getiterator(file::ParquetFile)
    p = ParFile(file.filename)

    T_row_name = Symbol("RCType$(String(gensym())[3:end])")

    schema(JuliaConverter(ParquetFiles), p, T_row_name)

    T_row = eval(T_row_name)

    col_names = fieldnames(T_row)
    col_types = [i <: Vector{UInt8} ? String : i for i in T_row.types]

    T = NamedTuple{(col_names...,),Tuple{col_types...}}

    rc = RecCursor(p, 1:nrows(p), colnames(p), JuliaBuilder(p, T_row))

    it = ParquetNamedTupleIterator{T,T_row}(rc, nrows(p))

    return it
end

end # module


text:
using Compat
using Test

@test isempty(detect_ambiguities(Base, Core, Compat))

# julia#29679
@test !isnothing(1)
@test isnothing(nothing)

# https://github.com/JuliaLang/julia/pull/29749
@testset "row/column/slice iterators" begin
    # Simple ones
    M = [1 2 3; 4 5 6; 7 8 9]
    @test collect(eachrow(M)) == collect(eachslice(M, dims = 1)) == [[1, 2, 3], [4, 5, 6], [7, 8, 9]]
    @test collect(eachcol(M)) == collect(eachslice(M, dims = 2)) == [[1, 4, 7], [2, 5, 8], [3, 6, 9]]
    @test_throws DimensionMismatch eachslice(M, dims = 4)

    # Higher-dimensional case
    M = reshape([(1:16)...], 2, 2, 2, 2)
    @test_throws MethodError collect(eachrow(M))
    @test_throws MethodError collect(eachcol(M))
    @test collect(eachslice(M, dims = 1))[1][:, :, 1] == [1 5; 3 7]
end

# Support for positional `stop`
@test range(0, 5, length = 6) == 0.0:1.0:5.0
@test range(0, 10, step = 2) == 0:2:10
@test_throws ArgumentError range(0, 10)

mutable struct TLayout
    x::Int8
    y::Int16
    z::Int32
end
tlayout = TLayout(5,7,11)
@test hasfield(TLayout, :y)
@test !hasfield(TLayout, :a)
@test hasproperty(tlayout, :x)
@test !hasproperty(tlayout, :p)

@test merge((a=1,b=1)) == (a=1,b=1)
@test merge((a=1,), (b=2,), (c=3,)) == (a=1,b=2,c=3)

@testset "only" begin
    @test only([3]) === 3
    @test_throws ArgumentError only([])
    @test_throws ArgumentError only([3, 2])

    @test @inferred(only((3,))) === 3
    @test_throws ArgumentError only(())
    @test_throws ArgumentError only((3, 2))

    @test only(Dict(1=>3)) === (1=>3)
    @test_throws ArgumentError only(Dict{Int,Int}())
    @test_throws ArgumentError only(Dict(1=>3, 2=>2))

    @test only(Set([3])) === 3
    @test_throws ArgumentError only(Set(Int[]))
    @test_throws ArgumentError only(Set([3,2]))

    @test @inferred(only((;a=1))) === 1
    @test_throws ArgumentError only(NamedTuple())
    @test_throws ArgumentError only((a=3, b=2.0))

    @test @inferred(only(1)) === 1
    @test @inferred(only('a')) === 'a'
    @test @inferred(only(Ref([1, 2]))) == [1, 2]
    @test_throws ArgumentError only(Pair(10, 20))

    @test only(1 for ii in 1:1) === 1
    @test only(1 for ii in 1:10 if ii < 2) === 1
    @test_throws ArgumentError only(1 for ii in 1:10)
    @test_throws ArgumentError only(1 for ii in 1:10 if ii > 2)
    @test_throws ArgumentError only(1 for ii in 1:10 if ii > 200)
end

# https://github.com/JuliaLang/julia/pull/32628
@testset "mod with ranges" begin
    for n in -10:10
        @test mod(n, 0:4) == mod(n, 5)
        @test mod(n, 1:5) == mod1(n, 5)
        @test mod(n, 2:6) == 2 + mod(n-2, 5)
        @test mod(n, Base.OneTo(5)) == mod1(n, 5)
    end
    @test mod(Int32(3), 1:5) == 3
    @test mod(big(typemax(Int))+99, 0:4) == mod(big(typemax(Int))+99, 5)
    @test_throws MethodError mod(3.141, 1:5)
    @test_throws MethodError mod(3, UnitRange(1.0,5.0))
    @test_throws MethodError mod(3, 1:2:7)
    @test_throws DivideError mod(3, 1:0)
end

# https://github.com/JuliaLang/julia/pull/33568
@testset "function composition" begin
    @test ∘(x -> x-2, x -> x-3, x -> x+5)(7) == 7
    fs = [x -> x[1:2], uppercase, lowercase]
    @test ∘(fs...)("ABC") == "AB"
end

nothing


text:
using GroundEffects
using Test

statements(block) = filter(ex -> ex isa Expr, block.args)

expressions_to_be_lowered = quote
    [1, 2, 3, 4]
    [1; 2]
    [1 2; 3 4]
    [1 2]
    Any[1; 2]
    Any[1 2; 3 4]
    Any[1 2]
    ones(2) .+ 2
    identity.(ones(2) .+ 2)
    [1, 2] .== [1, 3]
    [1, 2] .!= [1, 3]
    [1, 2][1]
    [1, 2][end]
    [1, 2][end ÷ 2]
    begin
        x = [1, 2, 3]
        ans = begin
            x[1] = 10
        end
        (ans, x)
    end
    begin
        x = [1, 2, 3]
        ans = begin
            x[end] = 30
        end
        (ans, x)
    end
    begin
        x = [1, 2, 3]
        ans = begin
            x[end÷2] = 10
        end
        (ans, x)
    end
    begin
        x = zeros(3)
        x .= 1
    end
    begin
        x = [1, 2, 3]
        x .+= 2
    end
    im.re
    identity() do; end()
    identity(identity)() do; end()
end |> statements

@testset for ex in expressions_to_be_lowered
    lex = GroundEffects.lower(ex)
    @test ex != lex
    actual = @eval $lex
    desired = @eval $ex
    @test actual == desired
end

mutable = Text("old")
@testset "setproperty!" begin
    ex = :(mutable.content = "new")
    lex = GroundEffects.lower(ex)
    @test ex != lex
    result = @eval $lex
    @test result == "new"
    @test mutable.content == "new"
end


text:
grid = [
	 [1/24, 1/24, 1/24, 1/24, 1/24]
	,[1/24, 0.0, 0.0, 1/24, 1/24]
	,[1/24, 0.0, 1/24, 1/24, 1/24]
	,[1/24, 0.0, 0.0, 1/24, 1/24]
	,[1/24, 0.0, 1/24, 1/24, 1/24]
	,[1/24, 1/24, 1/24, 1/24, 1/24]

	]

const HEIGHT=6
const WIDTH=5

@enum SquareType OPEN CLOSED
@enum Direction WEST NORTH EAST SOUTH
@enum DriftType STRAIGHT LEFT RIGHT

const Drift=Dict(STRAIGHT=>.7, LEFT=>.15, RIGHT=>.15)
const AllDirects =(WEST, NORTH, EAST, SOUTH)

# Actual => Evidence
const Sense=Dict(OPEN=>Dict(OPEN=>.8, CLOSED=>.2), CLOSED=>Dict(OPEN=>.25, CLOSED=>.75))

function print_grid(grid::Array{Array{Float64,1},1})
	for i = 1:HEIGHT
		for j = 1:WIDTH
			if(!notblocked(grid, (i,j)))
				print("####","\t")
			else
				print(round(grid[i][j]*100, digits=2), "\t")
		   	end
		end
		println()
	end
end

# Check if the posistion is blocked (either out of boudns/obstacles) or a free space
function notblocked(grid::Array{Array{Float64,1},1}, pos::Tuple{Int64, Int64})
	if( (pos[1]>0 && pos[1]<=HEIGHT) && (pos[2]>0 && pos[2]<=WIDTH) )
		return (grid[pos[1]][pos[2]] != 0.0)
	end
	return false
end

# Given direction adn posistions return psostion that is pos+direction
function move(pos::Tuple{Int64, Int64}, dir::Direction)

	if(dir==WEST)
		return (pos[1], pos[2]-1)
	end

	if(dir==NORTH)
		return (pos[1]-1, pos[2])
	end

	if(dir==EAST)
		return (pos[1], pos[2]+1)
	end

	if(dir==SOUTH)
		return (pos[1]+1, pos[2])
	end
end

# Given posistion and evidecne return probabilty of being in that posistion
# PRECONDITION: Grid, position and evidecne
# POSTCONDITION: Conditonal Evidence Probabilty at that point
function evidence_Probabilty(grid::Array{Array{Float64,1},1},
							pos::Tuple{Int64, Int64},
							evidence::Tuple{SquareType, SquareType, SquareType, SquareType})
	prod = 1
	for i in 1:4
		tmp_pos = move(pos,AllDirects[i])
		block = notblocked(grid, tmp_pos)
		if (block)
			prod*= Sense[OPEN][evidence[i]]
		else
			prod*= Sense[CLOSED][evidence[i]]
		end
	end
	prod
end

# Get the evidecne contional probabilty of each posistoin*Pos(s_i)
# Then divide each posistion with the evidnece conditonal probabilty
# PRECONDITION: Given grid and evidence
# POSTCONDITION: Return filter of the grid, usign conditonal evdiecne probabilty
function filter(grid::Array{Array{Float64,1},1}, evidence::Tuple{SquareType, SquareType, SquareType, SquareType})
	tmp_grid = deepcopy(grid)
	for row in 1:6
		for col in 1:5
			if(notblocked(grid, (row,col)))
				tmp_grid[row][col]*=evidence_Probabilty(grid, (row, col), evidence)
			end
		end
	end
	total_sum = sum(sum(tmp_grid))
	tmp_grid / total_sum
end

# Gets the transitional probabilty of posistions going into a point,
# or transitional probablites from the point
# the first is used for prediction, the other for smoothing
# PRECONDITION: Given grid, position, a direction and type of transtional probabitly
# POSTCONDITION: Return probabilty of points going into pos, or probabitly of points going away from pos.
function transprob( grid::Array{Array{Float64,1},1}, pos::Tuple{Int64, Int64}, dir::Direction, getforward=false) # => array of  (pos::Tuple{Int64, Int64}, Drift[Direction] * P(s), Drift[Direction])
	arr = []
	function _gen_parts(straightDir::Direction ,
						leftDir::Direction,
						rightDir::Direction,
						behindDir::Direction,
						grid::Array{Array{Float64,1},1},
						pos::Tuple{Int64,Int64} )
		parent_pos = []
		behind = move(pos, behindDir)
		straight = move(pos, straightDir)
		left = move(pos, leftDir)
		right = move(pos, rightDir)

		if( !notblocked(grid, straight)) # Bounce
			push!(parent_pos, (pos,Drift[STRAIGHT]*grid[pos[1]][pos[2]], Drift[STRAIGHT] ))
		elseif (getforward) # If this is smoothing, want the probabilty in front
			push!(parent_pos, (straight,Drift[STRAIGHT]*grid[straight[1]][straight[2]], Drift[STRAIGHT] ))
		end
		if( notblocked(grid, behind) && !getforward) # Prob of square behind current ot move to current
			push!(parent_pos, (behind,Drift[STRAIGHT]*grid[behind[1]][behind[2]], Drift[STRAIGHT]))
		end
		if(notblocked(grid, left)) # Get probabitly of the left pos coming to curretn
			push!(parent_pos, (left,Drift[LEFT]*grid[left[1]][left[2]], Drift[LEFT] ))
		else #Bounce from left
			push!(parent_pos, (pos,Drift[LEFT]*grid[pos[1]][pos[2]], Drift[LEFT] ))
		end
		if(notblocked(grid, right)) #Get probabilty of right pos coming to current
			push!(parent_pos, (right,Drift[RIGHT]*grid[right[1]][right[2]],Drift[RIGHT]))
		else #Bounce
			push!(parent_pos, (pos,Drift[RIGHT]*grid[pos[1]][pos[2]], Drift[RIGHT]))
		end
	end

	if(dir==WEST)
		arr=_gen_parts(WEST, SOUTH, NORTH, EAST, grid, pos)
	elseif(dir==NORTH)
		arr=_gen_parts(NORTH, EAST, WEST, SOUTH, grid, pos)
	elseif(dir==SOUTH)
		arr=_gen_parts(SOUTH, EAST, WEST, NORTH, grid, pos)
	elseif(dir==EAST)
		arr=_gen_parts(EAST, SOUTH, NORTH, WEST, grid, pos)
	end

	arr
end

# Prediction is the sum of possiple transitional probabilties
# that can reach pos_i, * P(pos).j
# PRECONDITON: Given a grid and direction
# POSTCONDITION: Prediction of where the agent will go after this movement.
function predict(grid::Array{Array{Float64,1}}, dir::Direction)
	tmp_grid = deepcopy(grid)
	for row in 1:6
		for col in 1:5
			if(notblocked(grid, (row,col)))
				x= transprob(grid, (row, col), dir)
				val=sum([x[2] for x in transprob(grid, (row, col), dir)])
				tmp_grid[row][col]=val
			end
		end
	end

	tmp_grid
end

# Get the transitional probabilty of a point going OUT, not in
# an it's conditional probabilty, with it's inital probabilty
# returns 2 things: B at pos, and B*p(s)
# PRECONDITON: Grid
function smoothpart( grid::Array{Array{Float64,1}},
					Bgrid::Array{Array{Float64,1}},
					evidence::Tuple{SquareType, SquareType, SquareType, SquareType},
					dir::Direction, pos::Tuple{Int64, Int64})

	parent_pos=transprob(grid, pos, dir, true)
         x=(sum([evidence_Probabilty(grid, i[1], evidence)* Bgrid[i[1][1]][i[1][2]]* i[3] for i in parent_pos]))
	(x, x * grid[pos[1]][pos[2]])
end

# Get the smoothing part for each posistion in grid
# Then divide the whoel grid by the sum of it's parts
function smooth(                grid::Array{Array{Float64,1}},
				Bgrid::Array{Array{Float64,1}},
				evidence::Tuple{SquareType, SquareType, SquareType, SquareType},
				dir::Direction)

	NewGrid = deepcopy(grid)
	B = deepcopy(Bgrid)

	for row in 1:6
		for col in 1:5
			if(notblocked(grid, (row,col)))
				val=smoothpart(grid, Bgrid, evidence, dir, (row,col))
				B[row][col] = val[1]
				NewGrid[row][col] = val[2]
			else
				B[row][col]  = 0
				NewGrid[row][col] = 0
			end
		end
	end
	NewGrid/=sum(sum(NewGrid))
	(NewGrid, B)
end

tmp  = filter(grid, (OPEN, OPEN, OPEN, OPEN))
tmp2 = predict(tmp, WEST)
tmp3 = filter(tmp2, (CLOSED, CLOSED, OPEN, CLOSED))
tmp4 = predict(tmp3, NORTH)
tmp5 = filter(tmp4, (CLOSED, CLOSED, OPEN, CLOSED))

Bgrid = [
	 [1.0, 1.0, 1.0, 1.0, 1.0]
	 ,[1.0, 1.0, 1.0, 1.0, 1.0]
	 ,[1.0, 1.0, 1.0, 1.0, 1.0]
	 ,[1.0, 1.0, 1.0, 1.0, 1.0]
	 ,[1.0, 1.0, 1.0, 1.0, 1.0]
	 ,[1.0, 1.0, 1.0, 1.0, 1.0]

	 ]

(tmp6,b)=smooth(tmp3,Bgrid,(CLOSED,CLOSED,OPEN,CLOSED),NORTH)
(tmp7,c)=smooth(tmp,b,(CLOSED,CLOSED,OPEN,CLOSED),WEST)

println("Initial Location Probabilities")
print_grid(grid)
println("\nFiltering after Evidence [0, 0, 0, 0]")
print_grid(tmp)
println("\nPrediction after Action W")
print_grid(tmp2)
println("\nFiltering after Evidence [1, 1, 0, 1]")
print_grid(tmp3)
println("\nPrediction after Action N")
print_grid(tmp4)
println("\nFiltering after Evidence [1, 1, 0, 1]")
print_grid(tmp5)
println("\nLast position Smoothing with Evidence [1, 1, 0, 1] and north")
print_grid(tmp6)
println("\nSecond Last posistion smoothing with Evidence [1, 1, 0, 1] And west")
print_grid(tmp7)


text:
using LoopVectorization
using Test

@testset "Utilities" begin
    @test LoopVectorization.isscopedname(:(Base.OneTo), :Base, :OneTo)
    @test LoopVectorization.isscopedname(:(A.B.C.D), (:A, :B, :C), :D)
    @test !LoopVectorization.isscopedname(:(A.B.D),  (:A, :B, :C), :D)
    @test !LoopVectorization.isscopedname(:(A.B.C.D), (:A, :B, :C), :E)
    @test !LoopVectorization.isscopedname(:hello,  :Base, :OneTo)

    @test LoopVectorization.isglobalref(GlobalRef(Base, :getindex), Base, :getindex)
    @test !LoopVectorization.isglobalref(GlobalRef(Base, :getindex), Base, :setindex!)
    @test !LoopVectorization.isglobalref(GlobalRef(Core, :getindex), Base, :getindex)
    @test !LoopVectorization.isglobalref(:getindex, Base, :getindex)
end


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

using Markdown
using InteractiveUtils

# ╔═╡ d09dec61-ab12-4424-bfaf-30f155c90ef9
using Revise

# ╔═╡ 89fcb3a0-db48-11eb-30ac-b501d7fcd55f
using PolaronMobility

# ╔═╡ 4a7ed3cb-88e4-4e86-90c1-285849994ab9
using Plots

# ╔═╡ 8c54eedc-7c8f-4c5d-879e-62743e302801
using DataFrames

# ╔═╡ 3edfd42b-4b8a-4080-b617-740ef3fbf76b
using CSV

# ╔═╡ 27af9453-8044-4b1f-85fd-9e2c5dd104eb
using QuadGK

# ╔═╡ 97aa8c27-f2fe-418c-a662-906dbc126fb6
begin
	ħ = 1.05457162825e-34
	eV = 1.602176487e-19
	me = 9.10938188e-31
	kB =  1.3806504e-23
	ϵ_0 = 8.854E-12
	amu = 1.660_539_066_60e-27
	Ha = 4.35974820e-18
	Bohr = 5.29177249e-11
end

# ╔═╡ d08e66f9-636b-4377-b653-389ee7b78a24
T_range = 100:100:400

# ╔═╡ 50c1a877-807e-498e-87f9-0e909080a452
MAPI= [
# 96.20813558773261 0.4996300522819191
# 93.13630357703363 1.7139631746083817
# 92.87834578121567 0.60108592692181
# 92.4847918585963 0.0058228799414729
# 92.26701437594754 0.100590086574602
# 89.43972834606603 0.006278895133832249
# 46.89209141511332 0.2460894564364346
# 46.420949316788 0.14174282581124137
# 44.0380222871706 0.1987196948553428
# 42.89702947649343 0.011159939465770681
# 42.67180170168193 0.02557751102757614
# 41.46971205834201 0.012555230726601503
# 37.08982543385215 0.00107488277468418
# 36.53555265689563 0.02126940080871224
# 30.20608114002676 0.009019481779712388
# 27.374810898415028 0.03994453721421388
# 26.363055017011728 0.05011922682554448
# 9.522966890022039 0.00075631870522737
4.016471586720514 0.08168931020200264
3.887605410774121 0.006311654262282101
3.5313112232401513 0.05353548710183397
2.755392921480459 0.021303020776321225
2.4380741812443247 0.23162784335484837
2.2490917637719408 0.2622203718355982
2.079632190634424 0.23382298607799906
2.0336707697261187 0.0623239656843172
1.5673011873879714 0.0367465760261409
1.0188379384951798 0.0126328938653956
1.0022960504442775 0.006817361620021601
0.9970130778462072 0.0103757951973341
0.9201781906386209 0.01095811116040592
0.800604081794174 0.0016830270365341532
0.5738689505255512 0.00646428491253749
# 0.022939578929507105 8.355742795827834e-06   # Acoustic modes!
# 0.04882611767873102 8.309858592685e-06
# 0.07575149723846182 2.78248540373041e-05
]

# ╔═╡ e349bf0a-554f-46ae-8e88-8a7b80d0c4cd
# Multiple Phonon Theory

# ╔═╡ 4cd925a7-048a-4696-b2b9-6a6b394e0528
B_scheme = PolaronMobility.HellwarthBScheme(MAPI)

# ╔═╡ b04ccf74-853d-4ab2-886f-5b6c21cc9a51
phonon_mode_freqs = MAPI[:, 1]

# ╔═╡ 4e349b51-0590-4f98-b777-f13db8050b6d
ir_activities = MAPI[:, 2]

# ╔═╡ 1c6fff5f-839a-43e1-acfe-040df2ed4822
ϵ_ionic = [PolaronMobility.ϵ_ionic_mode(f, r, (6.29e-10)^3) for (f, r) in zip(phonon_mode_freqs, ir_activities)]

# ╔═╡ 2f524fa1-c08b-4c11-b05a-e2d016ce6ec0
ϵ_total = sum(ϵ_ionic)

# ╔═╡ 8510da63-37d7-45d0-9012-8dd17745a9e8
α_j = [PolaronMobility.frohlich_α_j(4.5, ϵ_i, ϵ_total, f, 0.12) for (ϵ_i, f) in zip(ϵ_ionic, phonon_mode_freqs)]

# ╔═╡ d9fca06f-ee9c-4790-a059-98667dcd8428
begin
	multi_data = DataFrames.DataFrame(
			alpha = α_j,
			phonon_freqs = phonon_mode_freqs,
			ir_activities = ir_activities,
			ionic = ϵ_ionic
	    )
	CSV.write("multi_data.csv", multi_data)
end

# ╔═╡ 04e16581-ca4f-4456-bea2-2e8bf44a08f8
α_eff = sum(α_j)

# ╔═╡ ba626ebc-d7da-4eca-8149-907d724a2ef5
var_params = [PolaronMobility.multi_variation(i, 4.5, 0.12, (6.29e-10)^3, MAPI; N = 1) for i in T_range]

# ╔═╡ af26a6c8-9c97-4130-bd54-c9fe40e58879
[var_params[i][1] for i in 1:length(T_range)]

# ╔═╡ dd41fbf9-b4b9-46b9-abbf-96eddf59b0ad
 v_j = [var_params[i][1][j] for j in 1:length(MAPI[:, 1]), i in 1:length(T_range)]

# ╔═╡ 2975a6b8-bcec-4e4e-be94-5a3897c650db
begin
	multi_v = DataFrames.DataFrame([[0.0, [i for i in T_range]...]'; [phonon_mode_freqs v_j]], :auto)
	CSV.write("multi_v.csv", multi_v)
end

# ╔═╡ 748d6916-6125-4f9a-adbf-cee3b25dc3ba
w_j = [var_params[i][2][j] for j in 1:length(MAPI[:, 1]), i in 1:length(T_range)]

# ╔═╡ a304f008-b3fa-4a92-a8c1-9ce2f7b4b3df
begin
	multi_w = DataFrames.DataFrame([[0.0, [i for i in T_range]...]'; [phonon_mode_freqs w_j]], :auto)
	CSV.write("multi_w.csv", multi_w)
end

# ╔═╡ 5652b186-0e8d-469c-bfea-bf665f75c06a
F_j = [PolaronMobility.multi_free_energy(var_params[i][1][j, :], var_params[i][2][j, :], T_range[i], 4.5, 0.12, (6.29e-10)^3, MAPI, j) for j in 1:length(MAPI[:, 1]), i in 1:length(T_range)]

# ╔═╡ 2e0eb8f9-3b6b-4464-b715-258c96e22422
begin
	multi_energy = DataFrames.DataFrame([[0.0, [i for i in T_range]...]'; [phonon_mode_freqs F_j]], :auto)
	CSV.write("multi_energy.csv", multi_energy)
end

# ╔═╡ 143b1c45-ce4b-4f3c-9674-ff181c64fe90
F_total = [sum(F_j[:, i]) for i in 1:length(T_range)]

# ╔═╡ 3900999b-fdfa-48ae-afd4-bca5d03008e3
Ω_range = 0.01:0.1:20.01

# ╔═╡ 30d37afa-7b70-479f-8dcc-860af723a96e
β_red = [ħ / kB / T * f * 2π * 1e12 for f in MAPI[:, 1], T in T_range]

# ╔═╡ 9c140cb6-286b-4e25-a3dc-829b1eea9dfb
begin
	multi_beta = DataFrames.DataFrame([[0.0, [i for i in T_range]...]'; [phonon_mode_freqs β_red]], :auto)
	CSV.write("multi_beta.csv", multi_beta)
end

# ╔═╡ 65a76ecd-1dcd-4d5c-88be-1b57d31369da
function h_i(i, v, w) # some vector relating to harmonic eigenmodes
    h = v[i]^2 - w[i]^2
    if length(v) > 1
        for j in 1:length(v)
            if j != i
                h *= (w[j]^2 - v[i]^2) / (v[j]^2 - v[i]^2)
            end
        end
    end
    return h
end

# ╔═╡ 6a793986-844a-4a82-84b0-b71a48aa8242
function D_j(τ, β, v, w) # log of dynamic structure factor for polaron
    D = τ * (1 - τ / β)
    for i in 1:length(v)
        if v[i] != w[i]
        D += (h_i(i, v, w) / v[i]^2) * (2 * sinh(v[i] * τ / 2) * sinh(v[i] * (β - τ) / 2) / (v[i] * sinh(v[i] * β / 2)) - τ * (1 - τ / β))
        end
    end
    return D
end

# ╔═╡ 146bcc43-15fc-45de-bb9b-f16a68d84e41
function multi_memory_function(Ω, β, α, v, w, f, m_eff)

    m_e = 9.10938188e-31; # Electron Mass (kg)
    eV = 1.602176487e-19; # Electron Volt (kg m^2 s^{-2})

    # FHIP1962, page 1009, eqn (36).
    S(t, β_j, v_j, w_j) = cos(t - 1im * β_j / 2) / sinh(β_j / 2) / D_j(-1im * t, β_j, v_j, w_j)^(3 / 2)

    # FHIP1962, page 1009, eqn (35a). Scale Frequency Ω by phonon branch frequency f_j.
    integrand(t, β_j, v_j, w_j, Ω) = (1 - exp(-1im * Ω * t)) * imag(S(t, β_j, v_j, w_j))

    impedence = 0.0
    for j in 1:length(f) # sum over phonon branches
        impedence += 1im * (Ω - 2 * α[j] * (f[j])^2 * quadgk(t -> integrand(t, β[j], v[j], w[j], Ω * 2π / f[j]), 0.0, Inf)[1] / (3 * √π * Ω))
    end
    return impedence / eV * m_eff * m_e * 1e12 / 100^2
end

# ╔═╡ 3c7564b0-b9fe-49e8-9ad7-eaa80e9841a7
begin
	Z = Array{ComplexF64}(undef, length(Ω_range), length(T_range))
	for j in 1:length(T_range), i in 1:length(Ω_range)
		println("Ω: $(Ω_range[i]), T: $(T_range[j])")
		Z[i, j] = multi_memory_function(Ω_range[i], β_red[:, j], α_j, var_params[j][1], var_params[j][2], [i for i in phonon_mode_freqs], 0.12)
		CSV.write("multi_conductivity_data.csv", DataFrames.DataFrame([[0.0, [i for i in T_range]...]'; [Ω_range 1 ./ Z]], :auto))
	end
end

# ╔═╡ 0c424ddb-aecd-4fc6-9f26-54bb1fc9f099
Z

# ╔═╡ fff9cb8d-00aa-426c-bb66-4b7b21871c36
Plots.plot(Ω_range, real.(1 ./ Z), minorgrid = true, label = hcat(["T = $i K" for i in T_range]...), yaxis = :log)

# ╔═╡ 8ace5afe-29da-4bd1-9599-b51c03b67b6e
Plots.plot(Ω_range, imag.(1 ./ Z), minorgrid = true, label = hcat(["T = $i K" for i in T_range]...))

# ╔═╡ fab3ad29-61f6-45dd-b63f-4ea40c960afb
Plots.plot(Ω_range, abs.(1 ./ Z), minorgrid = true, label = hcat(["T = $i K" for i in T_range]...))

# ╔═╡ 3fb0f337-c6dd-4f56-a694-00d1c2f403bd
# Hellwarth Theory

# ╔═╡ 53117290-63b5-469c-93a9-bab8a40b2658
# B Scheme

# ╔═╡ a6b7c896-04d8-4704-8173-e8895ff84cbd
α_hellwarth_B = PolaronMobility.frohlichalpha(4.5, 24.1, B_scheme * 1e12, 0.12)

# ╔═╡ 6d322516-d247-4bd0-bcd5-9b8371b47750
β_red_hellwarth_B = [ħ * 2π / kB / T * B_scheme * 1e12 for T in T_range]

# ╔═╡ ae756cc9-f4ca-4073-abff-c56ecd1fdade
var_B = PolaronMobility.feynmanvw.(α_hellwarth_B, β_red_hellwarth_B)

# ╔═╡ 4c90b122-1bb9-460e-ac68-608ca8074a50
v_B = [i[1] for i in var_B]

# ╔═╡ 83068836-36d9-42aa-9867-7555d89e0d42
w_B = [i[2] for i in var_B]

# ╔═╡ 8082688f-ae0f-4aad-9bb9-71852be0fe00
Hellwarth_energy_B = F.(v_B, w_B, β_red_hellwarth_B, α_hellwarth_B) .* 1e3 .* ħ .* 2π .* B_scheme .* 1e12 ./ eV

# ╔═╡ 466354b5-7d93-4e2b-ab6d-f7a01112fbfb
begin
	B_data = DataFrames.DataFrame(
			alpha = [α_hellwarth_B for i in 1:length(T_range)],
			efffreq = [B_scheme for i in 1:length(T_range)],
			temp = T_range,
			beta = β_red_hellwarth_B,
			v = v_B,
			w = w_B,
			F = Hellwarth_energy_B
		)
	CSV.write("B_data.csv", B_data)
end

# ╔═╡ b1cf4fe2-8dda-4882-b39f-f23adb2738d1
begin
	ZB = Array{ComplexF64}(undef, length(Ω_range), length(T_range))
	for j in 1:length(T_range), i in 1:length(Ω_range)
		println("Ω: $(Ω_range[i]), T: $(T_range[j])")
		ZB[i, j] = multi_memory_function(Ω_range[i], β_red_hellwarth_B[j], [α_hellwarth_B], v_B[j], w_B[j], [B_scheme / 2π], 0.12) * 2π
		CSV.write("B_conductivity_data.csv", DataFrames.DataFrame([[0.0, [i for i in T_range]...]'; [Ω_range 1 ./ ZB]], :auto))
	end
end

# ╔═╡ a15e4d30-9ed2-492a-bfa4-0608180738c1
ZB

# ╔═╡ 2015d16c-e9f8-4072-834b-a317223b854f
Plots.plot(Ω_range, real.(1 ./ ZB), minorgrid = true, label = hcat(["T = $i K" for i in T_range]...), yaxis = :log)

# ╔═╡ 313a5d22-9a70-4c20-acb2-402c1683fbf8
Plots.plot(Ω_range, imag.(1 ./ ZB), minorgrid = true, label = hcat(["T = $i K" for i in T_range]...))

# ╔═╡ b8061ceb-1610-481b-986b-c5b82b9678db
Plots.plot(Ω_range, abs.(1 ./ ZB), minorgrid = true, label = hcat(["T = $i K" for i in T_range]...))

# ╔═╡ c7cdee43-9ac6-490b-9466-367a1f136766
# A Scheme

# ╔═╡ bb33a5c2-555b-4078-b278-90e5ae864fed
A_scheme = reverse([PolaronMobility.HellwarthAScheme(MAPI, T = i) for i in T_range])

# ╔═╡ 44604f2b-e5fd-438d-87df-0f10f874d00d
α_hellwarth_A = PolaronMobility.frohlichalpha.(4.5, 24.1, A_scheme .* 1e12, 0.12)

# ╔═╡ e858791b-e869-490f-9af3-2826f480d722
β_red_hellwarth_A = [ħ * 2π / kB / T_range[i] * A_scheme[i] * 1e12 for i in 1:length(T_range)]

# ╔═╡ 1f47a9e8-c2de-4963-ad58-aca71dbbb2ac
var_A = PolaronMobility.feynmanvw.(α_hellwarth_A, β_red_hellwarth_A)

# ╔═╡ 105665d5-ba53-4ddd-86df-dc40d5fefa8f
v_A = [i[1] for i in var_A]

# ╔═╡ 4468bf47-7668-4309-8f67-59fee2a1cf43
w_A = [i[2] for i in var_A]

# ╔═╡ 4bae3659-666b-4ae7-924a-b20d98a43148
Hellwarth_energy_A = F.(v_A, w_A, β_red_hellwarth_A, α_hellwarth_A) .* 1e3 .* ħ .* 2π .* A_scheme .* 1e12 ./ eV

# ╔═╡ 6e03cbc9-c6cd-440b-8652-fa1be4d7934a
begin
	A_data = DataFrames.DataFrame(
		alpha = α_hellwarth_A,
		efffreq = A_scheme,
		temp = T_range,
		beta = β_red_hellwarth_A,
		v = v_A,
		w = w_A,
		F = Hellwarth_energy_A
    )
	CSV.write("A_data.csv", A_data)
end

# ╔═╡ 53ac6ee6-f5d1-4923-9352-9d1cd96817fc
begin
	ZA = Array{ComplexF64}(undef, length(Ω_range), length(T_range))
	for j in 1:length(T_range), i in 1:length(Ω_range)
		println("Ω: $(Ω_range[i]), T: $(T_range[j])")
		ZA[i, j] = multi_memory_function(Ω_range[i], β_red_hellwarth_A[j], α_hellwarth_A[j], v_A[j], w_A[j], A_scheme[j] / 2π, 0.12) * 2π
		CSV.write("A_conductivity_data.csv", DataFrames.DataFrame([[0.0, [i for i in T_range]...]'; [Ω_range 1 ./ ZA]], :auto))
	end
end

# ╔═╡ 9ec88954-f8a6-4ae0-a800-3723de537d10
ZA

# ╔═╡ e129c967-858e-46cb-b6ce-29c238f45277
Plots.plot(Ω_range, real.(1 ./ ZA), label = hcat(["T = $i K" for i in T_range]...), minorgrid = true)

# ╔═╡ c8cdaa0e-05cb-4b9a-99dc-6932312ec56e
Plots.plot(Ω_range, imag.(1 ./ ZA), label = hcat(["T = $i K" for i in T_range]...), minorgrid = true)

# ╔═╡ 0f1066e0-d13a-4e1d-b9c0-c51509f25eb6
Plots.plot(Ω_range, abs.(1 ./ ZA), label = hcat(["T = $i K" for i in T_range]...), minorgrid = true)

# ╔═╡ Cell order:
# ╠═d09dec61-ab12-4424-bfaf-30f155c90ef9
# ╠═89fcb3a0-db48-11eb-30ac-b501d7fcd55f
# ╠═4a7ed3cb-88e4-4e86-90c1-285849994ab9
# ╠═8c54eedc-7c8f-4c5d-879e-62743e302801
# ╠═3edfd42b-4b8a-4080-b617-740ef3fbf76b
# ╠═97aa8c27-f2fe-418c-a662-906dbc126fb6
# ╠═d08e66f9-636b-4377-b653-389ee7b78a24
# ╠═50c1a877-807e-498e-87f9-0e909080a452
# ╠═e349bf0a-554f-46ae-8e88-8a7b80d0c4cd
# ╠═4cd925a7-048a-4696-b2b9-6a6b394e0528
# ╠═b04ccf74-853d-4ab2-886f-5b6c21cc9a51
# ╠═4e349b51-0590-4f98-b777-f13db8050b6d
# ╠═1c6fff5f-839a-43e1-acfe-040df2ed4822
# ╠═2f524fa1-c08b-4c11-b05a-e2d016ce6ec0
# ╠═8510da63-37d7-45d0-9012-8dd17745a9e8
# ╠═d9fca06f-ee9c-4790-a059-98667dcd8428
# ╠═04e16581-ca4f-4456-bea2-2e8bf44a08f8
# ╠═ba626ebc-d7da-4eca-8149-907d724a2ef5
# ╠═af26a6c8-9c97-4130-bd54-c9fe40e58879
# ╠═dd41fbf9-b4b9-46b9-abbf-96eddf59b0ad
# ╠═2975a6b8-bcec-4e4e-be94-5a3897c650db
# ╠═748d6916-6125-4f9a-adbf-cee3b25dc3ba
# ╠═a304f008-b3fa-4a92-a8c1-9ce2f7b4b3df
# ╠═5652b186-0e8d-469c-bfea-bf665f75c06a
# ╠═2e0eb8f9-3b6b-4464-b715-258c96e22422
# ╠═143b1c45-ce4b-4f3c-9674-ff181c64fe90
# ╠═3900999b-fdfa-48ae-afd4-bca5d03008e3
# ╠═30d37afa-7b70-479f-8dcc-860af723a96e
# ╠═9c140cb6-286b-4e25-a3dc-829b1eea9dfb
# ╠═27af9453-8044-4b1f-85fd-9e2c5dd104eb
# ╠═65a76ecd-1dcd-4d5c-88be-1b57d31369da
# ╠═6a793986-844a-4a82-84b0-b71a48aa8242
# ╠═146bcc43-15fc-45de-bb9b-f16a68d84e41
# ╠═3c7564b0-b9fe-49e8-9ad7-eaa80e9841a7
# ╠═0c424ddb-aecd-4fc6-9f26-54bb1fc9f099
# ╠═fff9cb8d-00aa-426c-bb66-4b7b21871c36
# ╠═8ace5afe-29da-4bd1-9599-b51c03b67b6e
# ╠═fab3ad29-61f6-45dd-b63f-4ea40c960afb
# ╠═3fb0f337-c6dd-4f56-a694-00d1c2f403bd
# ╠═53117290-63b5-469c-93a9-bab8a40b2658
# ╠═a6b7c896-04d8-4704-8173-e8895ff84cbd
# ╠═6d322516-d247-4bd0-bcd5-9b8371b47750
# ╠═ae756cc9-f4ca-4073-abff-c56ecd1fdade
# ╠═4c90b122-1bb9-460e-ac68-608ca8074a50
# ╠═83068836-36d9-42aa-9867-7555d89e0d42
# ╠═8082688f-ae0f-4aad-9bb9-71852be0fe00
# ╠═466354b5-7d93-4e2b-ab6d-f7a01112fbfb
# ╠═b1cf4fe2-8dda-4882-b39f-f23adb2738d1
# ╠═a15e4d30-9ed2-492a-bfa4-0608180738c1
# ╠═2015d16c-e9f8-4072-834b-a317223b854f
# ╠═313a5d22-9a70-4c20-acb2-402c1683fbf8
# ╠═b8061ceb-1610-481b-986b-c5b82b9678db
# ╠═c7cdee43-9ac6-490b-9466-367a1f136766
# ╠═bb33a5c2-555b-4078-b278-90e5ae864fed
# ╠═44604f2b-e5fd-438d-87df-0f10f874d00d
# ╠═e858791b-e869-490f-9af3-2826f480d722
# ╠═1f47a9e8-c2de-4963-ad58-aca71dbbb2ac
# ╠═105665d5-ba53-4ddd-86df-dc40d5fefa8f
# ╠═4468bf47-7668-4309-8f67-59fee2a1cf43
# ╠═4bae3659-666b-4ae7-924a-b20d98a43148
# ╠═6e03cbc9-c6cd-440b-8652-fa1be4d7934a
# ╠═53ac6ee6-f5d1-4923-9352-9d1cd96817fc
# ╠═9ec88954-f8a6-4ae0-a800-3723de537d10
# ╠═e129c967-858e-46cb-b6ce-29c238f45277
# ╠═c8cdaa0e-05cb-4b9a-99dc-6932312ec56e
# ╠═0f1066e0-d13a-4e1d-b9c0-c51509f25eb6


text:
using OrnsteinZernike
using Test

function oz_sw(clsr)
    rho = 0.8
    nrho = 30
    ktemp = 1.5
    p = OrnsteinZernike.Parameters(1.0, 6.0, rho, 2^11, nrho, ktemp)
    pot = OrnsteinZernike.SquareWell(1.5)
    st = OrnsteinZernike.Structure(p, pot)

    inter = OrnsteinZernike.Interaction(p, st, clsr)

    return OrnsteinZernike.solve(inter)
end

function oz_hs(br)
    rho = 0.8
    nrho = 30
    p = OrnsteinZernike.Parameters(1.0, 8.0, rho, 2^11, nrho, 1.0)
    pot = OrnsteinZernike.HardSphere()
    st = OrnsteinZernike.Structure(p, pot)
    inter = OrnsteinZernike.Interaction(p, st, br)
    result = OrnsteinZernike.solve(inter)

    return result
end

@testset "OrnsteinZernike.jl" begin
    include("physicstests.jl")
    include("potentials.jl")
    include("closures.jl")
    include("typetests.jl")
end


text:
num_node = 5
A = rand([0, 1], 5, 5)
A = A .| A'
num_edge = sum(A)
nf = rand(6, num_node)
ef = rand(7, num_edge)
gf = rand(8)

fg = FeaturedGraph(A, nf, ef, gf)

@testset "selector" begin
    fs = FeatureSelector(:node)
    @test fs(fg) == nf

    fs = FeatureSelector(:edge)
    @test fs(fg) == ef

    fs = FeatureSelector(:global)
    @test fs(fg) == gf

    @test_throws ArgumentError FeatureSelector(:foo)
end


text:
include("stream.jl")
# using Plots; pyplot()
using JLD2
using Parameters
using LaTeXStrings
using OffsetArrays
using BrownianDynamics
using PyPlot
using SparseArrays
rc("text",usetex=true)
rc("font",size=10)
rc("contour",negative_linestyle="solid")
rc("font",family="DejaVu Sans")
## Read in data.
cd(@__DIR__)
@load "$(pwd())/../steady_state_heat_engine.jld" u_steady params xx yy T0 V_equilibrium
# V_equilibrium -= minimum(V_equilibrium)
@unpack Pmat,Tmat,Jx,Jy,V,Vxshift,Vyshift,V_x,V_y,dx,dy,A,coupling = params
nx,ny = size(V_x)[1]-1,size(V_x)[2]-1
T = reshape(u_steady[(nx*ny+1:end)],nx,ny)
P = reshape(u_steady[1:nx*ny],nx,ny)
Tsol = fill(T0,0:nx+1,0:ny+1)
Tsol[1:nx,1:ny] .= reshape(u_steady[(nx*ny+1):end],nx,ny)
# Periodicity in the x direction.
Tsol[0,1:ny] .= Tsol[nx-1,1:ny]
Tsol[nx+1,1:ny] .= Tsol[2,1:ny]

nablaT_x = [(Tsol[i,j]-Tsol[i-1,j])/dx for i in 2:nx, j in 1:ny]
nablaT_y = [(Tsol[i,j]-Tsol[i,j-1])/dy for i in 2:nx, j in 1:ny]
Jx,Jy = density_current((xx,yy),u_steady,params)
Jx = rotr90(Jx) # fixed this so it is doing the right thing - MJ
Jy = -rotr90(Jy)
T = rotr90(T)
nablaT_x = -rotr90(nablaT_x)
nablaT_y = rotr90(nablaT_y)
# xx = reverse(xx)
# yy = reverse(yy)
# Plotting
## Potential.
close("all")
figure(1,figsize=(3,3))
contour(xx,yy,rotl90(V_equilibrium),linewidths=0.6)
hlines([1],0,1,linewidth=2,color="blue")
annotate("Cold bath: "*L"T=T_{C}",xy=(0.25,1.05))
hlines([-1],0,1,linewidth=2,color="red")
annotate("Hot bath: "*L"T=T_{H}",xy=(0.25,-1.18))
ylim(-1.2,1.2)
ylabel(L"y/l")
xlabel(L"x/L")
annotate("1",xy=(0.15,0.6))
annotate("2",xy=(0.5,0.0))
annotate("3",xy=(0.85,-0.7))
savefig("../figures/heat_engine_potential.pdf",bbox_inches="tight")
##
# Create a dummy contourf plot, so that we can steal its colorbar.
figure(3)
cbar_ax = contourf(V_equilibrium)
figure(1)
colorbar(cbar_ax)
# savefig("../../TwoDPaper/figures/heat_engine_potential.pdf",bbox_inches="tight")
## First do the heat engine.
## Probability current.
figure(2,figsize=(3.1,9))
plt = subplot(4,1,1)
subplots_adjust(hspace=0.22,wspace=0.08)
stride = 10
quiver(xx[3:stride:end-1],yy[3:stride:end-1],Jx[3:stride:end-1,2:stride:end-1],
       Jy[3:stride:end-1,2:stride:end-1],zorder=2,scale=6,width=0.005,headwidth=4)
ψ = flowfun(Jx[2:end-1,2:end-1],Jy[2:end-1,2:end-1])

#cs = contour(xx[2:end-1],yy[2:end-1],ψ,colors="g",
 #       levels=linspace(minimum(ψ),maximum(ψ),6),
 #       linewidths=0.8,linestyle="solid")

cs = contour(xx[2:end-1],yy[2:end-1],ψ,colors="g",
        levels=[-1.842,0.526,2.805],
       linewidths=0.8,linestyle="solid")
# clabel(cs,inline=1)
##
# close("all")
# contour(ψ,levels=[1.0,5.0,8.3])
# colorbar()
##
# Make an invisible colorbar, so that everything lines up.
cbar_ax = colorbar()
cbar_ax.ax.set_visible(false)
ylabel(L"y/l")
xticks([])
# Quiver plot.
ylabel(L"y/l")
yticks([-1,-0.5,0,0.5,1])
xticks([])
## Temperature and heat.
subplot(4,1,2)
stride = 14
quiver(xx[2:stride:end],yy[2:stride:end],-nablaT_x[2:stride:end,1:stride:end],
       -nablaT_y[2:stride:end,1:stride:end],zorder=2,scale=25,width=0.005,headwidth=4)
contour(xx,yy,T,linewidths=1,colorbar=true,zorder=2)
# Create a dummy contourf plot, so that we can steal its colorbar.
figure(3)
cbar_ax = contourf(T)
figure(2)
colorbar(cbar_ax,ticks=round.(range(1,stop=3,length=4),digits=2))
xticks([])
# yticks([-1,0])
ylabel(L"y/l")
## Now the heat pump.
@load "$(pwd())/../steady_state_heat_pump.jld" u_steady params xx yy T0 V_equilibrium
@unpack Pmat,Tmat,Jx,Jy,V,Vxshift,Vyshift,V_x,V_y,dx,dy,A,coupling = params
nx,ny = size(V_x)[1]-1,size(V_x)[2]-1
T = reshape(u_steady[(nx*ny+1:end)],nx,ny)
P = reshape(u_steady[1:nx*ny],nx,ny)
Tsol = fill(T0,0:nx+1,0:ny+1)
Tsol[1:nx,1:ny] .= reshape(u_steady[(nx*ny+1):end],nx,ny)
# Periodicity in the x direction.
Tsol[0,1:ny] .= Tsol[nx-1,1:ny]
Tsol[nx+1,1:ny] .= Tsol[2,1:ny]

nablaT_x = [(Tsol[i,j]-Tsol[i-1,j])/dx for i in 1:nx, j in 1:ny]
nablaT_y = [(Tsol[i,j]-Tsol[i,j-1])/dy for i in 1:nx, j in 1:ny]
Jx,Jy = density_current((xx,yy),u_steady,params)
Jx = rotr90(Jx)
Jy = -rotr90(Jy)
T = rotr90(T)
nablaT_x = -rotr90(nablaT_x)
nablaT_y = rotr90(nablaT_y)
# xx = reverse(xx)
# yy = reverse(yy)
# @load "$(pwd())/../steady_state_heat_pump2018-07-03.jld" u_steady params xx yy T0 V_equilibrium
# nx,ny = length(xx),length(yy)
# T = Tmat[1:nx,1:ny]
# nablaT_x = ∇Tx
# nablaT_y = ∇Ty
# T = rotr90(T)
# nablaT_x = -rotr90(nablaT_x)
# nablaT_y = rotr90(nablaT_y)
# Jx = -rotr90(Jx)
# Jy = rotr90(Jy)
## Probability current.
subplot(4,1,3)
ψ = flowfun(Jx[2:end-1,2:end-1],Jy[2:end-1,2:end-1])
contour(xx[2:end-1],yy[2:end-1],ψ,colors="g",levels=range(minimum(ψ),stop=maximum(ψ),length=6),
        linewidths=0.6,linestyle="solid")
# Make an invisible colorbar, so that everything lines up.
cbar_ax = colorbar()
cbar_ax.ax.set_visible(false)
ylabel(L"y/l")
# Quiver plot.
stride = 10
quiver(xx[3:stride:end-1],yy[3:stride:end-1],Jx[3:stride:end-1,2:stride:end-1],
    Jy[3:stride:end-1,2:stride:end-1],zorder=2,scale=0.5,width=0.004,headwidth=4)
xticks([])
# xlabel(L"x")
ylabel(L"y/l")
yticks([-1,-0.5,0,0.5,1])
## Temperature and heat.
subplot(4,1,4)
stride = 14
quiver(xx[2:stride:end],yy[2:stride:end],-nablaT_x[2:stride:end,2:stride:end],
       -nablaT_y[2:stride:end,2:stride:end],zorder=2,scale=0.5,width=0.005,headwidth=4)
contour(xx,yy,T,linewidths=1,colorbar=true,zorder=2,vmin=1.00)
# Create a dummy contourf plot, so that we can steal its colorbar.
figure(3)
# cbar_ax = contourf(T)
cbar_ax = contourf(xx,yy,T,range(1.00,stop=1.02,length=6))
figure(2)
colorbar(cbar_ax,ticks=range(1.00,stop=1.02,length=3))
# yticks([-1,0])
xlabel(L"x/L")
ylabel(L"y/l")
## Save the plot into the TeX directory.
savefig("../figures/heat_engine.pdf",bbox_inches="tight")


text:
using FileManager
using Visualization
file = raw"C:\Users\marte\Documents\GEOWEB\TEST\REGISTRATION\TEST_casaletto.las"

aabb = FileManager.las2aabb(file)

using BenchmarkTools
PC = FileManager.las2pointcloud(file) #2.681 ms

Visualization.VIEW(
    [
    Visualization.points(PC.coordinates, PC.rgbs)
    ]
)


text:
abstract StochasticProcess{F<:VariateForm, S<:ValueSupport} <: Sampleable{F,S}

typealias UnivariateStochasticProcess{S<:ValueSupport} StochasticProcess{Univariate, S}
typealias MultivariateStochasticProcess{S<:ValueSupport} StochasticProcess{Multivariate, S}

typealias DiscreteStochasticProcess{F<:VariateForm} StochasticProcess{F, Discrete}
typealias ContinuousStochasticProcess{F<:VariateForm} StochasticProcess{F, Continuous}

#typealias DiscreteUnivariateStochasticProcess StochasticProcess{Univariate, Discrete}
typealias ContinuousUnivariateStochasticProcess StochasticProcess{Univariate, Continuous}
#typealias DiscreteMultivariateStochasticProcess StochasticProcess{Multivariate, Discrete}
typealias ContinuousMultivariateStochasticProcess StochasticProcess{Multivariate, Continuous}


text:
using EngEconomics, Roots

# Given
# A = Landfill
# B = Ship
aInit = -1000000
aFinal = -100000
aLife = 30
aAnnual_Last20 = 30000
aAnnual_First10 = -20000 - aAnnual_Last20

bAnnual = -130000

MARR = 0.11

# Determine Annual Worth Method
aAW = (aInit + aAnnual_First10 * capitalRecoveryFactor(MARR, 10)) * capitalRecoveryFactor(MARR, aLife) +  aAnnual_Last20 + aLife * sinkingFundFactor(MARR, aLife)
bAW = bAnnual

if aAW >= bAW
	println("A")
else
	println("B")
end