text:
# Copyright (c) 2014: SCS.jl contributors
#
# Use of this source code is governed by an MIT-style license that can be found
# in the LICENSE.md file or at https://opensource.org/licenses/MIT.

abstract type LinearSolver end

abstract type AbstractSCSType end

Base.cconvert(::Type{Ptr{Cvoid}}, x::AbstractSCSType) = x

function Base.unsafe_convert(::Type{Ptr{Cvoid}}, x::AbstractSCSType)
    return pointer_from_objref(x)
end

mutable struct ScsMatrix{T} <: AbstractSCSType
    values::Ptr{Cdouble}
    rowval::Ptr{T}
    colptr::Ptr{T}
    m::T
    n::T
end

mutable struct ScsSettings{T} <: AbstractSCSType
    normalize::T # boolean, heuristic data rescaling
    scale::Cdouble # if normalized, rescales by this factor
    adaptive_scale::T # boolean, whether to adaptively update `scale`
    rho_x::Cdouble # x equality constraint scaling
    max_iters::T # maximum iterations to take
    eps_abs::Cdouble # absolute convergence tolerance
    eps_rel::Cdouble # relative convergence tolerance
    eps_infeas::Cdouble # infeasible convergence tolerance
    alpha::Cdouble # Douglas-Rachford relaxation parameter
    time_limit_secs::Cdouble # time limit in seconds
    verbose::T # boolean, write out progress
    warm_start::T # boolean, warm start (put initial guess in Sol struct)
    acceleration_lookback::T # memory for acceleration
    acceleration_interval::T # interval to apply acceleration
    write_data_filename::Cstring # if set dump raw problem data to the file
    log_csv_filename::Cstring # if set log solve data to the file
    function ScsSettings(linear_solver::Type{<:LinearSolver})
        settings = new{scsint_t(linear_solver)}()
        scs_set_default_settings(linear_solver, settings)
        return settings
    end
end

mutable struct ScsData{T} <: AbstractSCSType
    m::T
    n::T
    A::Ptr{ScsMatrix{T}}
    P::Ptr{Cvoid} # Ptr{ScsMatrix{T}} or C_NULL
    b::Ptr{Cdouble}
    c::Ptr{Cdouble}
end

mutable struct ScsCone{T} <: AbstractSCSType
    z::T
    l::T
    bu::Ptr{Cdouble}
    bl::Ptr{Cdouble}
    bsize::T
    q::Ptr{T}
    qsize::T
    s::Ptr{T}
    ssize::T
    ep::T
    ed::T
    p::Ptr{Cdouble}
    psize::T
end

mutable struct ScsSolution <: AbstractSCSType
    x::Ptr{Cdouble}
    y::Ptr{Cdouble}
    s::Ptr{Cdouble}
end

mutable struct ScsInfo{T} <: AbstractSCSType
    iter::T
    status::NTuple{128,Cchar}
    lin_sys_solver::NTuple{128,Cchar}
    status_val::T
    scale_updates::T
    pobj::Cdouble
    dobj::Cdouble
    res_pri::Cdouble
    res_dual::Cdouble
    gap::Cdouble
    res_infeas::Cdouble
    res_unbdd_a::Cdouble
    res_unbdd_p::Cdouble
    setup_time::Cdouble
    solve_time::Cdouble
    scale::Cdouble
    comp_slack::Cdouble
    rejected_accel_steps::T
    accepted_accel_steps::T
    lin_sys_time::Cdouble
    cone_time::Cdouble
    accel_time::Cdouble
    ScsInfo{T}() where {T} = new{T}()
end

mutable struct Solution{T}
    x::Array{Float64,1}
    y::Array{Float64,1}
    s::Array{Float64,1}
    info::ScsInfo{T}
    ret_val::T
end

"""
    _ScsDataWrapper

A type for wrapping all data inputs to SCS to prevent them from being freed by
the garbage collector during a solve.

You should not construct this manually. Call `scs_solve` instead.
"""
struct _ScsDataWrapper{S,T}
    linear_solver::S
    m::T
    n::T
    Avalues::Vector{Float64}
    Arowval::Vector{T}
    Acolptr::Vector{T}
    A::ScsMatrix{T}
    Pvalues::Vector{Float64}
    Prowval::Vector{T}
    Pcolptr::Vector{T}
    P::ScsMatrix{T}
    b::Vector{Cdouble}
    c::Vector{Cdouble}
    z::T
    l::T
    bu::Vector{Cdouble}
    bl::Vector{Cdouble}
    q::Vector{T}
    s::Vector{T}
    ep::T
    ed::T
    p::Vector{Cdouble}
    primal::Vector{Cdouble}
    dual::Vector{Cdouble}
    slack::Vector{Cdouble}
    settings::ScsSettings{T}
    options::Any
end

function _sanitize_options(options)
    option_dict = Dict{Symbol,Any}()
    for (key, value) in options
        if key == :linear_solver
            continue
        elseif !hasfield(ScsSettings, key)
            throw(
                ArgumentError(
                    "Unrecognized option passed to SCS solver: $(key)",
                ),
            )
        end
        option_dict[key] = value
    end
    return option_dict
end

function raw_status(info::ScsInfo)
    data = UInt8[info.status[i] for i in 1:findfirst(iszero, info.status)-1]
    return String(data)
end

function _to_sparse(::Type{T}, A::AbstractMatrix) where {T}
    return _to_sparse(T, SparseArrays.sparse(A))
end

function _to_sparse(::Type{T}, A::SparseArrays.SparseMatrixCSC) where {T}
    rowval = convert(Vector{T}, A.rowval .- T(1))
    colptr = convert(Vector{T}, A.colptr .- T(1))
    return A.nzval, rowval, colptr
end

"""
    scs_solve(linear_solver, args...)

The low-level interface to the SCS solver.

!!! warning
    This function is an advanced feature with a risk of incorrect usage. If you
    are a new user, we recommend that you use the JuMP or Convex interfaces
    instead.

## Problem definition

SCS solves a problem of the form
```
minimize        1/2 * x' * P * x + c' * x
subject to      A * x + s = b
                s in K
```
where K is a product cone of
- zero cone
- positive orthant `{ x | x ≥ 0 }`
- box cone `{ (t,x) | t*l ≤ x ≤ t*u }`
- second-order cone (SOC) `{ (t,x) | ||x||_2 ≤ t }`
- semi-definite cone (SDC) `{ X | X is psd }`
- exponential cone `{ (x,y,z) | y e^(x/y) ≤ z, y > 0 }`
- power cone `{ (x,y,z) | x^a * y^(1-a) ≥ |z|, x ≥ 0, y ≥ 0 }`
- dual power cone `{ (u,v,w) | (u/a)^a * (v/(1-a))^(1-a) ≥ |w|, u ≥ 0, v ≥ 0 }`

## Input arguments

The problem data are:
- `linear_solver`: a `LinearSolver` to use
- `m`: the number of affine constraints
- `n`: the number of variables
- `A`: an `AbstractMatrix` with `m` rows and `n` cols
- `P`: the upper-triangular component of a positive semidefinite symmetric
       matrix of size `(n, n)`
- `b`: a `Vector` of length `m`
- `c`: a `Vector` of length `n`

The rows of `A` correspond to cones in `K` and need to be specified in the order
above by the following arguments.

- `z`: the number of primal zero / dual free cones, i.e. primal equality
  constraints
- `l`: the number of linear cones
- `bu`: the `Vector` of upper bounds for the box cone
- `bl`: the `Vector` of lower bounds for the box cone
- `q`: the `Vector` of SOCs sizes
- `s`: the `Vector` of SDCs sizes
- `ep`: the number of primal exponential cones
- `ed`: the number of dual exponential cones
- `p`: the `Vector` of power cone parameters (±1, with negative values for the
  dual cone)

Provide a warm start to SCS by overriding:
- `primal_sol = zeros(n)`: a `Vector` to warmstart the primal variables,
- `dual_sol = zeros(m)`: a `Vector` to warmstart the dual variables,
- `slack = zeros(m)`: a `Vector` to warmstart the slack variables.

!!! note
    To successfully warmstart the solver `primal_sol`, `dual_sol` and `slack`
    must all be provided **and** `warm_start` option must be set to `true`.

Finally, a number of kwarg arguments may be passed as options to the solver.
The valid keywords are `fieldnames(ScsSettings)`. See the official SCS
documentation for information on each keyword.

## Output

This function eturns a `Solution` object, which contains the following fields:
```julia
mutable struct Solution{T}
    x::Vector{Float64}
    y::Vector{Float64}
    s::Vector{Float64}
    info::ScsInfo{T}
    ret_val::T
end
```
where `x` stores the optimal value of the primal variable, `y` stores the
optimal value of the dual variable, `s` is the slack variable, and `info`
contains various information about the solve step.

!!! warning
    SCS expects the semi-definite cones to be scaled by a factor of √2. That is,
    the off-diagonal elements in the A matrix and b vector should be multiplied
    by √2, and then the corresponding rows in the dual and slack solution
    vectors should be multiplied by 1/√2.

`SCS.raw_status(::ScsInfo)::String` describes the status, e.g. `"solved"`,
`"infeasible"`, `"unbounded"`, etc. For the precise return status, the value of
`ret_val` field should be compared with
https://github.com/cvxgrp/scs/blob/3aaa93c7aa04c7001df5e51b81f21b126dfa99b3/include/glbopts.h#L18.
"""
function scs_solve(
    linear_solver::Type{<:LinearSolver},
    m::Integer,
    n::Integer,
    A,
    P,
    b::Vector{Float64},
    c::Vector{Float64},
    z::Integer,
    l::Integer,
    bu::Vector{Float64},
    bl::Vector{Float64},
    q::Vector{<:Integer},
    s::Vector{<:Integer},
    ep::Integer,
    ed::Integer,
    p::Vector{Float64},
    primal_sol::Vector{Float64} = zeros(n),
    dual_sol::Vector{Float64} = zeros(m),
    slack::Vector{Float64} = zeros(m);
    warm_start::Bool = false,
    options...,
)
    if length(primal_sol) == n && length(dual_sol) == length(slack) == m
        if !warm_start
            fill!(primal_sol, 0.0)
            fill!(dual_sol, 0.0)
            fill!(slack, 0.0)
        end
    else
        if warm_start
            throw(ArgumentError("Warmstart doesn't match the problem size"))
        end
        primal_sol, dual_sol, slack = zeros(n), zeros(m), zeros(m)
    end
    T = scsint_t(linear_solver)
    m, n, z, l, ep, ed = T(m), T(n), T(z), T(l), T(ep), T(ed)
    Avalues, Arowval, Acolptr = _to_sparse(T, A)
    Pvalues, Prowval, Pcolptr = _to_sparse(T, P)
    option_dict = _sanitize_options(options)
    option_dict[:warm_start] = convert(Int, warm_start)
    model = _ScsDataWrapper(
        linear_solver,
        m,
        n,
        Avalues,
        Arowval,
        Acolptr,
        ScsMatrix(pointer(Avalues), pointer(Arowval), pointer(Acolptr), m, n),
        Pvalues,
        Prowval,
        Pcolptr,
        ScsMatrix(pointer(Pvalues), pointer(Prowval), pointer(Pcolptr), n, n),
        b,
        c,
        z,
        l,
        bu,
        bl,
        convert(Vector{T}, q),
        convert(Vector{T}, s),
        ep,
        ed,
        p,
        primal_sol,
        dual_sol,
        slack,
        ScsSettings(linear_solver),
        option_dict,
    )
    Base.GC.@preserve model begin
        return _unsafe_scs_solve(model)
    end
end

function _unsafe_scs_solve(model::_ScsDataWrapper{S,T}) where {S,T}
    scs_cone = ScsCone{T}(
        model.z,
        model.l,
        pointer(model.bu),
        pointer(model.bl),
        max(length(model.bu) - 1, 0),
        pointer(model.q),
        length(model.q),
        pointer(model.s),
        length(model.s),
        model.ep,
        model.ed,
        pointer(model.p),
        length(model.p),
    )
    scs_data = ScsData{T}(
        model.m,
        model.n,
        pointer_from_objref(model.A),
        iszero(model.Pvalues) ? C_NULL : pointer_from_objref(model.P),
        pointer(model.b),
        pointer(model.c),
    )
    for (key, value) in model.options
        if value isa AbstractString
            if value isa String
                cstr =
                    Base.unsafe_convert(Cstring, Base.cconvert(Cstring, value))
                setproperty!(model.settings, key, cstr)
            else
                error("You must pass a `String` as the value for $(key)")
            end
        else
            setproperty!(model.settings, key, value)
        end
    end
    scs_work = scs_init(model.linear_solver, scs_data, scs_cone, model.settings)
    scs_solution = ScsSolution(
        pointer(model.primal),
        pointer(model.dual),
        pointer(model.slack),
    )
    scs_info = ScsInfo{T}()
    status = scs_solve(
        model.linear_solver,
        scs_work,
        scs_solution,
        scs_info,
        model.options[:warm_start],
    )
    scs_finish(model.linear_solver, scs_work)
    return Solution(model.primal, model.dual, model.slack, scs_info, status)
end


text:
abstract type SamplerState end

mutable struct HMCState{Q, P, E, LL, MM, UU, KK, HH} <: SamplerState
    q::Q    # position, variables of interest
    p::P    # momentum, auxiliary variables
    ϵ::E    # step size
    L::LL   # num of steps
    M::MM   # mass matrix
    u::UU
    k::KK
    h::HH

    function HMCState(q, ϵ, L; M=nothing)
        if M === nothing
            M = I(length(q))
        end
        p = rand(MvNormal(zeros(length(q)), M))
        if length(p) == 1
            M = convert(Float64, M...)
            p = convert(Float64, p...)
        end
        u, k, h = 0.0, 0.0, 0.0
        new{typeof(q), typeof(p), typeof(ϵ), typeof(L), typeof(M), typeof(u), typeof(k), typeof(h)}(q, p, ϵ, L, M, u, k, h)
    end
end;

function leapfrog!(z::HMCState, ∂U::Function, ∂K::Function)
    ϵ = z.ϵ
    L = z.L
    q = z.q
    p = z.p

    # new state according to leapfrom for hamilton dynamics
    p -= ϵ .* ∂U(q) / 2
    for l in 1:1:L-1
        q += ϵ .* ∂K(p)
        p -= ϵ .* ∂U(q)
    end
    q += ϵ .* ∂K(p)
    p -= ϵ .* ∂U(q) / 2 # the second half for the kinetic variables

    z.q = q
    z.p = -p
end;

function hmcstep!(z::HMCState, U::Function, ∂U::Function, K::Function, ∂K::Function)

    current_q = copy(z.q)
    z.p = rand(MvNormal(zeros(length(current_q)), z.M))
    current_p = copy(z.p)

    # new state according to leapfrom for hamilton dynamics
    leapfrog!(z, ∂U, ∂K)

    # Evaluate potential and kinetic energies at start and end of trajectory
    current_U = U(current_q)
    current_K = K(current_p)
    proposed_U = U(z.q)
    proposed_K = K(z.p)

    if  rand(Float64) .< exp(current_U - proposed_U + current_K - proposed_K)
        z.q = z.q
        z.u = proposed_U
        z.k = proposed_K
        z.h = proposed_U + proposed_K
    else
        z.q = current_q
        z.u = current_U
        z.k = current_K
        z.h = current_U + current_K
    end
end;

ext: jl

max_stars_repo_name: cmerkatas/HMC

max_issues_repo_name: cmerkatas/HMC

max_forks_repo_name: cmerkatas/HMC


text:
macro _tvl_or(x, y)
    t1 = gensym()
    quote
        let $t1 = $x
            if $t1 === true
                true
            elseif $t1 === false
                $y
            else
                if $y === true
                    true
                else
                    missing
                end
            end
        end
    end
end

macro _tvl_and(x, y)
    t1 = gensym()
    quote
        let $t1 = $x
            if $t1 === false
                false
            elseif $t1 === true
                $y
            else
                if $y === false
                    false
                else
                    missing
                end
            end
        end
    end
end

replace_and_or(e::Any) = e

function replace_and_or(e::Expr)
    if e.head == :&&
        Expr(
            :macrocall,
            Symbol("@_tvl_and"),
            LineNumberNode(0, nothing),
            esc(replace_and_or(e.args[1])),
            esc(replace_and_or(e.args[2])),
        )
    elseif e.head == :||
        Expr(
            :macrocall,
            Symbol("@_tvl_or"),
            LineNumberNode(0, nothing),
            esc(replace_and_or(e.args[1])),
            esc(replace_and_or(e.args[2])),
        )
    else
        Expr(
            e.head,
            map(ex->esc(replace_and_or(ex)), e.args)...
        )
    end
end

macro tvl(e)
    replace_and_or(e)
end


text:
module AWSExceptions

using HTTP
using HTTP.MessageRequest: body_was_streamed
using JSON
using XMLDict
using XMLDict: XMLDictElement

export AWSException, ProtocolNotDefined, InvalidFileName, NoCredentials

struct ProtocolNotDefined <: Exception
    message::String
end
Base.show(io::IO, e::ProtocolNotDefined) = println(io, e.message)

struct InvalidFileName <: Exception
    message::String
end
Base.show(io::IO, e::InvalidFileName) = println(io, e.message)

struct NoCredentials <: Exception
    message::String
end
Base.show(io::IO, e::NoCredentials) = println(io, e.message)

struct AWSException <: Exception
    code::String
    message::String
    info::Union{XMLDictElement,Dict,String,Nothing}
    cause::HTTP.StatusError
    streamed_body::Union{String,Nothing}
end

function AWSException(code, message, info, cause)
    return AWSException(code, message, info, cause, nothing)
end

function Base.show(io::IO, e::AWSException)
    print(io, AWSException, ": ", e.code)
    !isempty(e.message) && print(io, " -- ", e.message)
    print(io, "\n\n", e.cause)

    # When the response is streamed then `e.cause` will not show the response body
    if e.streamed_body !== nothing
        print(io, "\n\n")
        if isempty(e.streamed_body)
            printstyled(io, "(empty body)"; bold=true)
        else
            print(io, e.streamed_body)
        end
    end

    println(io)
    return nothing
end

http_status(e::HTTP.StatusError) = e.status
content_type(e::HTTP.StatusError) = HTTP.header(e.response, "Content-Type")
is_valid_xml_string(str) = startswith(str, '<')

AWSException(e::HTTP.StatusError) = AWSException(e, String(copy(e.response.body)))

function AWSException(e::HTTP.StatusError, stream::IO)
    seekstart(stream)
    body = read(stream, String)
    return AWSException(e, body)
end

function AWSException(e::HTTP.StatusError, body::AbstractString)
    code = string(http_status(e))
    message = "AWSException"
    info = Dict{String,Dict}()

    if !isempty(body)
        # Extract API error code from Lambda-style JSON error message...
        if endswith(content_type(e), "json")
            info = JSON.parse(body)
        end

        # Extract API error code from JSON error message...
        if occursin(r"^application/x-amz-json-1\.[01]$", content_type(e))
            info = JSON.parse(body)
            if haskey(info, "__type")
                code = rsplit(info["__type"], '#'; limit=2)[end]
            end
        end

        # Extract API error code from XML error message...
        error_content_types = ["", "application/xml", "text/xml"]
        if content_type(e) in error_content_types && is_valid_xml_string(body)
            info = parse_xml(body)
        end
    elseif parse(Int, HTTP.header(e.response, "Content-Length", "0")) > 0
        # Should only occur streaming a response and error handling is improperly configured
        @error "Internal Error: provided body is empty while the reported content-length " *
            "is non-zero"
    end

    # There are times when Errors or Error are returned back
    info = get(info, "Errors", info)
    info = get(info, "Error", info)

    code = get(info, "Code", code)

    # There are also times when the response back is (M|m)essage
    message = get(info, "Message", message)
    message = get(info, "message", message)

    streamed_body = e.response.body == body_was_streamed ? body : nothing

    return AWSException(code, message, info, e, streamed_body)
end

end


text:
# This code calculates the lyapunov spectrum of the KS equation
# by transforming it in spectral space.
# It uses very efficient handling of Fourier transforms.
# It saves exponents to be used by the figure generating scripts.
using DrWatson
@quickactivate "NonlinearDynamicsTextbook"
include(srcdir("style.jl"))
using OrdinaryDiffEq
import FFTW, LinearAlgebra

# In the following `u` is the field in real space, while `y` is the Fourier
# transform of the field, which lives in inverse space.
# `w` is a matrix that has the first column the state `y` and the remaining the
# deviation vectors `δy`.

function kse_spectral_and_tangent!(dw, w, p, t)
    @unpack  ks, forward_plan, inverse_plan, yd, ud, ud2, ud3, yd2, ik2, k²_k⁴ = p
    y = view(w, :, 1)
    dy = view(dw, :, 1)

    # KS equation in spectral space
    LinearAlgebra.mul!(ud, inverse_plan, y) # create current u in x-space
    @. ud3 = ud*ud
    y² = LinearAlgebra.mul!(yd, forward_plan, ud3) # transform to k-space
    @. dy = - y²*ik2 + k²_k⁴*y

    # Equations for tangent space in spectral space (loop over each deviation vector)
    for i in 2:size(w, 2)
        δy = view(w, :, i)
        dδy = view(dw, :, i)
        ud2 .= ud # keeps track for tangent
        LinearAlgebra.mul!(ud2, inverse_plan, δy)
        ud2 .*= ud
        spectral_tan = LinearAlgebra.mul!(yd2, forward_plan, ud2)
        @. dδy = k²_k⁴*δy + im*ks*spectral_tan
    end
    return nothing
end

function orthonorm_spectral(D, M = D) # Orthonormal vectors converted to spectral space
    M > D && throw(ArgumentError("M must be ≤ D"))
    Δ = Matrix(LinearAlgebra.qr(rand(D, M)).Q)
    FFTW.rfft(Δ)
end

function tannorm(u::AbstractMatrix, t)
    s = size(u)[1]
    x = 0.0
    for i in 1:s; @inbounds x += abs2(u[i, 1]); end
    return sqrt(x/length(x))
end
tannorm(u, t) = abs(u)

# %% Lyapunov code
# This function uses some convenience methods from DynamicalSystems.jl
# based on the `tangent_integrator` format (state is a matrix. first column
# is actual system state, all other columns are deviation vectors)
import ProgressMeter
function kse_lyapunovs_spectral(integ, N, Δt::Real)
    progress = ProgressMeter.Progress(N; desc = "KSE Lyapunov Spectrum: ", dt = 1.0)
    M = size(integ.u, 2) - 1 # number of Ls
    λ = zeros(M)
    forward_plan = integ.p.forward_plan
    inverse_plan = integ.p.inverse_plan
    W = zeros(length(inverse_plan*integ.u[:, 1]), M) # for use in buffer
    t0 = integ.t
    ud = zeros(size(W, 1))
    yd = integ.u[:, 1]

    for n in 1:N
        step!(integ, Δt)
        # Get deviations in real space
        for i in 1:M
            LinearAlgebra.mul!(ud, inverse_plan, view(integ.u, :, i+1))
            W[:, i] .= ud
        end
        # Perform a (buffered) QR decomposition
        Q, R = LinearAlgebra.qr!(W)
        # Keep track of LEs
        for j in 1:M
            @inbounds λ[j] += log(abs(R[j,j]))
        end
        # Set the new deviations (also convert back to spectral)
        for i in 1:M
            ud .= Q[:, i]
            LinearAlgebra.mul!(yd, forward_plan, ud)
            integ.u[:, i+1] .= yd
        end
        u_modified!(integ, true) # Ensure that DiffEq knows we did something!
        ProgressMeter.update!(progress, n)
    end
    λ ./= (integ.t - t0)
    return λ
end

# %% Actually compute the exponents
function produce_kse_lyaps(config)
    @unpack b, dx, N, Δt, kport = config
    xs = range(0, b; step = dx) # space
    u0 = @. cos(xs) + 0.1*sin(xs/8) + 0.01*cos((2π/b)*xs)
    ks = Vector(FFTW.rfftfreq(length(u0))/dx) # conjugate space (wavenumbers)

    forward_plan = FFTW.plan_rfft(u0)
    y0 = forward_plan * u0
    inverse_plan = FFTW.plan_irfft(y0, length(u0))
    ik2 = -im .* ks ./ 2
    k²_k⁴ = @. ks^2 - ks^4

    ud = copy(u0)
    ud2 = copy(u0)
    ud3 = copy(u0)
    yd = copy(y0)
    yd2 = copy(y0)
    ksparams = @ntuple forward_plan inverse_plan ks ud ud2 ud3 yd yd2 k²_k⁴ ik2

    D = length(u0)
    M = round(Int, kport*D)
    δy0 = orthonorm_spectral(D, M)
    w0 = hcat(y0, δy0)
    prob = ODEProblem(kse_spectral_and_tangent!, w0, (0.0, 100.0), ksparams)
    integ = init(prob, Tsit5(); save_everystep = false, internalnorm = tannorm)

    @time λs = kse_lyapunovs_spectral(integ, N, Δt)
    return @strdict λs b dx N Δt kport
end

b = 20      # length of spatial domain
N = 4000    # steps to run the QR normalization
Δt = 1.0    # timestep between normalizations
dx = 0.2    # spatial discretization
kport = 0.5

for b in 20:20:60
    if b == 20
        kport = 0.5
    else
        kport = 0.3
    end
    @show b
    config = @strdict N dx Δt b kport
    file, path = produce_or_load(
        datadir("ksiva"), config, produce_kse_lyaps;
        prefix = "ksiva_spectrum", ignores = ["kport"], storepatch = false,
    )
end

b = 20

for dx in (0.2, 0.15, 0.1)
    kport = 0.5
    @show dx
    config = @strdict N dx Δt b kport
    file, path = produce_or_load(
        datadir("ksiva"), config, produce_kse_lyaps;
        prefix = "ksiva_spectrum", ignores = ["kport"], storepatch = false,
    )
end

max_stars_repo_name: JuliaDynamics/NonlinearDynamicsTextbook

max_issues_repo_name: JuliaDynamics/NonlinearDynamicsTextbook

max_forks_repo_name: JuliaDynamics/NonlinearDynamicsTextbook


text:
# Paul Rodriguez translation into Julia 0.3.5 of Monica Costa Dias and Corma O'Dea "Dynamic Economics in Practice"
# using Alex code for Getting Equiprobable Shock Points from a Normal Distributions
# March 2015, p.lesmes.11@ucl.ac.uk
# Please improve it!

# ------------------------------------------------------------------------
# DESCRIPTION
# This program solves and simulates a finite period consumption and saving
# problem. There is income that can be uncertain. The income program can be
# hardcoded by the user or can be set to follow a log normal autoregessive
# process

# IMPORTANT: Comment lines 188 and 189 if you want to change the number of
# periods or use shocks generated in Julia

# ------------------------------------------------------------------------
# PREAMBLE
# Ensure that all storage spaces variables, matrices, memory, globals are
# clean from information stored in past runs

using Distributions   # For drwaing distributions
using Optim           # Basic optimization routines
using Grid            # Interpolation tools

# Set here the relevant directory. Please check
cd("C://Users//Paul//Documents//GitHub//ucl-econ-julia//LifecycleCostaDias//v5_julia")

# ----------------------------------------------------------------------

mainDir=pwd()

# Numeric tools
cd("code_numericTools")
require("getEquiprobNormalDeviates.jl")
require("truncate.jl")
require("interp1lin.jl")
require("interp2D.jl")
cd(mainDir)

# Grids and other inputs
cd("code_gridsInputs")
require("checkInputs.jl")
require("getIncomeGrid.jl")
require("getMinandMaxAss.jl")
require("getGrid.jl")
cd(mainDir)

# Programs for the solution
cd("code_solution")
require("solveValueFunction.jl")
require("objectivefunc.jl")
require("utility.jl")
cd(mainDir)

# Programs for the simulation
cd("code_simulation")
require("simNoUncer.jl")
require("simWithUncer.jl")
cd(mainDir)

# Plot stuff
cd("code_plots")
require("plotPaths.jl")
require("plots.jl")
cd(mainDir)

tic()        # start the clock

# ------------------------------------------------------------------------
# DECLARE VARIABLES AND MATRICES THAT WILL BE 'GLOBAL'
# explicitly set variables and matrices to be shared throughout the routine
# as globals

#% ------------------------------------------------------------------------
# NUMERICAL METHODS
# select solution, interpolation and integration methods

#global interpMethod = 'pchip';      # interpolation method - IN JULIA, UP TO THE MOMENT, YOU CAN ONLY USE LINEAR WITH UNEVEN GRIDS
const linearise = 1;               # whether to linearise the slope of EV when using EE - set linearise=1 to do this, else = 0

# ------------------------------------------------------------------------
# NUMERICAL CONSTANTS
# set constants needed in numerical solution and simulation

# precision parameters
#--------------------------------------%
const tol = 1e-10;                 # max allowed error
const minCons = 1e-5;              # min allowed consumption

# where to truncate the normal distributions
#--------------------------------------%
const normBnd = 3;                 #Ignore draws less than -NormalTunc*sigma and greater than normalTrunc*sigma

# information for simulations
#--------------------------------------%
const numSims = 2;                #How many individuals to simulate

# ------------------------------------------------------------------------
# THE ECONOMIC ENVIRONMENT
# Set values of structural economic parameters

const T = 40;                      # Number of time period (If you want to reproduce matlab results, be sure that it has the same number of periods!)
const r = 0.01;                    # Interest rate
const beta = 0.98;                 # Discount factor
const gamma = 1.5;                 # Coefficient of relative risk aversion
const mu = 0;                      # mean of initial log income
const sigma = 0.25;                   # variance of log income
const rho = 0.75;                     # persistency of log income
const Tretire = 41;                # age after which there is no income earned
const borrowingAllowed = 0;        # Is borrowing allowed
const isUncertainty = 1;           # Is income uncertain?
const startA = 0;                  # How much asset do people start life with

# ------------------------------------------------------------------------
# GRIDS
# choose dimensions, set matrices and select methods to construct grids

#The grid for assets
#--------------------------------------%
const numPointsA = 20;             # number of points in the discretised asset grid
gridMethod = "3logsteps";    # method to construct grid. One of equalsteps, logsteps, 3logsteps, 5logsteps or 10logsteps

#The grid for income shocks
#--------------------------------------%
const numPointsY = 5;           #  points in grid for income (should be 2 if hard-coded)
const uncertaintyMethod = 1;    #  =0 if we enter shocks manually, =1 if put shocks on a grid using Tauchen (1986) method
hcIncome = [0.5941 ,   0.8179  ,  1.0000  ,  1.2227  ,  1.6832];  # hard-coded shocks, used if uncertaintyMethod = 1
hcIncPDF = [0.1016 ,   0.2492  ,  0.2983  ,  0.2492  ,  0.1017];  # and respective probabilities

# Check inputs
checkInputs()

## Get income grid
(Ygrid, incTransitionMrx, minInc, maxInc) = getIncomeGrid();

## ------------------------------------------------------------------------
# GET ASSET GRID
# populate grid for assets using 'gridMethod'
# populate matrix borrowingConstraints

(borrowCon, maxAss) = getMinAndMaxAss(borrowingAllowed, minInc, maxInc, startA);

Agrid = Array(Float64,T+1, numPointsA);
for ixt = 1:1:T+1
    Agrid[ixt, :] = getGrid(borrowCon[ixt], maxAss[ixt], numPointsA, gridMethod);
end

#% ------------------------------------------------------------------------
# SOLVE CONSUMER'S PROBLEM
# Get policy function and value function

(policyA1, policyC, val, exVal) = solveValueFunction();

	# Want to compare with Matlab version? In that way you can test how different are
	# simulations for the same policyFunc in Matlab and Julia
	# policyA1=reshape(readdlm("matlablObj\\policyA1matlab.csv", ',' ),40,20,5);

#% ------------------------------------------------------------------------
# SIMULATE CONSUMER'S PATHS
# start from initial level of assets and simulate optimal consumption and
# savings profiles over lifecycle

if isUncertainty == 0
    (ypath, cpath, apath, vpath) = simNoUncer(policyA1, exVal, startA);
else
	# Get the draws, and the feed the simulation with them (sligthly different than the original code)

	# For the discreate shocks of Income case... Julia makes our life way easier!
	srand(1223424) # Seed it!
	randYpath=rand( Categorical(vec(hcIncPDF))  ,(T,numSims));

	# For the normal distributed case
	srand(1223424) # Seed it!
 	e=rand(Normal(0, sigma),(T,numSims)); # normally distributed random draws for the innovation
 	srand(234636)  # Seed it!
 	sig_inc = sigma/ ((1-rho^2)^0.5);
 	logy1   =rand(Normal(mu, sig_inc),numSims);  # a random draw for the initial income

		# IMPORTANT: Want to compare with Matlab version? In that way you can test how different are
		# simulations using the same shocks. NOTICE THAT THIS OVERIDES THE GENERATED SHOCKS FROM
		# THE PREVIOUS LINES
		# Comment this lines if you want to change the number of periods
 		e      = readdlm("matlablObj\\ematlab.csv", ',' );
 		logy1  = readdlm("matlablObj\\logy1matlab.csv", ',' );

    (ypath, cpath, apath, vpath) = simWithUncer(policyA1,exVal, startA,e,logy1,randYpath);
end

toc();     # Stop the clock

#% ------------------------------------------------------------------------
# PLOTS
# Plots some features of the simulation and simulation
# Use a program that supports Gadfly, in these days, either JUNO, IJULIA

if 1==1  # In purspouse, Gadfly is just so slow!! But there are few good alternatives (March 2015)

	# Slow.........
	using Gadfly
	using DataFrames      # For making graphs or exporting data, is better as dataframes

	plotPaths()

	# Now plot value and policy functions
	whichyear = 20;
	plotNode1 = 3;
	plotNodeLast = numPointsA;

	plots();   #Plot Value Functions
end

# ------------------------------------------------------------------------
# ------------------------------------------------------------------------

max_stars_repo_name: floswald/ucl-econ-julia

max_issues_repo_name: floswald/ucl-econ-julia

max_forks_repo_name: floswald/ucl-econ-julia


text:
"""
    LRP(c::Chain, r::AbstractLRPRule)
    LRP(c::Chain, rs::AbstractVector{<:AbstractLRPRule})

Analyze model by applying Layer-Wise Relevance Propagation.

# Keyword arguments
- `skip_checks::Bool`: Skip checks whether model is compatible with LRP and contains output softmax. Default is `false`.
- `verbose::Bool`: Select whether the model checks should print a summary on failure. Default is `true`.

# References
[1] G. Montavon et al., Layer-Wise Relevance Propagation: An Overview
[2] W. Samek et al., Explaining Deep Neural Networks and Beyond: A Review of Methods and Applications
"""
struct LRP{R<:AbstractVector{<:AbstractLRPRule}} <: AbstractXAIMethod
    model::Chain
    rules::R

    # Construct LRP analyzer by manually assigning a rule to each layer
    function LRP(
        model::Chain,
        rules::AbstractVector{<:AbstractLRPRule};
        skip_checks=false,
        verbose=true,
    )
        model = flatten_model(model)
        if !skip_checks
            check_output_softmax(model)
            check_model(Val(:LRP), model; verbose=verbose)
        end
        if length(model.layers) != length(rules)
            throw(ArgumentError("Length of rules doesn't match length of Flux chain."))
        end
        return new{typeof(rules)}(model, rules)
    end
end

# Construct LRP analyzer by assigning a single rule to all layers
function LRP(model::Chain, r::AbstractLRPRule; kwargs...)
    model = flatten_model(model)
    rules = repeat([r], length(model.layers))
    return LRP(model, rules; kwargs...)
end
# Additional constructors for convenience:
LRP(model::Chain; kwargs...) = LRP(model, ZeroRule(); kwargs...)
LRPZero(model::Chain; kwargs...) = LRP(model, ZeroRule(); kwargs...)
LRPEpsilon(model::Chain; kwargs...) = LRP(model, EpsilonRule(); kwargs...)
LRPGamma(model::Chain; kwargs...) = LRP(model, GammaRule(); kwargs...)

# The call to the LRP analyzer.
function (analyzer::LRP)(
    input::AbstractArray{T}, ns::AbstractNeuronSelector; layerwise_relevances=false
) where {T}
    layers = analyzer.model.layers
    # Compute layerwise activations on forward pass through model:
    acts = [input, Flux.activations(analyzer.model, input)...]

    # Allocate array for layerwise relevances:
    rels = similar.(acts)

    # Mask output neuron
    output_indices = ns(acts[end])
    rels[end] .= zero(T)
    rels[end][output_indices] = acts[end][output_indices]

    # Backward pass through layers, applying LRP rules
    for (i, rule) in Iterators.reverse(enumerate(analyzer.rules))
        lrp!(rels[i], rule, layers[i], acts[i], rels[i + 1]) # inplace update rels[i]
    end

    return Explanation(
        first(rels),
        last(acts),
        output_indices,
        :LRP,
        ifelse(layerwise_relevances, rels, Nothing),
    )
end


text:
struct P2P2BearingReuse
  measvec::Vector{Float64}
  predvec::Vector{Float64}
  resid::Vector{Float64}
  P2P2BearingReuse() = new(zeros(2),zeros(2),zeros(2))
end

"""
    $TYPEDEF

Single dimension bearing constraint from Pose2 to Point2 variable.
"""
struct Pose2Point2Bearing{B <: IIF.SamplableBelief} <: IIF.AbstractRelativeMinimize
    bearing::B
    reuse::Vector{P2P2BearingReuse}
    Pose2Point2Bearing{B}() where B = new{B}()
    Pose2Point2Bearing{B}(x1::B) where {B <: IIF.SamplableBelief} = new{B}(x1, [P2P2BearingReuse() for i in 1:Threads.nthreads()])
end
Pose2Point2Bearing(x1::B) where {B <: IIF.SamplableBelief} = Pose2Point2Bearing{B}(x1)

function getSample(cfo::CalcFactor{<:Pose2Point2Bearing}, N::Int=1)
  return (reshape(rand(cfo.factor.bearing, N),1,N), )
end
# define the conditional probability constraint
function (cfo::CalcFactor{<:Pose2Point2Bearing})(meas, xi, lm)
  #
  reuse = cfo.factor.reuse[Threads.threadid()]
  reuse.measvec[1] = cos(meas[1] + xi[3])
  reuse.measvec[2] = sin(meas[1] + xi[3])

  # @simd for i in 1:2
    reuse.predvec[1:2] .= lm[1:2].-xi[1:2]
  # end
  normalize!(reuse.predvec)
  reuse.resid .= reuse.measvec
  reuse.resid .-= reuse.predvec

  # must return result of length zDim==1 in this case
  return sum(abs.(reuse.resid))
  # reuse.resid .^= 2
  # res[1] = reuse.resid[1]
  # res[1] += reuse.resid[2]

end

# Packing and Unpacking
mutable struct PackedPose2Point2Bearing <: IncrementalInference.PackedInferenceType
    bearstr::String
    PackedPose2Point2Bearing() = new()
    PackedPose2Point2Bearing(s1::AS) where {AS <: AbstractString} = new(string(s1))
end
function convert(::Type{PackedPose2Point2Bearing}, d::Pose2Point2Bearing{B}) where {B <: IIF.SamplableBelief}
  return PackedPose2Point2Bearing(convert(PackedSamplableBelief, d.bearing))
end
function convert(::Type{Pose2Point2Bearing}, d::PackedPose2Point2Bearing)
  Pose2Point2Bearing(convert(SamplableBelief, d.bearstr))
end


text:
using ColorBrewer
using Base.Test

### Test 1: ArgumentError thrown when the scheme is invalid.
@test_throws ArgumentError palette("NOTASCHEME", 5)

### Test 2: ArgumentError thrown when the number of categories for a scheme
#   is invalid.
@test_throws ArgumentError palette("Greens", 2)

### Test 3: Function runs correctly.
palette("Pastel1", 6)

### Test 4: Validate 9-color YlOrRd palette
palette("YlOrRd", 9)


text:
function FisherIndex(p0, p1, q0, q1)
    p1q0 = sum(p1.*q0)
    p0q0 = sum(p0.*q0)
    p1q1 = sum(p1.*q1)
    p0q1 = sum(p0.*q1)

    # fromula
    sqrt(p1q0/p0q0 * p1q1/p0q1)

end


text:
function getTarget(filename)
    values = readlines(filename*".txt")
    values = split(values[1], ':')
    values = values[2]
    values = split(values, ',')
    xstring = values[1]
    xstring = split(xstring, '=')
    xstring = xstring[2]
    xstring = split(xstring, "..")
    x1 = parse(Int, xstring[1])
    x2 = parse(Int, xstring[2])
    ystring = values[2]
    ystring = split(ystring, '=')
    ystring = ystring[2]
    ystring = split(ystring, "..")
    y1 = parse(Int, ystring[1])
    y2 = parse(Int, ystring[2])

    return (x1,x2), (y1, y2)
end

function step((x,y), (vx, vy))
    newx = x + vx
    newy = y + vy
    newvx = vx - sign(vx)
    newvy = vy - 1
    return (newx, newy), (newvx, newvy)
end

function checkTrajectory(startpos, startv, target)
    v = startv
    pos = startpos
    maxy = pos[2]
    valid = false
    while (pos[2] >= target[2][1] && pos[1] <= target[1][2])
        pos, v = step(pos, v)
        if pos[2] > maxy
            maxy = pos[2]
        end
        if pos[1] >= target[1][1] && pos[1] <= target[1][2] && pos[2] >= target[2][1] && pos[2] <= target[2][2]
            valid = true
        end
    end
    return valid, maxy
end

function solve(filename)
    xtarget, ytarget = getTarget(filename)
    pos = (0,0)
    maxy = 0
    optv = (1,1)
    valcount = 0
    for vx in -1000:1000
        for vy in -1000:1000
            v = (vx, vy)
            valid, my = checkTrajectory(pos, v, (xtarget, ytarget))
            if valid == true && my > maxy
                maxy = my
                optv = v
            end
            if valid == true
                valcount += 1
            end
        end
    end
    println(maxy)
    println(optv)
    println(valcount)
end

solve("input")
#solve("example")


text:
using HighDimPDE, Random, Test, CUDA

if CUDA.functional()
    CUDA.allowscalar(false)
    @testset "GPU MCSampling" begin
        # uniform sampling
        mc_sample = UniformSampling(-1f0,1f0)
        X = CUDA.zeros(Float32,10,8000)
        @test size(mc_sample(X)) == size(X)
        @test typeof(mc_sample(X)) == typeof(X)
        @test all( -1f0 .< mc_sample(X) .< 1f0 )

        # normal sampling
        mc_sample = NormalSampling(1f0,false)
        @test size(mc_sample(X)) == size(X)
        @test typeof(mc_sample(X)) == typeof(X)

        # normal sampling - true
        mc_sample = NormalSampling(1.,true)
        @test size(mc_sample(X)) == size(X)
        @test typeof(mc_sample(X)) == typeof(X)
    end
end

@testset "CPU MCSampling" begin
    # uniform sampling
    mc_sample = UniformSampling(0.,1.)
    X = zeros(10,8000)
    @test size(mc_sample(X)) == size(X)
    @test typeof(mc_sample(X)) == typeof(X)
    @test all( 0f0 .< mc_sample(X) .< 1f0 )

    # normal sampling
    mc_sample = NormalSampling(1.,false)
    @test size(mc_sample(X)) == size(X)
    @test typeof(mc_sample(X)) == typeof(X)

    # normal sampling - true
    mc_sample = NormalSampling(1.,true)
    @test size(mc_sample(X)) == size(X)
    @test typeof(mc_sample(X)) == typeof(X)
end


text:
using COSMO, Test, LinearAlgebra, Statistics, Random

tol = 1e-3
rng = Random.MersenneTwister(41)
mutable struct TestProblem
    P
    q
    constraints
end

function simpleQP()
    P = [4. 1;1 2]
    q = [1; 1.]
    A = [1. 1;1 0; 0 1]
    l = [1.;0;0]
    u = [1.;0.7;0.7]

    constraint1 = COSMO.Constraint(-A, u, COSMO.Nonnegatives)
    constraint2 = COSMO.Constraint(A, -l, COSMO.Nonnegatives)
    constraints = [constraint1; constraint2]
    return TestProblem(P, q, constraints)
end

@testset "Simple Tests" begin

    @testset "Simple QP" begin
        p = simpleQP()
        model = COSMO.Model()
        assemble!(model, p.P, p.q, p.constraints, settings = COSMO.Settings(verbose=true))

        res = COSMO.optimize!(model);

        @test res.status == :Solved
        @test isapprox(norm(res.x - [0.3; 0.7]), 0., atol=tol)
        @test isapprox(res.obj_val, 1.8800000298331538, atol=tol)

    end

    @testset "update_max_iter" begin
        p = simpleQP()
        settings = COSMO.Settings(max_iter=20)
        model = COSMO.Model()
        assemble!(model, p.P, p.q, p.constraints, settings = settings)

        res = COSMO.optimize!(model);
        @test res.status == :Max_iter_reached
    end

    @testset "update_check_termination" begin
         p = simpleQP()
        settings = COSMO.Settings(check_termination=100000)
        model = COSMO.Model()
        assemble!(model, p.P, p.q, p.constraints, settings = settings)

        res = COSMO.optimize!(model);

        @test res.status == :Max_iter_reached

    end

    @testset "update_rho" begin

        # TODO: write problem where rho is updated, results should be <exactly> the same

    end

    @testset "timelimit" begin
        p = simpleQP()
        settings = COSMO.Settings(time_limit=1, check_termination=100000000,max_iter=10000000)
        model = COSMO.Model()
        assemble!(model, p.P, p.q, p.constraints, settings = settings)

        res = COSMO.optimize!(model);
        @test res.status == :Time_limit_reached
    end

     @testset "warm_starting" begin
        p = simpleQP()
        settings = COSMO.Settings(check_termination = 1)
        model = COSMO.Model()
        assemble!(model, p.P, p.q, p.constraints, settings = settings)

        res1 = COSMO.optimize!(model);
        n = 2
        m = 6
        x0 = res1.x
        y0 = res1.y
        s0 = res1.s

        model = COSMO.Model()
        assemble!(model, p.P, p.q, p.constraints, settings = settings)
        COSMO.warm_start_primal!(model, x0)
        COSMO.warm_start_dual!(model, y0)
        COSMO.warm_start_slack!(model, s0)
        res2 = COSMO.optimize!(model);

        @test res1.status == :Solved && res2.status == :Solved && res2.iter < res1.iter

        @test_throws DimensionMismatch COSMO.warm_start_primal!(model, rand(4))
        @test_throws DimensionMismatch COSMO.warm_start_dual!(model, rand(2))
    end

end
nothing


text:
@testset "BasicEMD" begin
    @testset "BasicEMD Test 1" begin
        n = 200
        s = EMDSetting(n, 50, 10, 2)

        t = collect(range(0,1, length=n))
        x = zeros(n)

        for i=1:n
            x[i] = cos(22*pi*t[i]^2) + 6*t[i]^2
        end

        output = emd(x, s)

        @test size(output, 2) == s.m

        recovered = zeros(n)

        for i=1:size(output, 2)
            recovered += output[:, i]
        end

        @test recovered == x || isapprox(recovered, x)
    end

    @testset "BasicEMD Test 2" begin
        N = 1024
        s = EMDSetting(N, 50, 10, 10)

        function input_signal(x::Float64)
            omega = 2*pi/(N-1)
            return sin(17*omega*x)+0.5*(1.0-exp(-0.002*x))*sin(51*omega*x+1)
        end

        input = zeros(N)

        for i=1:N
            input[i] = input_signal(float(i))
        end

        output = emd(input, s)

        @test size(output, 2) == s.m

        recovered = zeros(N)

        for i=1:size(output, 2)
            recovered += output[:, i]
        end

        @test recovered == input || isapprox(recovered, input)
    end
end


text:
if basename(pwd()) == "aoc"
    cd("2018/15")
end

max_stars_repo_name: stellartux/advent-of-code

max_issues_repo_name: stellartux/advent-of-code

max_forks_repo_name: stellartux/advent-of-code


text:
using Random, Test, FFTW
using FourierTools
using ImageTransformations
using IndexFunArrays
using Zygote
using NDTools

Random.seed!(42)

include("fft_helpers.jl")
include("fftshift_alternatives.jl")
include("utils.jl")
include("fourier_shifting.jl")
include("fourier_shear.jl")
include("fourier_rotate.jl")
include("resampling_tests.jl")
include("convolutions.jl")
include("custom_fourier_types.jl")
include("damping.jl")
include("czt.jl")

return


text:
module SerialManipulators

# robot arm kinematics and dynamics; dh params, fk and ik, plots, robot definitions ...

end # (sub)module SerialManipulators


text:
module nfxp

    # theta is the parameter vector.
    # it contains
    # theta[1] = RC
    # theta[2] = variable cost parameter
    # theta[3:end] = transition probabilities on mileage grid

    using Logging, FreqTables
    using Plots, NLopt , Optim
    using ShowItLikeYouBuildIt, DataFrames, DataFramesMeta
    using JuMP,Ipopt

    pyplot()
    Logging.configure(level=DEBUG)

        # function P = statetransition(p, n)
        #     p=[p; (1-sum(p))];
        #     P=0;
        #     for i=0:numel(p)-1;
        #         P=P+sparse(1:n-i,1+i:n,ones(1,n-i)*p(i+1), n,n);
        #         P(n-i,n)=1-sum(p(1:i));
        #     end
        #     P=sparse(P);

    type Param
        # opts
        max_fxpiter :: Int
        min_cstp    :: Int
        max_cstp    :: Int
        ctol        :: Float64
        rtol        :: Float64
        nstep       :: Int
        ltol0       :: Float64
        tol_vfi     :: Float64
        rtolnk      :: Float64
        hessian     :: Symbol
        converged   :: Bool

        # spaces
        n :: Int  # num of grid points
        max :: Int # upper boudn of mileage grid

        # structural params
        pr :: Vector{Float64}  # transition probs
        RC :: Float64         # replacement cost
        c :: Float64  # cost parameter
        beta :: Float64 # discount factor

        # computation objects
        grid :: LinSpace{Float64}
        P :: SparseMatrixCSC{Float64,Int64}
        EV0 :: Vector{Float64}   # old value function
        EV1 :: Vector{Float64}   # new
        pk :: Vector{Float64}  # probability of keeping engine
        cost :: Vector{Float64}   #operating cost

        """
            create default param type
        """
        function Param(;max_fxpiter = 5,
                        min_cstp = 4,
                        max_cstp = 2000000,
                        ctol = 0.01,
                        rtol = 0.02,
                        nstep = 20,
                        ltol0 = 1.0e-10,
                        tol_vfi = 1.0e-6,
                        rtolnk = 0.5,
                        hessian = :bhhh)
            this = new()
            this.max_fxpiter=max_fxpiter
            this.min_cstp=min_cstp
            this.max_cstp =max_cstp
            this.ctol  =ctol
            this.rtol  =rtol
            this.nstep =nstep
            this.ltol0 =ltol0
            this.tol_vfi =tol_vfi
            this.rtolnk=rtolnk
            this.hessian =hessian
            this.converged = false

            this.n = 175
            this.max = 450
            this.grid = 0.0:(this.n-1)

            this.pr = [0.0937; 0.4475; 0.4459; 0.0127]
            this.RC = 11.7257
            this.c = 2.45569
            this.beta = 0.9999
            this.P = build_trans(this.pr,this.n)
            this.cost = 0.001*this.c*collect(this.grid)
            this.EV0 = zeros(this.n)
            this.EV1 = zeros(this.n)
            this.pk = zeros(this.n)
            return this
        end
    end
    function update!(p::Param,theta::Vector{Float64})
        p.RC = theta[1]
        p.c = theta[2]
        if length(theta) > 2
            p.pr = theta[3:end]
            p.P = build_trans(p.pr,p.n)
        end
        p.cost = 0.001*p.c*collect(p.grid)
        p.EV0[:] = 0.0
        p.EV1[:] = 0.0
        p.pk[:] = 0.0
        p.converged = false
        # p.theta = vcat(p.RC,p.c,p.pr)
    end

    """
        build_trans(p::Vector{Float64},n::Int)

    Build transition matrix for each of `n` states where each element `p[j]` of `p` is the probability of moving `j-1` states up in the grid defined by `1:n`.
    """
    function build_trans(p::Vector{Float64},n::Int)
        P = spzeros(n,n)
        p = vcat(p,(1-sum(p)))
        N = length(p)
        for i=0:(N-1)
            P =   P + sparse(1:(n-i),(1+i):n,ones(n-i)*p[i+1],n,n)
            P[n-i,n] = 1-sum(p[1:i])
        end
        return P
    end

    function read_busses(p::Param,num_types=4)
        # load data
        d = readtable(joinpath(dirname(@__FILE__),"buses.csv"),header=false)
        # setup data
        rename!(d,[:x1,:x2,:x5,:x7,:x9],[:id,:bus_type,:d1,:x,:dx1])

        # subset to bus type
        d = @where(d,:bus_type .<= num_types)

        # discretize odometer data
        d[:x] = ceil(Int,p.n/(p.max*1000) * d[:x])
        d[:dx1] = d[:x] .- vcat(0,d[:x][1:(end-1)])
        d[:dx1] = d[:dx1].*(1-d[:d1]) .+d[:d1].*d[:x]  # replace first diff of x by x if replaced

        # get replacement dummy
        d[:d] = vcat(d[:d1][2:end],0)

        # remove obs with missing lagged mileage
        # this is the first row for each obs
        remove = d[:id] .- vcat(0,d[:id][1:(end-1)])
        d = d[remove.==0,:]

        d = @select(d,:id,:bus_type,:d,:x,:dx1)

        # add a time index
        d = @by(d,:id, t = 1:length(:d),bus_type=:bus_type,d=:d,x=:x,dx1= :dx1 )

        return d
    end

    function bellman!(p::Param,probs=false)
        vk = -p.cost .+ p.beta .* p.EV0  # value of keep at each state
        vr = -p.RC - p.cost[1] + p.beta*p.EV0[1]
        maxv = maximum(vk)
        #recenter expected value function by substracting max value
        p.EV1[:] = p.P * (maxv +  log(exp(vk-maxv) + exp(vr-maxv)))

        # compute choice probs
        if probs
            p.pk[:] = 1.0 ./ (1.0+ (exp(vr-vk)))
        end
    end

    # get Frechet derivative of Bellman operator
    function dbellman(p::Param)
        tmp = p.P[:,2:p.n] .* repmat(p.pk[2:p.n,1]',p.n,1)
        p.beta * hcat((1.0 .- sum(tmp,2)), tmp)
    end
    function dbellman2(p::Param)
        dGamma_dEV = p.beta * (p.P .* repmat(p.pk,1,p.n))
        # println(size(dGamma_dEV))
        dGamma_dEV[:,1] = dGamma_dEV[:,1] .+ p.beta*p.P*(1.0-p.pk)
        # println(size(dGamma_dEV))
        return dGamma_dEV
    end

    update_ev!(p::Param) = copy!(p.EV0,p.EV1);

    """
        VFI(p::Param)

    Standard value function iteration
    """
    function VFI!(p::Param)
        iter = 0
        for j in 1:p.max_cstp  # do max_cstp contraction steps
            iter += 1
            bellman!(p)
            dist_vfi = maxabs(p.EV0.-p.EV1)
            update_ev!(p)
            # p.EV0[:] = p.EV1[:]
            if dist_vfi < p.tol_vfi
                p.converged = true
                info("VFI converged after $iter iterations")
                break
            end
        end
        if maxabs(p.EV0.-p.EV1) > p.tol_vfi
            info("not converged")
        end
        bellman!(p,true)    # get choice probs
        # plot(p.pk)
    end

    function likelihood!(new_theta::Vector{Float64},grad::Vector{Float64},data::DataFrame,p::Param)

        update!(p,new_theta)

        n_c = length(p.c)
        N = length(data[:x])

        # solve model
        VFI!(p)
        fval = 1e10

        # try

            # evaluate log likelihood of replacement decision
            # at each state
            ccps = p.pk[data[:x]]  # get ccp of keep at each state in the data
            log_like = log( ccps .+ (1.0-2.0*ccps) .* (data[:d]))

            # add log likelihood for exogenous mileage increases
            n_p = 0
            pr = Float64[]
            if length(new_theta) > 2
                pr = [p.pr ; 1.0 - sum(p.pr)]
                n_p = length(pr)-1
                # if any(pr[1+data[:dx1]] .< 0)
                #     # println(sum(pr[1+data[:dx1]] .< 0))
                # end
                # idx1 = clamp(data[:dx1],1,maximum(data[:dx1]))
                log_like[:] = log_like[:] .+ log(pr[1+data[:dx1]])
            end
            # println("n_p=$n_p")

            # ojbective function value
            fval = mean(-log_like)

            # save all of this

            # compute gradient of likelihood function

            N = size(data[:x],1)

            if length(grad)>0
                #1. derivative of bellman equation wrt parameters
                dc=0.001*p.grid;    # deriv of cost function
                dtdmp = zeros(p.n,1+n_c+n_p)  # gradient matrix of bellman operator
                dtdmp[:,1] = p.P * p.pk - 1   # wrt RC
                dtdmp[:,2:(1+n_c)] = -(p.P * dc) .*p.pk  # wrt c
                if length(new_theta) > 2   # wrt pr
                    vk = -p.cost + p.beta*p.EV1
                    vr = -p.RC - p.cost[1] + p.beta*p.EV1[1]
                    vmax = max(vk,vr)
                    dtp = vmax + log(exp(vk-vmax)+exp(vr-vmax))

                    for iP in 1:n_p
                        dtdmp[1:(p.n-iP),1+n_c+iP] = dtp[iP:(end-1)] .- vcat(dtp[(n_p+1):p.n],Float64[dtp[end] for i in 1:(n_p-iP)])
                    end
                    invp = exp(-log(pr))
                    invp = vcat(sparse(1:n_p,1:n_p,invp[1:n_p],n_p,n_p),
                               -ones(1,n_p)*invp[n_p+1])
                end
                #2. Derivative of ev wrt parameters:
                F=speye(p.n) - dbellman2(p)
                devdmp = F\dtdmp
                # print("devdmp=$(size(devdmp))")
                # print("devdmp=$(size(devdmp[data[:x],:]))")

                # step 3 get derivative of loglike wrt parameters
                score = (ccps-1 .+ data[:d] ) .* (hcat(-ones(N), dc[data[:x]], zeros(N,n_p) ) .+ (devdmp[ones(Int,N),:] .- devdmp[data[:x],:]))
                if length(new_theta) > 2   # wrt pr
                    # idx1 = clamp(data[:dx1],1,maximum(data[:dx1]))
                    for iP in 1:n_p
                        score[:,1+n_c+iP] = score[:,1+n_c+iP] + invp[1+data[:dx1],iP]
                    end
                end
                grad[:] = mean(-score,1)
            end

        # catch
        #     fval = 1e9
        # end

        # we don't do the hessian.
        # enough is enough. :-)

        info("fval=$fval")

        return fval

    end

    function simdata(N::Int,T::Int,p::Param)
        if !p.converged
            info("solving model by VFI")
            VFI!(p)
        end
        id = reshape(repeat(1:N,inner=1,outer=T),N,T)
        t  = reshape(repeat(1:T,inner=N,outer=1),N,T)
        shock_dx = rand(N,T)
        shock_d = rand(N,T)

        x = zeros(Int,N,T)
        x1 = zeros(Int,N,T)
        x[:,1] = ones(Int,N)
        d = zeros(Bool,N,T)

        # uncontrolled transition on odometer x
        csum_p = cumsum(p.pr) / maximum(cumsum(p.pr))
        dx1 = zeros(Int,N,T)
        for i in 1:length(csum_p)
            dx1 += (shock_dx.>csum_p[i])   # increases are random: notice that you can get dx1 = 0 here. need to add +1 to dx1 to get a valid array index.
        end

        for it in 1:T
            for i in 1:N
                if x[i,it] == 0
                    println(x[i,:])
                    println(d[i,:])
                end
                d[i,it] = shock_d[i,it] < (1.0-p.pk[x[i,it]])
                # make odometer progress
                x1[i,it] = min((1-d[i,it])*x[i,it] + d[i,it] + dx1[i,it] , p.n)   # this adds +1 for example
                if it < T
                    x[i,it+1] = x1[i,it]
                end
            end

        end
        return DataFrame(id = id[:],t = t[:], d = convert(Array{Int},d[:]), x = x[:],x1 = x1[:], dx1 = dx1[:])
    end

    function run_estim()
        p = Param()
        d = read_busses(p,4)  # select all bus types smaller than this number
        e = estimate(d,p)
        return e
    end

    function run_MC()
        p = Param()
        d = simdata(50,113,p)
        e = estimate(d,p)
        return e
    end

    function likelihood_single()
        p = Param()
        d = simdata(50,113,p)
        pr = freqtable(d[:dx1][d[:dx1].>0]) ./ sum(freqtable(d[:dx1][d[:dx1].>0]).array)
        # pr=probs[1:(end-1),:]
        update!(p,vcat(0,0,pr.array...))
        g = zeros(length(vcat(0,0,pr.array...)))
        f = nfxp.likelihood!(vcat(0,0,pr.array...),g,d,p)
    end
    function simulate_single_run()
        p = Param()
        d = simdata(50,113,p)
    end

    function estimate(d::DataFrame,p::Param,startv=zeros(2))
        # step 1: PML for transition probs
        probs = freqtable(d[:dx1][d[:dx1].>0]) ./ sum(freqtable(d[:dx1][d[:dx1].>0]).array)
        pr=probs[1:(end-1),:]

        info("starting values for probs = $(p.pr)")
        info("updated with $pr")
        update!(p,vcat(startv,pr.array...))

        # step 2: estimate structural params
        # f_closure(x,g) = likelihood!(x,g,d,p)
        f_closure(x) = likelihood!(x,Float64[],d,p)
        # opt = Opt(:LD_MMA,length(startv))
        opt = Opt(:LN_COBYLA,length(startv))
        # lower_bounds!(opt,[-Inf,-Inf,[0.0 for i in 1:4]...])
        # upper_bounds!(opt,[Inf,Inf,[1.0 for i in 1:4]...])
        min_objective!(opt,f_closure)
        xtol_rel!(opt,1e-4)
        ftol_rel!(opt,1e-6)
        maxeval!(opt,50)

        res = Optim.optimize(f_closure,startv)

        # inequality_constraint!(opt,myconstraint)

        # (minf,minx,ret) = NLopt.optimize(opt, startv )
        return res
    end

    function MC(p::Param;N=50,T=113)
        d = simdata(N,T,p)

        # estimate pr transition probs

        # put those values into the param type
        update!(p,[p.RC;p.c;probs.array])

        # maximize likelihood
    end

    function myconstraint(x::Vector,gr::Vector)
        if length(gr)>0
            gr[1] = 0
            gr[2] = 0
            gr[3:end] = 1.0
        end
        sum(x[3:end]) - 1.0 # < 0
    end

    function check_like()
        p = Param()
        d = simdata(10000,50,p)
        # pg = [Float64[x,p.theta[2:end]...] for x in linspace(1.0,20,20)]
        y = zeros(20)
        for i in eachindex(pg)
            y[i] = likelihood!(pg[i],zeros(6),d,p)
        end
        return y

    end

    function max_nlopt()
        p = Param()
        d = simdata(10000,13,p)
        # d = simdata(50,113,p)
        f_closure(x,g) = likelihood!(x,g,d,p)
        # opt = Opt(:LN_COBYLA,p.n_params)
        lower_bounds!(opt,[-Inf,-Inf,[0.0 for i in 1:4]...])
        upper_bounds!(opt,[Inf,Inf,[1.0 for i in 1:4]...])
        min_objective!(opt,f_closure)

        inequality_constraint!(opt,myconstraint)

        (minf,minx,ret) = NLopt.optimize(opt, [10.0,0.0001,0.1,0.1,0.1,0.1] )
        # (minf,minx,ret) = NLopt.optimize(opt, p.theta )
    end
    function max_optim()
        p = Param()
        d = simdata(50,113,p)
        f_closure(x) = likelihood!(x,Float64[],d,p)
        lower = [-Inf,-Inf,[0.0 for i in 1:4]...]
        upper = [Inf,Inf,[1.0 for i in 1:4]...]
        res = Optim.optimize(OnceDifferentiable(f_closure),[10.0,0.0001,0.1,0.1,0.1,0.1],lower,upper,Fminbox())
        return res
    end

    function dict_busses(df::DataFrame)
        N,T = (maximum(df[:id]),maximum(df[:t]))
        d = Dict()
        d[:d] = reshape(df[:d].data,N,T)
        d[:x] = reshape(df[:x].data,N,T)
        d[:dx1] = reshape(df[:dx1].data,N,T)
        return d
    end

    function mpec()
        p=Param()
        d = simdata(50,113,p)
        # d = read_busses(p)
        dd = dict_busses(d)
        m = Model(solver=IpoptSolver())
        T = maximum(d[:t])
        n = maximum(d[:id])
        M = 4 # number of slots in state transition: can move 0,1,2,3 slots up
        N = p.n  # number of states
        @variable(m,theta_cost >= 0)
        @variable(m,RC >= 0)
        @variable(m,theta_probs[1:M] >= 0)
        @variable(m,EV[1:N])
        cost = 0.001*theta_cost*collect(p.grid)

        # transformations

        @variable(m,CbEV[1:N])
        @constraint(m,constr_cbev[i=1:N], CbEV[i] == -cost[i] + p.beta*EV[i])
        @variable(m,PayoffDiff[1:N])
        @constraint(m,constr_payoff[i=1:N], PayoffDiff[i]== -CbEV[i] - RC + CbEV[1])
        @NLexpression(m,exp_probkeep[i=1:N],1/(1+exp(PayoffDiff[i])))
        # @variable(m,aux_probkeep[i]==[1:N])
        # @constraint(m,constr_aux_probkeep[i=1:N],aux_probkeep[i]==1.0 / (1.0 + exp_payoffdiff[i]))
        @variable(m,ProbKeep[1:N])
        @NLconstraint(m,const_keep[i=1:N],ProbKeep[i] == exp_probkeep[i])

        # BellmanViolation = sum()

        @NLobjective(m, Max,sum(log( dd[:d][i,it]*(1.0 - ProbKeep[dd[:x][i,it]]) + (1.0 - dd[:d][i,it])*(ProbKeep[dd[:x][i,it]])) for i=1:n, it=2:T) +
            sum( log( theta_probs[dd[:dx1][i,it]+1] ) for i=1:n,it=2:T))

        # bellman equation for states 1:(N-M+1) i.e. where all state progressions are possible
        @NLconstraint(m, constr_EV[i=1:(N-M+1)],
            EV[i] == sum(log(exp(CbEV[i+j]) + exp(-RC + CbEV[1])) * theta_probs[j+1] for j in 0:(M-1))
            )
        # bellman equation for states (N-M+1):(N-1) i.e. where not all state progressions are possible
        @NLconstraint(m, constr_EV_M[i=(N-M+2):(N-1)],
            EV[i] == sum(log(exp(CbEV[i+j]) + exp(-RC + CbEV[1])) * theta_probs[j+1] for j in 0:(N-i-1)) + (1- sum(theta_probs[k+1] for k in 0:(N-i-1))) * log(exp(CbEV[N])   + exp(-RC + CbEV[1])))

        # bellman equation for final state
        @NLconstraint(m, constr_EV_N,
            EV[N] == log(exp(CbEV[N]) + exp(-RC + CbEV[1])))

        # probabilities have to sum to 1
        @constraint(m,sum(theta_probs) == 1)
        # bound value function
        @constraint(m,constr_EV[i=1:N],EV[i] <= 50)

        solve(m)

    end

end  # module

max_stars_repo_name: snowdj/CoursePack

max_issues_repo_name: snowdj/CoursePack

max_forks_repo_name: snowdj/CoursePack


text:
#=

  Euler #48 in Julia.

  Problem 48:
  """
  The series, 1^(1) + 2^(2) + 3^(3) + ... + 10^(10) = 10405071317.

  Find the last ten digits of the series,
  1^(1) + 2^(2) + 3^(3) + ... + 1000^(1000).
  """

  This Julia program was created by Hakan Kjellerstrand, hakank@gmail.com
  See also my Julia page: http://www.hakank.org/julia/

=#

include("Euler.jl")

# 0.00162826s
function euler48a()
    sum = 0
    t = 10_000_000_000
    for i = 1:1_000
        n = i
        for j in 2:i
            n = (n * i) % t
        end
        sum = (sum + n) % t
    end

    return sum
end

run_euler(euler48a)

max_stars_repo_name: Wikunia/hakank

max_issues_repo_name: Wikunia/hakank

max_forks_repo_name: Wikunia/hakank


text:
module PCA

    using LinearAlgebra

    function findPCAs(
        X    :: AbstractArray{T, 2};
        num  :: Integer = 1
    ) where T <: AbstractFloat

        dim, N = size(X)

        Σ = Symmetric(X * X')

        println("Dimension of Σ : ", size(Σ))
        println("Any missing data? ", any(isnan.(Σ)))

        if any(isnan.(Σ))
            throw(ErrorException("Data contains NaN!"))
        end

        @time F = eigen(Σ)

        #println(F.values)
        #println(length(F.values))

        order = collect(1:dim)
        sort!(order; by = i -> F.values[i], rev=true)

        evs = zeros(T, dim, num)
        for i = 1:num
            evs[:, i] = normalize(F.vectors[:, order[i]])
        end

        return evs

    end

end

max_stars_repo_name: meteorologytoday/CESM-diagnostic

max_issues_repo_name: meteorologytoday/CESM-diagnostic

max_forks_repo_name: meteorologytoday/CESM-diagnostic


text:
# Julia wrapper for header: /usr/include/libavutil/spherical.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

export
    av_spherical_alloc,
    av_spherical_tile_bounds,
    av_spherical_projection_name,
    av_spherical_from_name

function av_spherical_alloc(size)
    ccall((:av_spherical_alloc, libavutil), Ptr{AVSphericalMapping}, (Ptr{Csize_t},), size)
end

function av_spherical_tile_bounds(map, width::Csize_t, height::Csize_t, left, top, right, bottom)
    ccall((:av_spherical_tile_bounds, libavutil), Cvoid, (Ptr{AVSphericalMapping}, Csize_t, Csize_t, Ptr{Csize_t}, Ptr{Csize_t}, Ptr{Csize_t}, Ptr{Csize_t}), map, width, height, left, top, right, bottom)
end

function av_spherical_projection_name(projection::AVSphericalProjection)
    ccall((:av_spherical_projection_name, libavutil), Cstring, (AVSphericalProjection,), projection)
end

function av_spherical_from_name(name)
    ccall((:av_spherical_from_name, libavutil), Cint, (Cstring,), name)
end


text:
############################
### google challenge    ####
### Author F. Clautiaux ####
############################

# Notes
# Julia uses indices from 1 to size while google data are indexed from 0 to size-1

function main(instanceFilename, assignmentFilename, verbose)

    #verbose = false # true for debug

    println("Received instance files ", instanceFilename, " and ", assignmentFilename)

    # reading the data
    data = parserGoogle(instanceFilename, assignmentFilename, verbose)

    if verbose
        printData(data)
    end

    # solve the problem
    solution = solveGoogle(data,verbose)

    if verbose # print the solution
        println("Solution of cost ", solution.cost)
        for i in 1:data.nbProcess
            print(solution.assignment[i], " ") # here the solution is written with our convention 1..n
        end
        println()
    end

    # test if the solution is feasible and compute its cost
    checkAndComputeCost(data, solution, true) # check if the solution is feasible
end


text:
using InducingPoints
using Distances
using KernelFunctions
using KernelFunctions: ColVecs
using Test
using Random: seed!
include("test_utils.jl")

@testset "InducingPoints.jl" begin
    @testset "Offline" begin
        for file in readdir(joinpath(@__DIR__, "offline"))
            include(joinpath(@__DIR__, "offline", file))
        end
    end
    @testset "Online" begin
        for file in readdir(joinpath(@__DIR__, "online"))
            include(joinpath(@__DIR__, "online", file))
        end
    end
end


text:
# This file is part of Intersections.jl

export bk

"""
    bk(onmatch::Function, L, findpos::Function=doublingsearch)
    bk(L, output=eltype(L[1])[], findpos::Function=doublingsearch)

Computes intersection of posting lists in `L` using `findpos` storing the intersection in `output`.

The method accepting `onmatch` calls `onmatch(L, P)` where `L` are the input lists (in some functions this could be a permutation of the input), and `P` and array of
indices where each match occurred in `L[i]`.

"""
function bk(onmatch::Function, L, findpos::Function=doublingsearch)
    P = ones(Int, length(L))
    n = length(L)
    _max = _get_key(L[1], 1)
    c = 0

    while true
        @inbounds for i in eachindex(P)
            P[i] = findpos(L[i], _max, P[i])
            P[i] > length(L[i]) && return
            pval = _get_key(L[i], P[i])
            if pval == _max
                c += 1
                if c == n
                    onmatch(L, P)
                    c = 0
                    P[i] += 1
                    P[i] > length(L[i]) && return
                    _max = _get_key(L[i], P[i])
                end
            else
                c = 0
                _max = pval
            end
        end
    end
end

function bk(L, output=eltype(L[1])[], findpos::Function=doublingsearch)
    bk(L, findpos) do L_, P
        push!(output, _get_key(L_[1], P[1]))
    end

    output
end


text:
"""
Convert a Franklin html string into a html string (i.e. replace `{{ ... }}`
blocks).
"""
function convert_html(hs::AS)::String
    isempty(hs) && return hs

    (:convert_html, "input hs: '$hs'") |> logger
    (:convert_html, "hasmath/code: $(locvar.((:hasmath, :hascode)))") |> logger

    # Tokenize
    tokens = find_tokens(hs, HTML_TOKENS, HTML_1C_TOKENS)

    # Find hblocks ({{ ... }})
    hblocks, tokens = find_all_ocblocks(tokens, HTML_OCB)
    deactivate_inner_blocks!(hblocks, (:COMMENT, :SCRIPT))
    filter!(hb -> hb.name ∉ (:COMMENT, :SCRIPT), hblocks)

    # Find qblocks (qualify the hblocks)
    qblocks = qualify_html_hblocks(hblocks)

    fhs = process_html_qblocks(hs, qblocks)

    (:convert_html, "fhs (pre fix): '$fhs'") |> logger

    # See issue #204, basically not all markdown links are processed  as
    # per common mark with the JuliaMarkdown, so this is a patch that kind
    # of does
    if locvar(:reflinks)::Bool
        fhs = find_and_fix_md_links(fhs)
    end
    isempty(fhs) && return ""
    return String(fhs)
end

"""
$SIGNATURES

Return the HTML corresponding to a Franklin-Markdown string as well as all the
page variables. See also [`fd2html`](@ref) which only returns the html.
"""
function fd2html_v(st::AS; internal::Bool=false,
                   dir::String="", nop::Bool=false)::Tuple{String,Dict}
    isempty(st) && return (st, LittleDict())
    if !internal
        FD_ENV[:SOURCE] = "input string"
        empty!(ALL_PAGE_VARS)
        FOLDER_PATH[] = isempty(dir) ? mktempdir() : dir
        set_paths!()
        def_GLOBAL_LXDEFS!()
        def_GLOBAL_VARS!()
        FD_ENV[:CUR_PATH] = "index.md"
    end
    # corner case if `serve` is used and cleanup has emptied global vars
    !(@isdefined GLOBAL_VARS) || isempty(GLOBAL_VARS) && def_GLOBAL_VARS!()
    m = convert_md(st; isinternal=internal)
    h = convert_html(m)
    if nop && startswith(h, "
") && endswith(h, r"
\n?")
        # remove initial
 and final
\n
        # note that this is unsafe in that these may not be related.
        h = replace(h, r"(^
|
\n?$)" => "")
    end
    return h, LOCAL_VARS
end

"""
    fd2html(s; kw)

Take a Markdown string and return the HTML that Franklin would produce for it.

# Keywords

* `internal=false`: if set to true, the current scope (in terms of page variables) will
                    be passed to process the String, otherwise the string will be processed
                    in an isolated scope.
* `dir=""`: if given and if `internal=false`, sets up a temporary dir for Franklin
            processing.
* `nop=false`: if set to true, check if the returned HTML string has a starting `<p>` and
               ending `</p>` and, if so, strips them off. This should be used with care:
               for simple strings this is likely to be ok but for strings which would
               contain multiple paragraphs (e.g. with headers) this may not work well.
"""
fd2html(s; kw...)::String = fd2html_v(s; kw...)[1]

"""
    fd2text(s; kw)

Take a Markdown string and return a version with no markup (used in RSS circumstances).
This is bruteforce to say the least, it just discards all HTML tags. It should be safe
in that '<' and '>' should have been escaped by the fd2html pass.
"""
fd2text(s; kw...)::String = replace(fd2html(s; kw...), r"\<.*?\>" => "")

"""
$SIGNATURES

Take a qualified html block stack and go through it, with recursive calling.
"""
function process_html_qblocks(hs::AS, qblocks::Vector{AbstractBlock},
                              head::Int=1, tail::Int=lastindex(hs))::String
    htmls = IOBuffer()
    head  = head # (sub)string index
    i     = 1    # qualified block index
    while i ≤ length(qblocks)
        β = qblocks[i]
        # write what's before the block
        fromβ = from(β)
        (head < fromβ) && write(htmls, subs(hs, head, prevind(hs, fromβ)))

        if β isa HTML_OPEN_COND
            content, head, i = process_html_cond(hs, qblocks, i)
            write(htmls, content)
        elseif β isa HFor
            content, head, i = process_html_for(hs, qblocks, i)
            write(htmls, content)
        # should not see an HEnd by itself --> error
        elseif β isa HEnd
            throw(HTMLBlockError("I found a lonely {{end}}."))
        # it's a function block, process it
        else
            write(htmls, convert_html_fblock(β))
            head = nextind(hs, to(β))
        end
        i += 1
    end
    # write whatever is left after the last block
    head ≤ tail && write(htmls, subs(hs, head, tail))
    return String(take!(htmls))
end

"""
    match_url(base, cand)

Try to match two url indicators.
"""
function match_url(base::AS, cand::AS)
    sbase = base[1] == '/' ? base[2:end] : base
    scand = cand[1] == '/' ? cand[2:end] : cand
    # joker-style syntax
    if endswith(scand, "/*")
        return startswith(sbase, scand[1:prevind(scand, lastindex(scand), 2)])
    elseif endswith(scand, "/")
        scand = scand[1:prevind(scand, lastindex(scand))]
    end
    return splitext(scand)[1] == sbase
end


text:
module Qsosed

# Write your package code here.
include("constants.jl")
include("types.jl")
include("parameters.jl")
include("black_hole.jl")
include("disk.jl")
include("corona.jl")
include("warm.jl")
include("xspec_routines.jl")
include("compton.jl")
include("model.jl")
include("fractions.jl")

end


text:
export GraphShortestPathEnv

using Random
using SparseArrays
using LinearAlgebra

mutable struct GraphShortestPathEnv{G, R} <: AbstractEnv
    graph::G
    pos::Int
    goal::Int
    max_steps::Int
    rng::R
    reward::Int
    step::Int
end

"""
    GraphShortestPathEnv([rng]; n=10, sparsity=0.1, max_steps=10)

Quoted **A.3** in the the paper [Decision Transformer: Reinforcement Learning via Sequence Modeling](https://arxiv.org/abs/2106.01345).

> We give details of the illustrative example discussed in the introduction.
> The task is to find theshortest path on a fixed directed graph, which can be
> formulated as an MDP where reward is0whenthe agent is at the goal node
> and−1otherwise.  The observation is the integer index of the graphnode the
> agent is in. The action is the integer index of the graph node to move to
> next. The transitiondynamics transport the agent to the action’s node index if
> there is an edge in the graph, while theagent remains at the past node
> otherwise. The returns-to-go in this problem correspond to negativepath
> lengths and maximizing them corresponds to generating shortest paths.

"""
function GraphShortestPathEnv(rng=Random.GLOBAL_RNG; n=20, sparsity=0.1, max_steps=10)
    graph = sprand(rng, Bool, n, n, sparsity) .| I(n)

    goal = rand(rng, 1:n)
    pos = rand(rng, 1:n)
    while pos == goal
        pos = rand(rng, 1:n)
    end
    GraphShortestPathEnv(graph, pos, goal, max_steps, rng, 0, 0)
end

function (env::GraphShortestPathEnv)(action)
    env.step += 1
    if env.graph[action, env.pos]
        env.pos = action
    end
    env.reward = env.pos == env.goal ? 0 : -1
end

RLBase.state(env::GraphShortestPathEnv) = env.pos
RLBase.state_space(env::GraphShortestPathEnv) = axes(env.graph, 2)
RLBase.action_space(env::GraphShortestPathEnv) = axes(env.graph, 2)
RLBase.legal_action_space(env::GraphShortestPathEnv) = (env.graph[:, env.pos]).nzind
RLBase.reward(env::GraphShortestPathEnv) = env.reward
RLBase.is_terminated(env::GraphShortestPathEnv) = env.pos == env.goal || env.step >= env.max_steps

function RLBase.reset!(env::GraphShortestPathEnv)
    env.step = 0
    env.reward = 0
    env.goal = rand(env.rng, state_space(env))
    env.pos = rand(env.rng, state_space(env))
    while env.pos == env.goal
        env.pos = rand(env.rng, state_space(env))
    end
end

Random.seed!(env::GraphShortestPathEnv, seed) = Random.seed!(env.rng, seed)

function floyd_shortest_path(env::GraphShortestPathEnv)
    n = size(env.graph, 1)
    M = fill(Inf, n, n)
    for idx in CartesianIndices(M)
        if idx[1] == idx[2]
            M[idx] = 0
        elseif env.graph[idx]
            M[idx] = 1
        end
    end
    for k in 1:n
        for i in 1:n
            for j in 1:n
                if M[i, j] > M[i, k] + M[k, j]
                    M[i, j] = M[i, k] + M[k, j]
                end
            end
        end
    end
    M
end

#=

using ReinforcementLearning
using Random

rng = MersenneTwister(123)

env = GraphShortestPathEnv(rng)
policy = RandomPolicy(rng=rng)
M = RLEnvs.floyd_shortest_path(env)

Base.@kwdef struct ShortestPathCount <: AbstractHook
    shortest_paths::Vector{Float64} = []
end

(h::ShortestPathCount)(::PreEpisodeStage, policy, env) = push!(h.shortest_paths, M[env.goal, env.pos])

h = run(policy, env, StopAfterEpisode(1_000), ComposedHook(StepsPerEpisode(), ShortestPathCount()))

using UnicodePlots

barplot(x, [sum(h[2].shortest_paths .== i) for i in x])  # shortest path

#       ┌                                        ┐
#   1.0 ┤■■■■■■■■■■■ 113
#   2.0 ┤■■■■■■■■■■■■■ 134
#   3.0 ┤■■■■■■■■■■■■■■ 147
#   4.0 ┤■■■■■■■■■■■ 111
#   5.0 ┤■■■■■■■ 76
#   6.0 ┤■■■■ 40
#   7.0 ┤■■ 19
#   8.0 ┤ 3
#   9.0 ┤ 0
#   Inf ┤■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 357
#       └                                        ┘
#

barplot(1:10, [sum(h[1].steps .== i) for i in 1:10])  # random walk

#      ┌                                        ┐
#    1 ┤■■ 39
#    2 ┤■■ 38
#    3 ┤■ 17
#    4 ┤■ 23
#    5 ┤■ 28
#    6 ┤■ 24
#    7 ┤■ 25
#    8 ┤■ 21
#    9 ┤■ 16
#   10 ┤■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 769
#      └                                        ┘
#
=#


text:
@testset "Laplace" begin
    test_interface(LaplaceLikelihood(3.0), Laplace)
    @test LaplaceLikelihood().β == 1
    test_auglik(LaplaceLikelihood(1.0); rng=MersenneTwister(42))
    # Test the custom kl divergence
    λ = rand()
    μ = rand()
    @test kldivergence(InverseGaussian(μ, 2λ), InverseGamma(1//2, λ)) ≈
        (log(2λ) / 2 - log(2π) / 2 - log(λ) / 2 + loggamma(1//2) + λ / μ)
end


text:
# Licensed to the Apache Software Foundation (ASF) under one
# or more contributor license agreements.  See the NOTICE file
# distributed with this work for additional information
# regarding copyright ownership.  The ASF licenses this file
# to you under the Apache License, Version 2.0 (the
# "License"); you may not use this file except in compliance
# with the License.  You may obtain a copy of the License at
#
#   http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing,
# software distributed under the License is distributed on an
# "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied.  See the License for the
# specific language governing permissions and limitations
# under the License.

###############################################################################
#  Types
###############################################################################

"""
    AbstractOptimizer

Base type for all optimizers.
"""
abstract type AbstractOptimizer end

"""
    AbstractLearningRateScheduler

Base type for all learning rate scheduler.
"""
abstract type AbstractLearningRateScheduler end

"""
    AbstractMomentumScheduler

Base type for all momentum scheduler.
"""
abstract type AbstractMomentumScheduler end

"""
    OptimizationState

### Attributes
* `batch_size`: The size of the mini-batch used in stochastic training.
* `curr_epoch`:
  The current epoch count. Epoch 0 means no training yet, during the first
  pass through the data, the epoch will be 1; during the second pass, the
  epoch count will be 1, and so on.
* `curr_batch`:
  The current mini-batch count. The batch count is reset during every epoch.
  The batch count 0 means the beginning of each epoch, with no mini-batch
  seen yet. During the first mini-batch, the mini-batch count will be 1.
* `curr_iter`:
  The current iteration count. One iteration corresponds to one mini-batch,
  but unlike the mini-batch count, the iteration count does **not** reset
  in each epoch. So it track the *total* number of mini-batches seen so far.
"""
mutable struct OptimizationState
  batch_size :: Int
  curr_epoch :: Int
  curr_batch :: Int
  curr_iter  :: Int
end

OptimizationState(batch_size::Int) = OptimizationState(batch_size, 0, 0, 0)

###############################################################################
#  LearningRate module
###############################################################################

module LearningRate

using Markdown

import Base: get
import ..mx: AbstractLearningRateScheduler, OptimizationState, update!

export initlrsched

initlrsched(η::Real) = LearningRate.Fixed(η)

update!(a::AbstractLearningRateScheduler) = (isdefined(a, :t) && (a.t += 1))

"""
    get(sched::AbstractLearningRateScheduler)

Returns the current learning rate.
"""
get(::AbstractLearningRateScheduler) = nothing

"""
    LearningRate.Fixed(η)

Fixed learning rate scheduler always return the same learning rate.
"""
struct Fixed <: AbstractLearningRateScheduler
  η::Float64
end

get(f::Fixed) = f.η

@doc doc"""
    LearningRate.Exp(η₀; γ = 0.9)

```math
\eta_t = \eta_0\gamma^t
```

Where `t` is the epoch count, or the iteration count.
"""
mutable struct Exp <: AbstractLearningRateScheduler
  η₀::Float64
  γ ::Float64
  t ::Int
end

function Exp(η₀; γ = 0.9, t = 0)
  @assert 0 < γ < 1
  Exp(η₀, γ, t)
end

get(a::Exp) = a.η₀ * a.γ^a.t

@doc doc"""
    LearningRate.Inv(η₀; γ = 0.9, p = 0.5)

```math
\eta_t = \eta_0 (1 + \gamma t)^{-p}
```

Where `t` is the epoch count, or the iteration count.
"""
mutable struct Inv <: AbstractLearningRateScheduler
  η₀::Float64
  γ ::Float64
  p ::Float64
  t ::Int
end

function Inv(η₀; γ = 0.9, p = 0.5, t = 0)
  @assert 0 < γ < 1
  @assert 0 <= p
  Inv(η₀, γ, p, t)
end

get(i::Inv) = i.η₀ * (1 + i.γ*i.t)^(-i.p)

end  # module LearningRate

using .LearningRate

###############################################################################
#  Momentum module
###############################################################################

module Momentum

using Markdown

import Base: get
import ..mx: AbstractMomentumScheduler, OptimizationState

export initmomsched

"""
    get(sched)

* `sched::AbstractMomentumScheduler`: the momentum scheduler.

Returns the current momentum.
"""
get

initmomsched(μ::Real) = iszero(μ) ? Momentum.Null() : Momentum.Fixed(μ)

"""
    Momentum.Null

The null momentum scheduler always returns 0 for momentum. It is also used to
explicitly indicate momentum should not be used.
"""
struct Null <: AbstractMomentumScheduler
end

get(::Null) = 0.0

"""
    Momentum.Fixed

Fixed momentum scheduler always returns the same value.
"""
mutable struct Fixed <: AbstractMomentumScheduler
  μ::Float64
end

get(f::Fixed) = f.μ

@doc doc"""
    NadamScheduler(; μ = 0.99, δ = 0.004, γ = 0.5, α = 0.96)

Nesterov-accelerated adaptive momentum scheduler.

Description in [Incorporating Nesterov Momentum into Adam]
(http://cs229.stanford.edu/proj2015/054_report.pdf).

```math
\mu_t = \mu_0 * (1 - \gamma * \alpha^{t * \delta})
```

Where
* `t`: iteration count
* `μ`: default `0.99`, μ₀
* `δ`: default `0.004` is scheduler decay.
* `γ`: default `0.5`
* `α`: default `0.96`
"""
struct NadamScheduler <: AbstractMomentumScheduler
  μ::Float64
  δ::Float64
  γ::Float64
  α::Float64
end

function NadamScheduler(; μ = 0.99, δ = 0.004, γ = 0.5, α = 0.96)
  @assert 0.0 <= μ < 1.0
  @assert 0.0 <= δ
  @assert 0.0 <= γ <= 1.0
  @assert 0.0 <= α <= 1.0
  NadamScheduler(μ, δ, γ, α)
end

"""
    get(n::NadamScheduler, t)

Where `t` is the iteration count.
"""
get(n::NadamScheduler, t) =
  n.μ * (1.0 - n.γ * n.α^( t      * n.δ)),
  n.μ * (1.0 - n.γ * n.α^((t + 1) * n.δ))

end  # module Momentum

using .Momentum

###############################################################################
# Public APIs
###############################################################################

"""
    getupdater(optimizer)

A utility function to create an updater function of `KVStore`,
that uses its closure to store all the states needed for each weights.

Ther returned function has following signature:

```julia
decend!(index::Int, ∇::NDArray, x::NDArray)
```

If the optimizer is stateful and need access/store states during updating,
`index` will be the key to access/store states.
"""
function getupdater(optimizer::AbstractOptimizer)
  states = Dict{Int,Any}()
  function updater(index::Int, ∇::NDArray, x::NDArray)
    if !haskey(states, index)
      states[index] = create_state(optimizer, index, x)
    end
    update!(optimizer, index, x, ∇, states[index])
  end
  updater
end

"""
    normgrad(optimizer, W, ∇)

Get the properly normalized gradient (re-scaled and clipped if necessary).

* `optimizer`: the optimizer,
  should contain the field `scale`, `clip` and `λ`.
* `W::NDArray`: the trainable weights.
* `∇::NDArray`: the original gradient of the weights.
"""
function normgrad!(opt::AbstractOptimizer, W::NDArray, ∇::NDArray)
  # rescaling
  s = opt.scale
  !iszero(s) && @inplace ∇ .*= s
  # gradient clipping
  c = opt.clip
  c > 0 && clamp!(∇, -c, c)
  # weight decay
  λ = opt.λ
  λ > 0 && @inplace ∇ += λ .* W

  ∇
end

###############################################################################
# Builtin Optimizers
###############################################################################

include("optimizers/sgd.jl")
include("optimizers/adam.jl")
include("optimizers/adagrad.jl")
include("optimizers/adadelta.jl")
include("optimizers/adamax.jl")
include("optimizers/rmsprop.jl")
include("optimizers/nadam.jl")


text:
using Test, MLToolkit.MonteCarlo

@testset "Monte Carlo" begin
    tests = [
        "russian_roulette.jl",
    ]
    foreach(include, tests)

    @testset "estimator_stats" begin
        n_mc = 10_000
        estimates = zeros(n_mc)
        times = zeros(n_mc)
        ground = 0.0
        for i = 1:10_000
            t = @elapsed e = randn()
            estimates[i] = e
            times[i] = t
        end
        stats = estimator_stats(estimates; times=times, ground=ground)
    end
end


text:
# julia 1.2
push!(LOAD_PATH,pwd())

using PyPlot
# import src
include("src.jl")

# init a model
MC = (
    L = exp.([-0.4, 0.4]), # labor endow states
    P = [0.8 0.2 ; 0.2 0.8],  #transition matrix
    π = [0.5,0.5],  # stationary distribution
)
M = src.ModelBasic( N = 1.0,
BOTTOM_VAL = -6.66E36, AGG_SEED = 20200417,
AGG_SIMU_SIZE = 10000, AGG_DROP_FIRST_SIZE = 1000,
Na = 200, alim = (0.0,4.5),
lgrid = MC.L, P = MC.P, π = MC.π,
θ = 0.36, δ = 0.08, β = 0.96, μ = 2.0 )

# search ss
@time src.solve_ss_gs!(M, 1.0 , max_iter_ss = 500, max_iter_hh = 500,
rtol_ss = 1E-3, rtol_hh = 1E-3, step_size_ss = 0.7, step_size_hh = 0.7,
print_err_ss = true, print_err_hh = false, if_golden = true )

# visualization
plot(M.agrid, M.HHPolicy.aOpt[:,1], linewidth = 0.75)
plot(M.agrid, M.HHPolicy.aOpt[:,2], "--", linewidth = 0.75)
plot(M.agrid, M.agrid, "-.", linewidth = 0.75)
legend(["l - low","l - high", L"$\pi/4$"])


text:
function _density(data; xlims = (-Inf, Inf), trim = false, kwargs...)
    k = kde(data; kwargs...)
    x, y = k.x, k.density
    xmin, xmax = xlims
    xmin = max(xmin, minimum(data))
    xmax = min(xmax, maximum(data))
    if trim
        for i in eachindex(x, y)
            xmin ≤ x[i] ≤ xmax || (y[i] = NaN)
        end
    end
    return spec(:Lines) * style(x, y)
end

function _density(datax, datay; xlims = (-Inf, Inf), ylims = (-Inf, Inf), trim = false, kwargs...)
    k = kde((datax, datay); kwargs...)
    x, y, z = k.x, k.y, k.density
    xmin, xmax = xlims
    xmin = max(xmin, minimum(datax))
    xmax = min(xmax, maximum(datax))
    ymin, ymax = ylims
    ymin = max(ymin, minimum(datay))
    ymax = min(ymax, maximum(datay))
    if trim
        for i in eachindex(x, y)
            xmin ≤ x[i] ≤ xmax && ymin ≤ y[i] ≤ ymax || (z[i] = NaN)
        end
    end
    return spec(:Heatmap) * style(x, y, z)
end

const density = Analysis(_density)


text:
using NLsolve, BenchmarkTools

# number of sectors
numSectors = 3

# test initial capital
k = ones(numSectors)

# test initial output statespace
y = ones(numSectors)

# linear production function parameter
ϕ = 0.5
# time discount factor
β = 0.9
# elasticity parameter
ρ = -10
# depreciation
δ = 1
# savings scaling parameters
σ = ones(numSectors, numSectors)

# utility function of consumption vector
u(c) = sum(log.(c))
# derivative of utility w.r.t consumption in one element
u′(c) = 1 / c

# production function of capital
f(k) = ϕ * k
# derivative of production w.r.t capital
f′(k) = ϕ

# future capital from current parameters
kp(s,k) = (1 - δ) * k + s

# this would be replaced by our interpolators in each sector in the iterations
# a g function of the current below form is equivalent to our "first guess" of an interpolator
g(y) = ones(2*numSectors-1)

# consumption as function of S_ij ∀ i,j
c(S,y) = y .- sum(S, 2)

# S_ij as function of S39
function S(S39)
    S = ones(numSectors, numSectors)

    # fills in the known values first (the diagonal and end collumn)
    for i in 1:numSectors
        for j in 1:numSectors
            if i == j
                S[i,j] = S39[numSectors - 1 + i]
            elseif j == numSectors
                S[i,j] = S39[i]
            end
        end
    end
    # make partial Sb from the know values of S above
    Sb = S ./ σ

    # fills in the empty values in S using our partial Sb
    for i in 1:numSectors
        for j in 1:numSectors
            if i != j && j != numSectors
                S[i,j] = sigma[i,j] * Sb[jj] * Sb[i,numSectors] / Sb[j,numSectors]
            end
        end
    end

    # make full Sbb
    for i in 1:numSectors
        for j in 1:numSectors
            Sbb[i,j] = Sb[i,j] / Sb[j,j]
        end
    end

    # return savings matrix
    return S
end

Sb(S) = S ./ σ

# s_j from full S_ij
s(S) = (sum(σ[l,:] .* S[l,:].^ρ for l in 1:numSectors)).^(1/ρ)

# NLsolve approach
function f!(F, S39)
    for i in 1:2*numSectors-1
        if i < numSectors
            # Nth (last) collumn of double blackboard S matrix, excluding diagonal corner
            F[i] = u′(c(S, y)[numSectors]) / u′(c(S, y)[i]) - Sb(S(S39))[i,numSectors]/Sb(S(S39))[numSectors,numSectors]
        else
            # Diagonals of S matrix
            j = i - numSectors + 1
            F[i] = (f′(β) * β * u′(c(S(g(y)), kp(s(S), k))[j])) / (u′(c(S, y)))^(-rho/(1 - rho)) - S39[j] / (σ[j,j] * s(S))
        end
    end
end

# warm start for savings values
initial_x = 1.5*ones(2*numSectors-1)

# solve
nlsolve(f!, initial_x)

max_stars_repo_name: uq-aibe/aim

max_issues_repo_name: uq-aibe/aim

max_forks_repo_name: uq-aibe/aim


text:
using Test, BookBayes

X = Variable(:x, 2);
Y = Variable(:y, 2);
Z = Variable(:z, 2);
ϕxyz = Factor([X, Y, Z], FactorTable(
    (x=1, y=1, z=1) => 0.08,
    (x=1, y=1, z=2) => 0.31,
    (x=1, y=2, z=1) => 0.09,
    (x=1, y=2, z=2) => 0.37,
    (x=2, y=1, z=1) => 0.01,
    (x=2, y=1, z=2) => 0.05,
    (x=2, y=2, z=1) => 0.02,
    (x=2, y=2, z=2) => 0.07,
));

ϕxz_m = marginalize(ϕxyz, :y)
ϕx_m = marginalize(ϕxz_m, :z)
ϕyz_m = marginalize(ϕxyz, :x)
ϕy_m = marginalize(ϕyz_m, :z)
ϕz_m = marginalize(ϕyz_m, :y)

ϕyz = ϕy_m * ϕz_m

ϕxz_y = condition(ϕxyz, :y, 2)
normalize!(ϕxz_y)

ψxz_y = infer(ExactInference(), )
ϕy = Factor([Y], FactorTable(
    (y=1,) => 0.20,
    (y=2,) => 0.80,
));
ϕz = Factor([Z], FactorTable(
    (z=1,) => 0.40,
    (z=2,) => 0.60,
));
g = SimpleDiGraph(3);
add_edge!(g,1,3);
add_edge!(g,2,3);
bn = BayesianNetwork([Y, Z, X], [ϕy, ϕz, ϕ], g);
assignment = (x=1, y=2, z=1);
probability(bn, assignment)

k = assignments(ϕ.vars);
@show k

@show k[2,1,1]
@show k[2,1,1].x
@show k[2,1,1].y
@show k[2,1,1].z

display(ϕ.table);
display(ϕ.vars);
display(variablenames(ϕ));

A = Variable(:A, 3);
ϕA = Factor([A], FactorTable(
    (a=1,) => 0.40,
    (a=2,) => 0.50,
    (a=3,) => 0.1,
));

ϕAy = ϕA * ϕy
ϕyz = ϕy * ϕz;
ϕAyz = ϕAy * ϕyz

[v for (k,v) in ϕAyz.table]
sum([v for (k,v) in ϕAyz.table])

sum([v for (k,v) in ϕAy.table])
sum([v for (k,v) in ϕyz.table])

@test true


text:
using Distributions
using SplineHazard
using SplineHazard.Spline

function ph_data(n=200, beta=0.0, tmax=36.0)
    z = rand(Binomial(1, 0.5), n) .- 0.5

    ## 25 patients per month
    max_entry = n / 25

    entry = rand(Uniform(0, max_entry), n)
    U = -log.(1 .-rand(Uniform(0, 1), n))

    ## 12 months median survival time
    T = U.*exp.(-beta*z)./(log(2)/12)

    ## approx: 10% admin cens, 10% drop-out
    C = min.(rand(Uniform(0, 3*36), n), tmax)

    time = min.(T, C)
    status = T .< C
    id = 1:n

    (id=id, entry=entry, time=time, status=status, z=z)
end

function ph_data_llog(n=200, beta=0.0, tmax=36.0)
    z = rand(Binomial(1, 0.5), n) .- 0.5

    ## 25 patients per month
    max_entry = n / 25

    entry = rand(Uniform(0, max_entry), n)
    U = 1 .- rand(Uniform(0, 1), n)
    T = Vector{Float64}(undef, n)
    for i in 1:n
        ##T[i] = exp(quantile(Logistic(log(12), 1/6), 1 - exp(log(U[i])*exp(-beta*z[i]))))
    end

    C = min.(rand(Uniform(0, 3*36), n), tmax)

    time = min.(T, C)
    status = T .< C
    id = 1:n

    (id=id, entry=entry, time=time, status=status, z=z)
end

function llog_hazard(time, a, b)
    d = Logistic(a, b)
    DataFrame(time=time, hazard=pdf.(d, time) ./ (1 .- cdf.(d, time)))
end

function ph_loglik(s, beta, data)
    h = hazard(data.time[data.status], s)
    ch = cumhaz(data.time, s)

    sum(log.(h) .+ beta.*data.z[data.status]) - sum(ch.*exp.(beta.*data.z))
end

function sample_beta(beta::Float64, s::Sampler, data, prior::UnivariateDistribution)
    X = data.z
    time = data.time
    status = data.status

    ch = extract_cumhaz(s, time, s.t)

    function log_target(x::Float64)
        logpdf(prior, x) + sum(status.*x.*X - ch.*exp.(x.*X))
    end

    prop = rand(Normal(beta, 1.0), 1)[1]
    r = log_target(prop) - log_target(beta)

    if log(rand()) < r
        return prop, true
    else
        return beta, false
    end
end

function ph(data, M=10000)
    s = setup_sampler(M, data.time[data.status], (sp, beta) -> ph_loglik(sp, beta, data))

    beta = Vector{Float64}(undef, M+1)
    beta[1] = 0.0
    beta_ac = 0

    prior_beta = Normal(0, 10)

    for t in 1:M
        beta[t+1], ac = sample_beta(beta[t], s, data, prior_beta)
        beta_ac += ac

        update!(s, beta[t+1])
    end

    println("Acceptance rate for beta: $(beta_ac/M)")
    println("Posterior mean: $(mean(beta[(s.warmup+1):end]))")
    println("Posterior variance: $(var(beta[(s.warmup+1):end]))")

    ##summary(s, data.time)

    beta[(s.warmup+1):end], beta_ac, s
end

function ph_example()
    data = ph_data(200, log(0.67))
    ph(data, 10000)
end