text:
# Uniforms are OpenGL variables that stay the same for the entirety of a drawcall.
# There are a lot of functions, to upload them, as OpenGL doesn't rely on multiple dispatch.
# here is my approach, to handle all of the uniforms with one function, namely gluniform
# For uniforms, the Vector and Matrix types from ImmutableArrays should be used, as they map the relation almost 1:1

GLSL_COMPATIBLE_NUMBER_TYPES = (GLfloat, GLint, GLuint, GLdouble)
const NATIVE_TYPES = Union{FixedArray, GLSL_COMPATIBLE_NUMBER_TYPES..., GLBuffer, GPUArray, Shader, GLProgram, NativeMesh}

opengl_prefix(T)  = error("Object $T is not a supported uniform element type")
opengl_postfix(T) = error("Object $T is not a supported uniform element type")

opengl_prefix{T <: Union{FixedPoint, Float32, Float16}}(x::Type{T})  = ""
opengl_prefix{T <: Float64}(x::Type{T})                     = "d"
opengl_prefix(x::Type{Cint})                                = "i"
opengl_prefix{T <: Union{Cuint, UInt8, UInt16}}(x::Type{T}) = "u"

opengl_postfix(x::Type{Float64}) = "dv"
opengl_postfix(x::Type{Float32}) = "fv"
opengl_postfix(x::Type{Cint})    = "iv"
opengl_postfix(x::Type{Cuint})   = "uiv"

function uniformfunc(typ::DataType, dims::Tuple{Int})
    Symbol(string("glUniform", first(dims), opengl_postfix(typ)))
end
function uniformfunc(typ::DataType, dims::Tuple{Int, Int})
    M, N = dims
    Symbol(string("glUniformMatrix", M==N ? "$M":"$(M)x$(N)", opengl_postfix(typ)))
end

function gluniform{FSA <: Union{FixedArray, Colorant}}(location::Integer, x::FSA)
    x = [x]
    gluniform(location, x)
end

Base.size(p::Colorant) = (length(p),)
Base.size{T <: Colorant}(p::Type{T}) = (length(p),)
Base.ndims{T <: Colorant}(p::Type{T}) = 1

@generated function gluniform{FSA <: Union{FixedArray, Colorant}}(location::Integer, x::Vector{FSA})
    func = uniformfunc(eltype(FSA), size(FSA))
    if ndims(FSA) == 2
        :($func(location, length(x), GL_FALSE, pointer(x)))
    else
        :($func(location, length(x), pointer(x)))
    end
end

#Some additional uniform functions, not related to Imutable Arrays
gluniform(location::Integer, target::Integer, t::Texture)   = gluniform(GLint(location), GLint(target), t)
gluniform(location::Integer, target::Integer, t::GPUVector) = gluniform(GLint(location), GLint(target), t.buffer)
gluniform(location::Integer, target::Integer, t::Signal)    = gluniform(GLint(location), GLint(target), t.value)
gluniform(location::Integer, target::Integer, t::TextureBuffer) = gluniform(GLint(location), GLint(target), t.texture)
function gluniform(location::GLint, target::GLint, t::Texture)
    activeTarget = GL_TEXTURE0 + UInt32(target)
    glActiveTexture(activeTarget)
    glBindTexture(t.texturetype, t.id)
    gluniform(location, target)
end
gluniform(location::Integer, x::Enum) 	 						     = gluniform(GLint(location), GLint(x))
gluniform(location::Integer, x::Signal)                              = gluniform(GLint(location), value(x))
gluniform(location::Integer, x::Union{GLubyte, GLushort, GLuint}) 	 = glUniform1ui(GLint(location), x)
gluniform(location::Integer, x::Union{GLbyte, GLshort, GLint, Bool}) = glUniform1i(GLint(location),  x)
gluniform(location::Integer, x::GLfloat)                             = glUniform1f(GLint(location),  x)
gluniform(location::Integer, x::GLdouble)                            = glUniform1d(GLint(location),  x)

#Uniform upload functions for julia arrays...
gluniform(location::GLint, x::Vector{Float32}) = glUniform1fv(location,  length(x), pointer(x))
gluniform(location::GLint, x::Vector{GLdouble}) = glUniform1dv(location,  length(x), pointer(x))
gluniform(location::GLint, x::Vector{GLint})   = glUniform1iv(location,  length(x), pointer(x))
gluniform(location::GLint, x::Vector{GLuint})  = glUniform1uiv(location, length(x), pointer(x))

glsl_typename{T}(x::T)           = glsl_typename(T)
glsl_typename(t::Type{Void})     = "Nothing"
glsl_typename(t::Type{GLfloat})  = "float"
glsl_typename(t::Type{GLdouble}) = "double"
glsl_typename(t::Type{GLuint})   = "uint"
glsl_typename(t::Type{GLint})    = "int"
glsl_typename{T<:Union{FixedVector,Colorant}}(t::Type{T}) = string(opengl_prefix(eltype(t)), "vec", length(t))
glsl_typename{T}(t::Type{TextureBuffer{T}}) = "$(opengl_prefix(eltype(T)))samplerBuffer"

function glsl_typename{T, D}(t::Texture{T, D})
    str = string(opengl_prefix(eltype(T)), "sampler", D, "D")
    t.texturetype == GL_TEXTURE_2D_ARRAY && (str *= "Array")
    str
end
function glsl_typename{T<:FixedMatrix}(t::Type{T})
    M,N = size(t)
    string(opengl_prefix(eltype(t)),"mat", M==N ? "$M" : "$(M)x$(N)")
end
toglsltype_string(t::Signal) = toglsltype_string(t.value)
toglsltype_string{T<:Union{Real, FixedArray, Texture, Colorant, TextureBuffer, Void}}(x::T) = "uniform $(glsl_typename(x))"
#Handle GLSL structs, which need to be addressed via single fields
function toglsltype_string{T}(x::T)
    if isa_gl_struct(x)
        "uniform $(T.name.name)"
    else
        error("can't splice $T into an OpenGL shader. Make sure all fields are of a concrete type and isbits(FieldType)-->true")
    end
end
toglsltype_string{T}(t::Union{GLBuffer{T}, GPUVector{T}}) = "in $(glsl_typename(T))"
# Gets used to access a
function glsl_variable_access{T,D}(keystring, t::Texture{T, D})
    t.texturetype == GL_TEXTURE_BUFFER && return "texelFetch($(keystring), index)."*"rgba"[1:length(T)]*";"
    return "getindex($(keystring), index)."*"rgba"[1:length(T)]*";"
end
glsl_variable_access(keystring, ::Union{Real, GLBuffer, GPUVector, FixedArray, Colorant}) = keystring*";"
glsl_variable_access(keystring, s::Signal) = glsl_variable_access(keystring, s.value)
glsl_variable_access(keystring, t::Any)    = error("no glsl variable calculation available for : ", keystring, " of type ", typeof(t))

function uniform_name_type(program::GLuint)
    uniformLength = glGetProgramiv(program, GL_ACTIVE_UNIFORMS)
    Dict{Symbol, GLenum}(ntuple(uniformLength) do i # take size and name
        name, typ = glGetActiveUniform(program, i-1)
    end)
end
function attribute_name_type(program::GLuint)
    uniformLength = glGetProgramiv(program, GL_ACTIVE_ATTRIBUTES)
    Dict{Symbol, GLenum}(ntuple(uniformLength) do i
    	name, typ = glGetActiveAttrib(program, i-1)
    end)
end
function istexturesampler(typ::GLenum)
    return (
        typ == GL_SAMPLER_BUFFER || typ == GL_INT_SAMPLER_BUFFER || typ == GL_UNSIGNED_INT_SAMPLER_BUFFER ||
    	typ == GL_IMAGE_2D ||
        typ == GL_SAMPLER_1D || typ == GL_SAMPLER_2D || typ == GL_SAMPLER_3D ||
        typ == GL_UNSIGNED_INT_SAMPLER_1D || typ == GL_UNSIGNED_INT_SAMPLER_2D || typ == GL_UNSIGNED_INT_SAMPLER_3D ||
        typ == GL_INT_SAMPLER_1D || typ == GL_INT_SAMPLER_2D || typ == GL_INT_SAMPLER_3D ||
        typ == GL_SAMPLER_1D_ARRAY || typ == GL_SAMPLER_2D_ARRAY ||
        typ == GL_UNSIGNED_INT_SAMPLER_1D_ARRAY || typ == GL_UNSIGNED_INT_SAMPLER_2D_ARRAY ||
        typ == GL_INT_SAMPLER_1D_ARRAY || typ == GL_INT_SAMPLER_2D_ARRAY
    )
end

gl_promote{T <: Integer}(x::Type{T})       = Cint
gl_promote(x::Type{Union{Int16, Int8}})    = x

gl_promote{T <: Unsigned}(x::Type{T})      = Cuint
gl_promote(x::Type{Union{UInt16, UInt8}})  = x

gl_promote{T <: AbstractFloat}(x::Type{T}) = Float32
gl_promote(x::Type{Float16})               = x

gl_promote{T <: UFixed}(x::Type{T})        = UFixed32
gl_promote(x::Type{UFixed16})              = x
gl_promote(x::Type{UFixed8})               = x

typealias Color3{T} Colorant{T, 3}
typealias Color4{T} Colorant{T, 4}

gl_promote(x::Type{Bool})                  = GLboolean
gl_promote{T <: Gray}(x::Type{T})          = Gray{gl_promote(eltype(T))}
gl_promote{T <: Color3}(x::Type{T})        = RGB{gl_promote(eltype(T))}
gl_promote{T <: Color4}(x::Type{T})        = RGBA{gl_promote(eltype(T))}
gl_promote{T <: BGRA}(x::Type{T})          = BGRA{gl_promote(eltype(T))}
gl_promote{T <: BGR}(x::Type{T})           = BGR{gl_promote(eltype(T))}

gl_promote{T <: FixedVector}(x::Type{T}) = similar_type(T, gl_promote(eltype(T)))

gl_promote{T <: HomogenousMesh}(x::Type{T}) = NativeMesh{T}

gl_convert{T <: Number}(x::T) = gl_promote(T)(x)
gl_convert{T <: Colorant}(x::T) = gl_promote(T)(x)
gl_convert{T <: AbstractMesh}(x::T) = gl_convert(convert(GLNormalMesh, x))
gl_convert{T <: HomogenousMesh}(x::T) = gl_promote(T)(x)
gl_convert{T <: HomogenousMesh}(x::Signal{T}) = gl_promote(T)(x)

gl_convert{T<:Colorant}(s::Vector{Matrix{T}}) = Texture(s)
gl_convert(s::AABB) = s
gl_convert(s::Void) = s

isa_gl_struct(x::Array) = false
isa_gl_struct(x::NATIVE_TYPES) = false
isa_gl_struct(x::Colorant) = false
function isa_gl_struct{T}(x::T)
    !isleaftype(T) && return false
    if T <: Tuple
        return false
    end
    fnames = fieldnames(T)
    !isempty(fnames) && all(name -> isleaftype(fieldtype(T, name)) && isbits(getfield(x, name)), fnames)
end
function gl_convert_struct{T}(x::T, uniform_name::Symbol)
    if isa_gl_struct(x)
        return Dict{Symbol, Any}(map(fieldnames(x)) do name
            (Symbol("$uniform_name.$name") => gl_convert(getfield(x, name)))
        end)
    else
        error("can't convert $x to a OpenGL type. Make sure all fields are of a concrete type and isbits(FieldType)-->true")
    end
end

# native types don't need convert!
gl_convert{T <: NATIVE_TYPES}(a::T) = a

gl_convert{T <: NATIVE_TYPES}(s::Signal{T}) = s
gl_convert{T}(s::Signal{T}) = const_lift(gl_convert, s)

for N=1:4
    @eval gl_convert{T}(x::FixedVector{$N, T}) = map(gl_promote(T), x)
end
for N=1:4, M=1:4
    @eval gl_convert{T}(x::Mat{$N, $M, T}) = map(gl_promote(T), x)
end

gl_convert{T <: Face}(a::Vector{T}) = indexbuffer(s)
gl_convert{T <: NATIVE_TYPES}(::Type{T}, a::NATIVE_TYPES; kw_args...) = a
function gl_convert{T <: GPUArray, X, N}(::Type{T}, a::Array{X, N}; kw_args...)
    T(map(gl_promote(X), a); kw_args...)
end
function gl_convert{T <: Texture, X}(::Type{T}, a::Vector{Array{X, 2}}; kw_args...)
    T(a; kw_args...)
end

function gl_convert{T <: GPUArray, X, N}(::Type{T}, a::Signal{Array{X, N}}; kw_args...)
    TGL = gl_promote(X)
    s = (X == TGL) ? a : const_lift(map, TGL, a)
    T(s; kw_args...)
end

gl_convert(f::Function, a) = f(a)


text:
module Covid_Forecast

# Write your package code here.
#using WebIO
using IJulia
#using Conda #https://juliagizmos.github.io/WebIO.jl/latest/providers/ijulia/
#using Pluto
#using Interact
#using CSV
#using DataFrames
#Conda.pip_interop(true)
#Conda.pip("install", "webio_jupyter_extension")
notebook()
#Pluto.run()

end

ext: jl


text:
#=

# 102: 1D Stationary convection-diffusion equation
([source code](SOURCE_URL))

Solve the equation

```math
-\nabla ( D \nabla u - v u) = 0
```
in $\Omega=(0,1)$ with boundary condition $u(0)=0$ and $u(1)=1$.
$v$ could be e.g. the velocity
of a moving medium or the gradient of an electric field.

This is a convection dominant second order boundary value problem which obeys
a local and a global maximum principle:
the solution which is bounded by the values at the boundary and has no local extrema in the
interior.
If $v$ is large compared to $D$, a boundary layer is observed.

The maximum principle of the solution can only be guaranteed it the discretization is
performed accordingly: the flux function must monotonically increase in the first argument
and monotonically decrease in the second argument.

The example describes three possible ways to define the flux function and demonstrates
the impact on the qualitative properties of the solution.

=#

module Example102_StationaryConvectionDiffusion1D
using Printf
using VoronoiFVM
using ExtendableGrids
using GridVisualize

## Central difference flux. The velocity term is discretized using the
## average of the solution in the endpoints of the grid. If the local Peclet
## number v*h/D>1, the monotonicity property is lost.  Grid refinement
## can fix this situation by decreasing $h$.

function central_flux!(f,u0,edge,data)
    u=unknowns(edge,u0)
    f_diff=data.D*(u[1,1]-u[1,2])
    vh=project(edge,data.v)
    f[1]=f_diff+vh*(u[1,1]+u[1,2])/2
end

## The simple upwind flux corrects the monotonicity properties essentially
## via brute force and loses one order of convergence for small $h$ compared
## to the central flux.

function upwind_flux!(f,u0,edge,data)
    u=unknowns(edge,u0)
    fdiff=data.D*(u[1,]-u[1,2])
    vh=project(edge,data.v)
    if vh>0
        f[1]=fdiff+vh*u[1,1]
    else
        f[1]=fdiff+vh*u[1,2]
    end
end

## The exponential fitting flux has the proper monotonicity properties and
## kind of interpolates in a clever way between central
## and upwind flux. It can be derived by solving the two-point boundary value problem
## at the grid interval analytically.

## Bernoulli function used in the exponential fitting discretization
function bernoulli(x)
    if abs(x)<nextfloat(eps(typeof(x)))
        return 1
    end
    return x/(exp(x)-1)
end

function exponential_flux!(f,u0,edge,data)
    u=unknowns(edge,u0)
    vh=project(edge,data.v)
    Bplus= data.D*bernoulli(vh/data.D)
    Bminus=data.D*bernoulli(-vh/data.D)
    f[1]=Bminus*u[1,1]-Bplus*u[1,2]
end

function calculate(grid,data,flux,verbose)

    sys=VoronoiFVM.System(grid,VoronoiFVM.Physics(flux=flux, data=data))

    ## Add species 1 to region 1
    enable_species!(sys,1,[1])

    ## Set boundary conditions
    boundary_dirichlet!(sys,1,1,0.0)
    boundary_dirichlet!(sys,1,2,1.0)

    ## Create a solution array
    inival=unknowns(sys,inival=0.5)
    solution=unknowns(sys)

    ## Create solver control info
    control=VoronoiFVM.NewtonControl()
    control.verbose=verbose

    ## Stationary solution of the problem
    solve!(solution,inival,sys, control=control)
    return solution
end

function main(;n=10,Plotter=nothing,verbose=false,D=0.01,v=1.0)

    ## Create a one-dimensional discretization
    h=1.0/convert(Float64,n)
    grid=VoronoiFVM.Grid(collect(0:h:1))

    data=(v=[v],D=D)

    # Calculate three stationary solutions with different ways to calculate flux
    solution_exponential=calculate(grid,data,exponential_flux!,verbose)
    solution_upwind=calculate(grid,data,upwind_flux!,verbose)
    solution_central=calculate(grid,data,central_flux!,verbose)

    # Visualize solutions using GridVisualize.jl
    p=GridVisualizer(Plotter=Plotter,layout=(3,1))
    scalarplot!(p[1,1],grid,solution_exponential[1,:],title="exponential")
    scalarplot!(p[2,1],grid,solution_upwind[1,:],title="upwind")
    scalarplot!(p[3,1],grid,solution_central[1,:],title="centered",show=true)

    ## Return test value
    return sum(solution_exponential)+sum(solution_upwind)+sum(solution_central)
end

function test()
    testval=2.523569744561089
    main() ≈ testval
end

end


text:
using Test
using CalculatedABC

using DelimitedFiles
using Random, Distributions
Random.seed!(123)

include("test_curve.jl")
include("test_analysis.jl")
include("test_gini.jl")


text:
struct StatelessNormal{T} <: AbstractStatelessPolicy where {T}
    μ::T
    σ::T

    function StatelessNormal(::Type{T}, num_actions::Int, std=1.0) where {T}
        μ = zeros(T, num_actions)
        σ = ones(T, num_actions) .* convert.(T, std)
        return new{typeof(μ)}(μ, σ)
    end

    function StatelessNormal(num_actions::Int, std=1.0)
        return StatelessNormal(Float64, num_actions, std)
    end

end

function params(π::StatelessNormal)
    return (π.μ, π.σ)
end

function paramsvec(π::StatelessNormal)
    return vcat(copy(vec(π.μ)), copy(vec(π.σ)))
end

function paramsfromvec!(π::StatelessNormal, θ)
    nμ = length(π.μ)
    @. π.μ = @view θ[1:nμ]
    @. π.σ = @view θ[nμ+1:end]
    return π
end

function rrule(::typeof(paramsfromvec!), π::StatelessNormal, θ)
    nμ = length(π.μ)
    @. π.μ = @view θ[1:nμ]
    @. π.σ = @view θ[nμ+1:end]
    function paramsfromvec_statelessnormal!(ȳ)
        dθ = vcat(vec(ȳ.μ), vec(ȳ.σ))
        return NoTangent(), ȳ, dθ
    end
    return π, paramsfromvec_statelessnormal!
end

struct LinearNormal{T,TS} <: AbstractPolicy  where {T,TS}
    W::T
    σ::TS

    function LinearNormal(::Type{T}, num_features::Int, num_actions::Int, std=1.0) where {T}
        W = zeros(T, (num_features, num_actions))
        σ = ones(T, num_actions) .* convert.(T, std)
        return new{typeof(W),typeof(σ)}(W, σ)
    end

    function LinearNormal(num_features::Int, num_actions::Int, std=1.0)
        return LinearNormal(Float64, num_features, num_actions, std)
    end
end

function params(π::LinearNormal)
    return (π.W, π.σ)
end

function paramsvec(π::LinearNormal)
    return vcat(copy(vec(π.W)), copy(vec(π.σ)))
end

function paramsfromvec!(π::LinearNormal, θ)
    nW = length(π.W)
    θW = @view θ[1:nW]
    π.W .= reshape(θW, size(π.W))
    @. π.σ = @view θ[nW+1:end]
    return π
end

function rrule(::typeof(paramsfromvec!), π::LinearNormal, θ)
    nW = length(π.μ)
    θW = @view θ[1:nW]
    π.W .= reshape(θW, size(π.W))
    @. π.σ = @view θ[nW+1:end]
    function paramsfromvec_normal!(ȳ)
        dθ = vcat(vec(ȳ.W), vec(ȳ.σ))
        return NoTangent(), ȳ, dθ
    end
    return π, paramsfromvec_normal!
end

struct NormalBuffer{TA, T, TS} <: Any
    action::TA
    μ::T
    ψ::TS

    function NormalBuffer(π::StatelessNormal)
        T = eltype(π.μ)
        n = length(π.μ)
        μ = zeros(T, n)
        ψ = (zero(π.μ), zero(π.σ))
        action = [zeros(T, n)]
        return new{typeof(action),typeof(μ),typeof(ψ)}(action, μ, ψ)
    end

    function NormalBuffer(π::LinearNormal)
        T = eltype(π.W)
        n = size(π.W,2)
        μ = zeros(T, n)
        ψ = (zero(π.W), zero(π.σ))
        action = [zeros(T, n)]
        return new{typeof(action),typeof(μ), typeof(ψ)}(action, μ, ψ)
    end
end

function logpdf_normal(z::T, μ::T, σ::T) where {T}
    tmp = -log(sqrt(T(2.0) * T(π)))
    logp = tmp - (z - μ)^2 / (T(2.0) * σ^2) - log(σ)
    return logp
end

function normal_grad_mu!(G, z, μ, σ)
    @. G = (z - μ) / σ^2
end

function normal_grad_sigma!(G, z, μ, σ)
    T = eltype(z)
    @. G = (T(-1) + ((z - μ) / σ)^2) / σ
end

function linear_normal_gradmu!(G,z,μ,σ,x)
    N = length(z)

    for i in 1:N
        @. G[i, :] = x * (z[i] - μ[i]) / σ[i]^2
    end
end

function (π::StatelessNormal)()
    μ, σ = π.μ, π.σ
    N = length(μ)
    if N > 1
        return MvNormal(μ, σ)
    else
        return Normal(μ[1], σ[1])
    end
end

function (π::StatelessNormal)(buff::NormalBuffer)
    μ, σ = π.μ, π.σ
    N = length(μ)
    if N > 1
        return MvNormal(μ, σ)
    else
        return Normal(μ[1], σ[1])
    end
end

function logpdf(π::StatelessNormal, a)
    μ, σ = π.μ, π.σ
    logp = sum(logpdf_normal.(a, μ, σ))
    return logp
end

function logpdf!(buff::NormalBuffer, π::StatelessNormal, a)
    μ, σ = π.μ, π.σ
    @. buff.μ = logpdf_normal(a, μ, σ)
    logp = sum(buff.μ)
    return logp
end

function grad_logpdf!(ψ, π::StatelessNormal, a)
    μ, σ = π.μ, π.σ
    ψμ, ψσ = ψ[1], ψ[2]
    normal_grad_mu!(ψμ, a, μ, σ)
    normal_grad_sigma!(ψσ, a, μ, σ)

    logp = sum(logpdf_normal.(a, μ, σ))

    return logp, buff.ψ
end

function grad_logpdf!(buff::NormalBuffer, π::StatelessNormal, a)
    μ, σ = π.μ, π.σ
    ψμ, ψσ = buff.ψ[1], buff.ψ[2]
    normal_grad_mu!(ψμ, a, μ, σ)
    normal_grad_sigma!(ψσ, a, μ, σ)
    @. buff.μ = logpdf_normal(a, μ, σ)
    logp = sum(buff.μ)
    return logp, buff.buff.ψ
end

function sample_with_trace!(ψ, action, π::StatelessNormal)
    μ, σ = π.μ, π.σ

    T = eltype(μ)
    @. action[1] = μ + σ * randn((T,))
    logp = grad_logpdf!(ψ, π, action)

    return action[1], logp, buff.ψ
end

function sample_with_trace!(buff::NormalBuffer, π::StatelessNormal)
    μ, σ = π.μ, π.σ

    T = eltype(μ)
    @. buff.action[1] = μ + σ * randn((T,))
    ψ = buff.ψ
    logp = grad_logpdf!(ψ, π, buff.action[1])

    return buff.action[1], logp, buff.ψ
end

function (π::LinearNormal)(s)
    W, σ = π.W, π.σ
    μ = W's
    N = length(σ)
    if N > 1
        return MvNormal(μ, σ)
    else
        return Normal(μ[1], σ[1])
    end
end

function (π::LinearNormal)(buff::NormalBuffer, s)
    W, σ = π.W, π.σ
    mul!(buff.μ, W', s)
    N = length(σ)
    if N > 1
        return MvNormal(buff.μ, σ)
    else
        return Normal(buff.μ[1], σ[1])
    end
end

function logpdf(π::LinearNormal, s, a)
    W, σ = π.W, π.σ
    μ = W's
    logp = sum(logpdf_normal.(a, μ, σ))
    return logp
end

function logpdf(π::LinearNormal, s::AbstractMatrix, a)
    W, σ = π.W, π.σ
    μ = W's
    logp = sum(logpdf_normal.(a, μ, σ), dims=1)
    return logp
end

function logpdf!(buff::NormalBuffer, π::LinearNormal, s, a)
    W, σ = π.W, π.σ
    mul!(buff.μ, W', s)
    T = eltype(σ)
    logp = T(0.0)
    for i in 1:length(buff.μ)
        logp += logpdf_normal(a[i], buff.μ[i], σ[i])
    end
    return logp
end

function grad_logpdf!(ψ, π::LinearNormal, s, a)
    W, σ = π.W, π.σ
    μ = W's

    ψw, ψσ = ψ[1]', ψ[2]

    linear_normal_gradmu!(ψw, a, μ, σ, s)
    normal_grad_sigma!(ψσ, a, μ, σ)

    logp = sum(logpdf_normal.(a, μ, σ))

    return logp, ψ
end

function grad_logpdf!(buff::NormalBuffer, π::LinearNormal, s, a)
    W, σ = π.W, π.σ
    mul!(buff.μ,W', s)
    ψ = buff.ψ
    ψw, ψσ = ψ[1]', ψ[2]

    linear_normal_gradmu!(ψw, a, buff.μ, σ, s)
    normal_grad_sigma!(ψσ, a, buff.μ, σ)

    T = eltype(σ)
    logp = T(0.0)
    for i in 1:length(buff.μ)
        logp += logpdf_normal(a[i], buff.μ[i], σ[i])
    end

    return logp, ψ
end

function sample_with_trace!(ψ, action, π::LinearNormal, s)
    W, σ = π.W, π.σ
    μ = W's

    T = eltype(σ)
    @. action[1] = μ + σ * randn((T,))

    ψw, ψσ = ψ[1]', ψ[2]

    linear_normal_gradmu!(ψw, action, μ, σ, s)
    normal_grad_sigma!(ψσ, action, μ, σ)

    logp = sum(logpdf_normal.(action, μ, σ))

    return action[1], logp, buff.ψ
end

function sample_with_trace!(buff::NormalBuffer, π::LinearNormal, s)
    μ, action, ψ = buff.μ, buff.action[1], buff.ψ
    W, σ = π.W, π.σ
    mul!(μ, W', s)

    T = eltype(σ)
    @. action = μ + σ * randn((T,))

    ψw, ψσ = ψ[1]', ψ[2]

    linear_normal_gradmu!(ψw, action, μ, σ, s)
    normal_grad_sigma!(ψσ, action, μ, σ)

    logp = T(0.0)
    for i in 1:length(μ)
        logp += logpdf_normal(action[i], μ[i], σ[i])
    end

    return action, logp, buff.ψ
end


text:
# Types for storing info about reduction tree
type ReductionTree
  Q
  leaves::Array{Any,1}
end

type TSQRQ
  Q::SharedArray
  T::RemoteRef
end

#
# Functions to map to leaf nodes and then reduce in binary tree
#
function tsqr_map( F, A::SharedArray )
  m,n = size(A)
  p = procs(A)
  refs = Array(RemoteRef, length(p))
  for (i, idx) in enumerate(Base.splitrange(m,length(p)))
    refs[i] = remotecall( p[i], F, A, idx )
  end
  return refs
end

function tsqr_reduce( F, A::SharedArray, refs )
  p = procs(A)
  np = length(p)
  for k = 1:ceil(Int64, log2(np))
    for i = 1:2^k:np
      if i + 2^(k-1) <= np
        j = i+2^(k-1)
        refs[i] = remotecall( p[i], F, A, refs[i], refs[j] )
      end
    end
  end
  return refs[1]
end

#
# TSQR, where A is overwritten by implicit Q
#
function tsqr!( A::SharedArray, I::UnitRange )
  m,n = size(A)
  # Factor this segment of A, overwriting with R and implicit Q
  F = qrfact!(sub(sdata(A), I, :))
  T = RemoteRef()
  put!( T, ReductionTree(F,[]) )
  # Return a COPY of R, since LAPACK's tpqrt appears to write over the lower
  # triangle
  return F[:R], T
end

function tsqr!( A::SharedArray, rr1::RemoteRef, rr2::RemoteRef )
  m,n = size(A)
  R1, T1 = fetch(rr1)
  R2, T2 = fetch(rr2)
  F = tpqrt!( R1, R2 )
  T = RemoteRef()
  put!( T, ReductionTree(F, [T1; T2]) )
  return R1, T
end

function tsqr!( A::SharedArray )
  rr = tsqr_map( tsqr!, A )
  R,T = fetch(tsqr_reduce( tsqr!, A, rr ))
  return TSQRQ(A,T), R
end

#
# TSQR, A is not overwritten and Q is discarded. Only R is returned
#
function tsqr( A::SharedArray, I::UnitRange )
  F = qrfact(sub(sdata(A),I,:))
  return F[:R]
end

function tsqr( A::SharedArray, rr1::RemoteRef, rr2::RemoteRef )
  m,n = size(A)
  R1, R2 = fetch(rr1), fetch(rr2)
  tpqrt!( R1, R2 )
  return R1
end

function tsqr( A::SharedArray )
  rr = tsqr_map( tsqr, A )
  R = fetch( tsqr_reduce( tsqr, A, rr ) )
  return R
end

#
# Construct full Q from explicit representation
#
function constructQ!( Q::TSQRQ )
  m,n = size(Q.Q)
  remotecall_wait( Q.T.where, constructQ!, Q.T, eye(n,n) )
  return Q.Q
end

function constructQ!( T::RemoteRef, A::Array{Float64,2} )
  T = fetch(T)

  # If at leaf nodes, use gemqrtr!
  if length(T.leaves) == 0
    m,n = size(T.Q.factors)
    R = zeros(m,n)
    for j = 1:n
      for i = 1:j
        @inbounds R[i,j] = A[i,j]
      end
    end
    # Multiply Q by the new array A
    gemqrtr!(T.Q,R)
    copy!(T.Q.factors,R)

  # Otherwise, use tpmqrt
  else
    # Multiply Q by the new array R
    Q1,Q2 = tpmqrt!( T.Q, A )
    # Perform on children
    @sync begin
      @async remotecall_wait( T.leaves[1].where, constructQ!, T.leaves[1], Q1)
      @async remotecall_wait( T.leaves[2].where, constructQ!, T.leaves[2], Q2)
    end
  end
end


text:
# We actually don't need the RigidBody simulation record.
struct ThrowAwaySink <: OdeIntegrators.OdeResultsSink
end
OdeIntegrators.initialize(::ThrowAwaySink, args...) = nothing
OdeIntegrators.process(::ThrowAwaySink, args...) = nothing


text:
# This file is a part of RRRMC.jl. License is MIT: http://github.com/carlobaldassi/RRRMC.jl/LICENCE.md

module RE

using ExtractMacro
using ..Interface
using ..Common

using ...RRRMC # this is silly but it's required for correct cross-linking in docstrings, apparently

export GraphRE, GraphRobustEnsemble, REenergies

import ..Interface: energy, delta_energy, neighbors, allΔE,
                    update_cache!, delta_energy_residual

import Base: iterate, length, eltype

function logcoshratio(a, b)
    #return log(cosh(a) / cosh(b))
    #log((exp(a) + exp(-a)) / (exp(b) + exp(-b)))
    a = abs(a)
    b = abs(b)
    return a - b + (log1p(exp(-2a)) - log1p(exp(-2b)))
end

fk(μ̄::Integer, γ::Real, β::Real) = logcoshratio(γ * (μ̄ + 1), γ * (μ̄ - 1)) / β

mutable struct GraphRE{M,γ,β} <: DiscrGraph{Float64}
    N::Int
    Nk::Int
    μ::IVec
    cache::LocalFields{Float64}

    function GraphRE{M,γ,β}(N::Integer, μ::Union{Nothing,IVec} = nothing) where {M,γ,β}
        isa(M, Int) || throw(ArgumentError("invalid parameter M, expected Int, given: $(typeof(M))"))
        M > 2 || throw(ArgumentError("M must be greater than 2, given: $M"))
        isa(β, Float64) || throw(ArgumentError("invalid parameter β, expected Float64, given: $(typeof(β))"))
        isa(γ, Float64) || throw(ArgumentError("invalid parameter γ, expected Float64, given: $(typeof(γ))"))
        N % M == 0 || throw(ArgumentError("N must be divisible by M, given: N=$N M=$M"))
        Nk = N ÷ M

        μ::IVec = (μ ≡ nothing ? zeros(Int, Nk) : μ)
        @assert length(μ) == Nk
        cache = LocalFields{Float64}(N)
        return new{M,γ,β}(N, Nk, μ, cache)
    end
end

@doc """
    GraphRE{M,γ,β}(N::Integer) <: DiscrGraph

An auxiliary `DiscrGraph` used to implement the interactions in the
Robust Ensemble.

It is only useful when implementing other graph types; see [`GraphRobustEnsemble`](@ref).
""" -> GraphRE{M,γ,β}(N::Integer)

GraphRE(X::GraphRE{M,oldγ,β}, newγ::Float64) where {M,oldγ,β} = GraphRE{M,newγ,β}(X.N, X.μ)

@generated function ΔElist(::Type{GraphRE{M,γ,β}}) where {M,γ,β}
    Expr(:tuple, ntuple(d->fk(2*(d - 1 - (M-1) >>> 0x1) - iseven(M), γ, β), M)...)
end
lstind(μ̄::Int, M::Int) = (μ̄ + M-1) >>> 0x1 + 1

function getk(X::GraphRE{M,γ,β}, μ̄::Int) where {M,γ,β}
    @inbounds k = ΔElist(GraphRE{M,γ,β})[lstind(μ̄, M)]
    return k
end

function energy(X::GraphRE{M,γ,β}, C::Config) where {M,γ,β}
    # @assert X.N == C.N
    @extract X : Nk μ cache
    @extract cache : lfields lfields_last
    @extract C : s

    fill!(μ, 0)
    j = 0
    for i = 1:Nk, k = 1:M
        j += 1
        # @assert j == k + (i-1) * M
        σj = 2s[j] - 1
        μ[i] += σj
    end

    n = 0.0
    for i = 1:Nk
        n -= log(2cosh(γ * μ[i])) / β
    end

    j = 0
    for i = 1:Nk, k = 1:M
        j += 1
        # @assert j == k + (i-1) * M
        σj = 2s[j] - 1
        μ̄ = μ[i] - σj
        # kj = fk(μ̄, γ, β)
        kj = getk(X, μ̄)
        # @assert kj == fk(μ̄, γ, β)
        lfields[j] = σj * kj
    end
    cache.move_last = 0
    fill!(lfields_last, 0.0)
    return n
end

function kinterval(move::Integer, M::Integer)
    j0 = move - ((move-1) % M)
    j1 = j0 + M - 1
    return j0:j1
end

function update_cache!(X::GraphRE{M,γ,β}, C::Config, move::Int) where {M,γ,β}
    # @assert X.N == C.N
    # @assert 1 ≤ move ≤ C.N
    @extract C : N s
    @extract X : Nk μ cache

    @inbounds σx = 2s[move] - 1
    i = (move-1) ÷ M + 1

    Ux = kinterval(move, M)

    @extract cache : lfields lfields_last move_last
    if move_last == move
        @inbounds begin
            #for y in neighbors(X, move)
            for y in Ux
                lfields[y], lfields_last[y] = lfields_last[y], lfields[y]
            end
            #lfields[move] = -lfields[move]
            #lfields_last[move] = -lfields_last[move]
        end
        μ[i] += 2σx
        return
    end

    @inbounds begin
        μnew = μ[i] + 2σx
        μ[i] = μnew
        #for y in neighbors(X, move)
            # @assert y ≠ move
        for y in Ux
            σy = 2s[y] - 1
            μ̄ = μnew - σy
            # @assert -M+1 ≤ μ̄ ≤ M-1    μ̄
            ky = getk(X, μ̄)
            # @assert ky == fk(μ̄, γ, β)
            lfields_last[y] = lfields[y]
            lfields[y] = σy * ky
        end
        #lfm = lfields[move]
        #lfields_last[move] = lfm
        #lfields[move] = -lfm
    end
    cache.move_last = move

    # lfields_bk = copy(lfields)
    # energy(X, C)
    # lfields_bk ≠ lfields && @show move hcat(lfields,lfields_bk)
    # @assert lfields_bk == lfields

    return
end

@inline function delta_energy(X::GraphRE, C::Config, move::Int)
    # @assert X.N == C.N
    # @assert 1 ≤ move ≤ C.N
    @extract X : cache
    @extract cache : lfields

    @inbounds Δ = lfields[move]
    return Δ
end

struct CavityRange
    j0::Int
    j1::Int
    jX::Int
    function CavityRange(j0::Integer, j1::Integer, jX::Integer)
        j0 ≤ jX ≤ j1 || throw(ArgumentError("invalid CavityRange parameters, expected j0≤jX≤j1, given: j0=$j0, j1=$j1, jX=$X"))
        return new(j0, j1, jX)
    end
end

function iterate(crange::CavityRange, j = crange.j0 + (crange.jX == crange.j0))
    j > crange.j1 && return nothing
    return j, j + 1 + (j == crange.jX - 1)
end

# start(crange::CavityRange) = crange.j0 + (crange.jX == crange.j0)
# done(crange::CavityRange, j) = j > crange.j1
# @inline function next(crange::CavityRange, j)
#     @extract crange : j0 j1 jX
#     # @assert j ≠ jX
#     nj = j + 1
#     nj += (nj == jX)
#     return (j, nj)
# end
length(crange::CavityRange) = crange.j1 - crange.j0
eltype(::Type{CavityRange}) = Int

@inline function neighbors(X::GraphRE{M}, j::Int) where {M}
    j0 = j - ((j-1) % M)
    j1 = j0 + M - 1
    return CavityRange(j0, j1, j)
end

@generated function allΔE(::Type{GraphRE{M,γ,β}}) where {M,γ,β}
    K = M - 1
    iseven(K) ? Expr(:tuple, ntuple(d->fk(2*(d-1), γ, β), K÷2+1)...) :
                Expr(:tuple, ntuple(d->fk(2d-1, γ, β), (K+1)÷2)...)
end

# Replicate an existsing graph

mutable struct GraphRobustEnsemble{M,γ,β,G<:AbstractGraph} <: DoubleGraph{DiscrGraph{Float64},Float64}
    N::Int
    Nk::Int
    X0::GraphRE{M,γ,β}
    X1::Vector{G}
    C1::Vector{Config}
    function GraphRobustEnsemble{M,γ,β,G}(N::Integer, g0::G, Gconstr, args...) where {M,γ,β,G}
        X0 = GraphRE{M,γ,β}(N)
        Nk = X0.Nk
        X1 = Array{G}(undef, M)
        X1[1] = g0
        for k = 2:M
            X1[k] = Gconstr(args...)
        end
        C1 = [Config(Nk, init=false) for k = 1:M]
        return new{M,γ,β,G}(N, Nk, X0, X1, C1)
    end
end

"""
    GraphRobustEnsemble(N::Integer, M::Integer, γ::Float64, β::Float64, Gconstr, args...) <: DoubleGraph

A `DoubleGraph` which implements a "Robust Ensemble" model, given any other Ising model previously defined.
This kind of graph can be simulated efficiently with [`rrrMC`](@ref).

`N` is the number of spins, `M` the number of replicas, `γ` the interaction strength,
`β` the inverse temperature. `GConstr` is the (original) graph constructor, and `args` the
arguments to the contructor.

This is similar to [`GraphLocalEntropy`](@ref RRRMC.GraphLocalEntropy), but the reference configuration is traced out.

See also [`REenergies`](@ref).
"""
function GraphRobustEnsemble(Nk::Integer, M::Integer, γ::Float64, β::Float64, Gconstr, args...)
    g0 = Gconstr(args...)
    G = typeof(g0)
    return GraphRobustEnsemble{M,γ,β,G}(Nk * M, g0, Gconstr, args...)
end

function update_cache!(X::GraphRobustEnsemble{M}, C::Config, move::Int) where {M}
    @extract X : X0 X1 C1
    k = mod1(move, M)
    i = (move - 1) ÷ M + 1

    spinflip!(X1[k], C1[k], i)

    update_cache!(X0, C, move)
end

function energy(X::GraphRobustEnsemble{M}, C::Config) where {M}
    # @assert X.N == C.N
    @extract X : Nk X0 X1 C1
    @extract C : s

    E = energy(X0, C)

    for k = 1:M
        s1 = C1[k].s
        for (i,j) = enumerate(k:M:(k + M * (Nk-1)))
            s1[i] = s[j]
        end
        E += energy(X1[k], C1[k])
    end

    return E
end

"""
    REenergies(X::GraphRobustEnsemble)

Returns a Vector with the individual energy (as defined by the original model)
of each replica in a [`GraphRobustEnsemble`](@ref) graph.
"""
function REenergies(X::GraphRobustEnsemble{M}) where {M}
    @extract X : X1 C1

    Es = zeros(M)

    for k = 1:M
        Es[k] = energy(X1[k], C1[k])
    end

    return Es
end

function delta_energy_residual(X::GraphRobustEnsemble{M}, C::Config, move::Int) where {M}
    @extract X : X1 C1

    k = mod1(move, M)
    i = (move - 1) ÷ M + 1

    return delta_energy(X1[k], C1[k], i)
end

function delta_energy(X::GraphRobustEnsemble, C::Config, move::Int)
    return delta_energy(X.X0, C, move) +
           delta_energy_residual(X, C, move)
end

# This is inefficient. On the other hand, it is
# basically only written for testing purposes...
function neighbors(X::GraphRobustEnsemble{M}, i::Int) where {M}
    @extract X : X0 X1
    jts = neighbors(X0, i)

    k = mod1(i, M)
    i1 = (i - 1) ÷ M + 1

    is = neighbors(X1[k], i1)

    return tuple(jts..., map(i->((i-1)*M+k), is)...)
end

end # module


text:
export impute_values

function logistic(x)
    return 1.0 ./ (1.0 .+ exp.(-x))
end

function compute_Z(model::MatFacModel;
                   covariates::Union{Nothing,Matrix}=nothing,
                   factor_set::Union{Nothing,Vector{Int64}}=nothing)

    if factor_set == nothing
        factor_set = collect(1:size(model.X,1))
    end
    X = model.X[factor_set,:]
    Y = model.Y[factor_set,:]

    Z = transpose(X)*Y
    Z .+= model.instance_offset
    Z .+= transpose(model.feature_offset)

    if covariates != nothing
        Z .+= covariates*model.instance_covariate_coeff
    end

    return Z
end

function apply_link_function(Z::Matrix, model::MatFacModel)

    type_vec = typeof.(model.losses)
    logistic_idx = (type_vec .== LogisticLoss)
    poisson_idx = (type_vec .== PoissonLoss)

    Z[:,logistic_idx] .= logistic(Z[:,logistic_idx])
    Z[:,poisson_idx] .= exp.(Z[:,poisson_idx])

    return Z
end

function impute_values(model::MatFacModel;
                       covariates::Union{Nothing,Matrix}=nothing,
                       factor_set::Union{Nothing,Vector{Int64}}=nothing)

    Z = compute_Z(model; covariates=covariates, factor_set=factor_set)
    Z = apply_link_function(Z, model)
    return Z
end


text:
using MPI
using StaticArrays
using LinearAlgebra
using ClimateMachine
using ClimateMachine.VariableTemplates
using ClimateMachine.Mesh.Topologies
using ClimateMachine.Mesh.Grids
using ClimateMachine.MPIStateArrays
using ClimateMachine.DGMethods
using Printf

using ClimateMachine.BalanceLaws

import ClimateMachine.BalanceLaws: vars_state, nodal_init_state_auxiliary!

using ClimateMachine.Mesh.Geometry: LocalGeometry

if !@isdefined integration_testing
    const integration_testing = parse(
        Bool,
        lowercase(get(ENV, "JULIA_CLIMA_INTEGRATION_TESTING", "false")),
    )
end

struct GradSphereTestModel{dir} <: BalanceLaw end

vars_state(m::GradSphereTestModel, ::Auxiliary, T) = @vars begin
    a::T
    ∇a::SVector{3, T}
end
vars_state(::GradSphereTestModel, ::Prognostic, T) = @vars()

function nodal_init_state_auxiliary!(
    ::GradSphereTestModel{dir},
    aux::Vars,
    tmp::Vars,
    g::LocalGeometry,
) where {dir}
    x, y, z = g.coord
    r = hypot(x, y, z)
    aux.a = r
    if !(dir isa HorizontalDirection)
        aux.∇a = g.coord / r
    else
        aux.∇a = SVector(0, 0, 0)
    end
end

using Test
function test_run(mpicomm, Ne_horz, Ne_vert, N, FT, ArrayType, direction)

    Rrange = range(FT(1 // 2); length = Ne_vert + 1, stop = FT(1))
    topl = StackedCubedSphereTopology(mpicomm, Ne_horz, Rrange)

    grid = DiscontinuousSpectralElementGrid(
        topl,
        FloatType = FT,
        DeviceArray = ArrayType,
        polynomialorder = N,
        meshwarp = Topologies.cubedshellwarp,
    )

    model = GradSphereTestModel{direction}()
    dg = DGModel(
        model,
        grid,
        nothing,
        nothing,
        nothing;
        state_gradient_flux = nothing,
    )

    exact_aux = copy(dg.state_auxiliary)

    continuous_field_gradient!(
        model,
        dg.state_auxiliary,
        (:∇a,),
        dg.state_auxiliary,
        (:a,),
        grid,
        direction,
    )

    err = euclidean_distance(exact_aux, dg.state_auxiliary)
    return err
end

let
    ClimateMachine.init()
    ArrayType = ClimateMachine.array_type()

    mpicomm = MPI.COMM_WORLD

    polynomialorder = 4
    base_Nhorz = 4
    base_Nvert = 2

    expected_result = [
        2.0924087890777517e-04
        1.3897932154337201e-05
        8.8256018429045312e-07
        5.5381072850485303e-08
    ]

    lvls = integration_testing || ClimateMachine.Settings.integration_testing ?
        length(expected_result) : 1

    @testset for FT in (Float64,)
        @testset for direction in (
            EveryDirection(),
            HorizontalDirection(),
            VerticalDirection(),
        )
            err = zeros(FT, lvls)
            @testset for l in 1:lvls
                Ne_horz = 2^(l - 1) * base_Nhorz
                Ne_vert = 2^(l - 1) * base_Nvert

                err[l] = test_run(
                    mpicomm,
                    Ne_horz,
                    Ne_vert,
                    polynomialorder,
                    FT,
                    ArrayType,
                    direction,
                )
                if !(direction isa HorizontalDirection)
                    @test err[l] ≈ expected_result[l]
                else
                    @test abs(err[l]) < FT(1e-13)
                end
            end
            @info begin
                msg = ""
                for l in 1:(lvls - 1)
                    rate = log2(err[l]) - log2(err[l + 1])
                    msg *= @sprintf("\n  rate for level %d = %e\n", l, rate)
                end
                msg
            end
        end
    end
end

nothing

max_stars_repo_name: mwarusz/CLIMA

max_issues_repo_name: mwarusz/CLIMA

max_forks_repo_name: mwarusz/CLIMA


text:
using LiftedMaps
using BlockArrays
using Test

import LiftedMaps.LinearMaps

@testset "LiftedMaps.jl" begin
    A = rand(10,10)
    B = rand(3,3)

    U = Base.OneTo(13)
    V = Base.OneTo(13)

    I = 1:10
    J = 11:13

    Alifted = LiftedMap(A, 1:10, 1:10, U, V)
    Blifted = LiftedMap(B, 11:13, 11:13, U, V)

    @test size(Alifted) == (13,13)
    @test size(Blifted) == (13,13)

    C = Alifted + Blifted

    @test size(C) == (13,13)
    @test C isa LinearMaps.LinearMap
    @test C isa LinearMaps.LinearCombination
    @test isblockdiagonal(C)

    D = Alifted + 2 * Blifted
    @test D isa LinearMaps.LinearMap
    @test D isa LinearMaps.LinearCombination
    @test isblockdiagonal(D)

    x = rand(13)
    x1 = x[I]
    x2 = x[J]

    y1 = A * x1
    y2 = 2 * B * x2
    y = [y1; y2]

    z = D * x
    @test y ≈ z

    Amatrix = Matrix(Alifted)
    @test Amatrix isa Matrix
    @test Amatrix[I,I] == A

    ax1 = blockedrange([10,3])
    ax2 = blockedrange([10,3])
    Alifted = @inferred LiftedMap(A,Block(1),Block(1),ax1,ax2)
    Blifted = @inferred LiftedMap(B,Block(2),Block(2),ax1,ax2)
    @test axes(Alifted) === (ax1,ax2)
    @test axes(Blifted) === (ax1,ax2)

    D = Alifted + 2 * Blifted
    z = @inferred D * x
end


text:
using Winston

x = linspace(0., 2pi, 30)
y = sin(x)

p = FramedArray(2, 2,
        title="title",
        aspect_ratio=0.75,
        xlabel="x label",
        ylabel="y label",
        uniform_limits=true,
        cellspacing=1.)

add(p, LineY(0, kind="dot"))

add(p[1,1], Curve(x, .25*y))
add(p[1,2], Curve(x, .50*y))
add(p[2,1], Curve(x, .75*y))
add(p[2,2], Curve(x, y))

file(p, "example5.png")


text:
if Sys.WORD_SIZE == 32
    @inline function first_cache_buffer(::Type{T}) where {T}
        reinterpret(Ptr{T}, ACACHEPTR[] + Threads.threadid() * first_cache_size() - first_cache_size())
    end
else
    @inline function first_cache_buffer(::Type{T}) where {T}
        MemoryBuffer{T}(undef, first_cache_size(T))
    end
end

BCache(i::Integer) = BCache(BCACHEPTR[]+cld_fast(second_cache_size()*i, Threads.nthreads()), i % UInt)
BCache(::Nothing) = BCache(BCACHEPTR[], nothing)

@inline Base.pointer(b::BCache) = b.p
@inline Base.unsafe_convert(::Type{Ptr{T}}, b::BCache) where {T} = Base.unsafe_convert(Ptr{T}, b.p)

function _use_bcache()
    while Threads.atomic_cas!(BCACHE_LOCK, zero(UInt), typemax(UInt)) != zero(UInt)
        pause()
    end
    return BCache(nothing)
end
@inline _free_bcache!(b::BCache{Nothing}) = reseet_bcache_lock!()

_use_bcache(::Nothing) = _use_bcache()
function _use_bcache(i)
    f = one(UInt) << i
    while (Threads.atomic_or!(BCACHE_LOCK, f) & f) != zero(UInt)
        pause()
    end
    BCache(i)
end
_free_bcache!(b::BCache{UInt}) = (Threads.atomic_xor!(BCACHE_LOCK, one(UInt) << b.i); nothing)

"""
  reset_bcache_lock!()

Currently not using try/finally in matmul routine, despite locking.
So if it errors for some reason, you may need to manually call `reset_bcache_lock!()`.
"""
@inline reseet_bcache_lock!() = (BCACHE_LOCK[] = zero(UInt); nothing)


text:
"""
The QNetwork type is a mutable structure that contains QNodes, QChannels, and
a dictionary of costs to weighted LightGraphs.
"""
mutable struct QNetwork <: QObject
    name::String
    nodes::Array{QNode}
    channels::Array{QChannel}
    # Dictionary between costs and corresponding LightGraphs weighted graphs.
    graph::Dict
    time::Float64

    # QNetwork() = new("QuNet", [], [], Dict())
    function QNetwork()
        network = new("QuNet", [], [], Dict(), 0.0)
        for cost_key in keys(zero_costvector())
            network.graph[cost_key] = SimpleWeightedDiGraph()
        end
        return network
    end
end

function QNetwork(graph::AbstractGraph; edge_costs=unit_costvector())
    network = QNetwork()

    vertexCount = 0
    for vertex in vertices(graph)
        vertexCount += 1
        add(network, BasicNode(string(vertexCount)))
    end

    edgeCount = 0
    for edge in edges(graph)
        edgeCount += 1
        add(network, BasicChannel(string(edgeCount), network.nodes[edge.src], network.nodes[edge.dst],
        edge_costs))
    end

    return network
end

function print(network::QNetwork)
    println("name: ", network.name)
    println("nodes: ", length(network.nodes))
    println("channels: ", length(network.channels))
end

# NOTE Outdated refresh_graph! below.
# We'll no longer use TemporalGraph structure to update graphs since it has
# A lot more bells and whistles now. Good practice to have something custom built.
# function refresh_graph!(network::QNetwork)
#    network.graph = TemporalGraph(network, 1).graph
# end

"""
    refresh_graph!(network::QNetwork)

Converts a QNetwork into several weighted LightGraphs (one
graph for each associated cost), then updates the QNetwork.graph attribute
with these new graphs.
"""
function refresh_graph!(network::QNetwork)

    refreshed_graphs = Dict{String, SimpleWeightedDiGraph}()

    for cost_key in keys(zero_costvector())
        refreshed_graphs[cost_key] = SimpleWeightedDiGraph()

        # Vertices
        add_vertices!(refreshed_graphs[cost_key], length(network.nodes))

        # Channels
        for channel in network.channels
            if channel.active == true
                src = findfirst(x -> x == channel.src, network.nodes)
                dest = findfirst(x -> x == channel.dest, network.nodes)
                weight = channel.costs[cost_key]
                add_edge!(refreshed_graphs[cost_key], src, dest, weight)
                add_edge!(refreshed_graphs[cost_key], dest, src, weight)
            end
        end
    end
    network.graph = refreshed_graphs
end

"""
    GridNetwork(dim::Int64, dimY::Int64)

Generates an X by Y grid network.
"""
function GridNetwork(dimX::Int64, dimY::Int64; edge_costs::Dict = unit_costvector())
    graph = LightGraphs.grid([dimX,dimY])
    net = QNetwork(graph; edge_costs=edge_costs)

    for x in 1:dimX
        for y in 1:dimY
            this = x + (y-1)*dimX
            net.nodes[this].location = Coords(x,y)
        end
    end

    refresh_graph!(net)
    return net
end

"""
```function update(network::QNetwork, new_time::Float64)```

The `update` function iterates through all objects in the network and updates
them according to a new global time.
"""
function update(network::QNetwork, new_time::Float64)
    old_time = network.time
    for node in network.nodes
        update(node, old_time, new_time)
    end

    for channel in network.channels
        update(channel, old_time, new_time)
    end
    network.time = new_time
end

"""
```function update(network::QNetwork)```

This instance of update iterates through all objects in the network and updates
them by the global time increment TIME_STEP defined in QuNet.jl
"""
function update(network::QNetwork)
    old_time = network.time
    new_time = old_time + TIME_STEP
    for node in network.nodes
        update(node, old_time, new_time)
    end

    for channel in network.channels
        update(channel, old_time, new_time)
    end
    network.time = new_time
end

"""
    getnode(network::QNetwork, id::Int64)

Fetch the node object corresponding to the given ID / Name
"""
function getnode(network::QNetwork, id::Int64)
    return network.nodes[id]
end

function getnode(network::QNetwork, name::String)
    for node in network.nodes
        if node.name == name
            return node
        end
    end
end

function getchannel(network::QNetwork, src::Union{Int64, String},
    dst::Union{Int64, String})
    src = getnode(network, src)
    dst = getnode(network, dst)
    for channel in network.channels
        if channel.src == src && channel.dest == dst
            return channel
        end
    end
end

"""
function add_qnode!(network::QNetwork; nodename::String="", nodetype::DataType=BasicNode)
    @assert nodetype in subtypes(QNode)
    new_node = nodetype(nodename)
    add(network, new_node)
end
"""

"""
function add_channel!(network::QNetwork, src::Union{Int64, String},
    dst::Union{Int64, String};
    name::string="", type=BasicChannel)
    src = getnode(src)
    dst = getnode(dst)
    @assert type in [BasicChannel, AirChannel]
    new_channel = type(name, src, dst)
    add(network, new_channel)
end
"""

max_stars_repo_name: FalafelGood/QuNet

max_issues_repo_name: FalafelGood/QuNet

max_forks_repo_name: FalafelGood/QuNet


text:
##################################
# Load Packages and Activate Env #
##################################
using Pkg
Pkg.activate(".")
Pkg.instantiate()

using MPIReco, JLD, PyPlot, LazyArtifacts

###############
# download data
###############
datadir = artifact"MDFStore"
store = MDFDatasetStore(datadir)

#####################
# run example scripts
#####################
include("smSampling.jl")
include("smRecovery.jl")


text:
# Constant folding

# Quick example
function constant_folding_example()
    a = 2 * 3
    b = a + 1
    return b > 1 ? 10 : 20
end
constant_folding_example()
@code_typed constant_folding_example()
#=
julia> @code_typed constant_folding_example()
CodeInfo(
1 ─     return 10
) => Int64
=#

# ------------------------------------------------------------
# Optional content below
# ------------------------------------------------------------

# Constant propagation
const constant_a = 1
function constant_propagation_example()
    a = constant_a
    b = a + 1
    return b > 1 ? 10 : 20
end
constant_propagation_example()
@code_typed constant_propagation_example()
#=
julia> @code_typed constant_propagation_example()
CodeInfo(
1 ─     return 10
) => Int64
=#

# If it's a variable, then all bets are off...
variable_a = 1
function propagate_my_variable()
    a = variable_a
    b = a + 1
    return b > 1 ? 10 : 20
end
propagate_my_variable()
@code_typed propagate_my_variable()
#=
julia> @code_typed propagate_my_variable()
CodeInfo(
1 ─ %1 = Main.variable_a::Any
│   %2 = (%1 + 1)::Any
│   %3 = (%2 > 1)::Any
└──      goto #3 if not %3
2 ─      return 10
3 ─      return 20
) => Int64
=#

# how bad is it?
@btime constant_propagation_example();
@btime propagate_my_variable();
#=
julia> @btime constant_propagation_example();
  0.030 ns (0 allocations: 0 bytes)

julia> @btime propagate_my_variable();
  40.749 ns (0 allocations: 0 bytes)
=#


text:
using HyrodiumHome
using Test

@testset "HyrodiumHome.jl" begin
    # Write your tests here.
end


text:
#sigtools.jl - functions that help with working with the sigmoid numbers.

@generated function clz{N, ES, mode}(x::Sigmoid{N, ES, mode})
  if (mode == :ubit)
    :(min(leading_zeros(@u(x) & (~@signbit)), N - 1))
  else
    :(min(leading_zeros(@u(x) & (~@signbit)), N))
  end
end

@generated function clo{N, ES, mode}(x::Sigmoid{N, ES, mode})
  if (mode == :ubit)
    :(min(leading_ones(@u(x) | (@signbit)), N - 1))
  else
    :(min(leading_ones(@u(x) | (@signbit)), N))
  end
end

@generated function regimebits{N, ES, mode}(x::Sigmoid{N, ES, mode})
  if mode == :ubit
    quote
      bitcount = ((@u(x) & (@invertbit) != 0) ? clo(x) : clz(x))
      min(bitcount, N - 2)
    end
  else
    quote
      bitcount = ((@u(x) & (@invertbit) != 0) ? clo(x) : clz(x))
      min(bitcount, N - 1)
    end
  end
end

max_stars_repo_name: ityonemo/SigmoidNumbers

max_issues_repo_name: ityonemo/SigmoidNumbers

max_forks_repo_name: ityonemo/SigmoidNumbers


text:
module Schema # Octo

tables = Dict{Core.TypeName, Dict{Symbol,Union{String,Symbol,Vector{Symbol}}}}()
validation_models = Dict{Core.TypeName, Function}()

struct Model
    typename::Core.TypeName
    info::Dict{Symbol,Union{String,Symbol,Vector{Symbol}}}
end

function _model(M::Type; table_name::String, primary_key::Union{Nothing,Symbol,Vector{Symbol}})::Model
    Tname = Base.typename(M)
    dict = Dict{Symbol,Union{String,Symbol,Vector{Symbol}}}(
        :table_name => table_name,
    )
    if primary_key !== nothing
        dict[:primary_key] = primary_key
    end
    tables[Tname] = dict
    Model(Tname, dict)
end

"""
    model(M::Type; table_name::String, primary_key::Union{Nothing,Symbol,String,Vector,Tuple}="id")::Model
"""
function model(M::Type; table_name::String, primary_key::Union{Nothing,Symbol,String,Vector,Tuple}="id")::Model
    if primary_key isa Vector{String}
        pk = Symbol.(primary_key)
    elseif primary_key isa String
        pk = Symbol(primary_key)
    elseif primary_key isa Tuple # legacy
        pk = collect(Symbol.(primary_key))
    else
        pk = primary_key
    end
    _model(M; table_name=table_name, primary_key=pk)
end

function Base.show(io::IO, mime::MIME"text/plain", model::Model)
    Base.show(io, mime, (; model.info...))
end

struct TableNameError <: Exception
    msg::String
end

struct PrimaryKeyError <: Exception
    msg::String
end

"""
    changeset(validations, M::Type)
"""
function changeset(validations, M::Type)
    Tname = Base.typename(M)
    validation_models[Tname] = validations
    Pair(Tname, validations)
end

include(joinpath("Schema", "validations.jl"))

end # module Octo.Schema


text:
import Base:haslength
import Serialization

export Workout, saveWorkout, loadWorkout, predict, fit!, hasmetric,
        freeze!, unfreeze!, validate

# Callback niceties from Flux.jl
call(f, xs...) = f(xs...)
runall(f) = f
runall(fs::AbstractVector) = () -> foreach(call, fs)

"""
The Workout keeps track of the progress of the training session. At least a model
and a loss function needs to be provided. Optional an optimizer and one or more
metrics can be specified.

If no optimizer is provided, SGD will be used. If no metrics are provided, only
the loss during training and validation will be registered (:loss and :val_loss).

# Usage

```julia
workout = Workout(model, L1Loss())

workout = Workout(model, CrossEntropy(), SGD())

workout = Workout(model, HingeLoss(), SGD(); acc=BinaryAccuracy())
```
"""
mutable struct Workout
    model::Layer
    loss::Union{Loss,Function}
    opt
    metrics::Vector{Pair}
    history::IdDict{Symbol,MetricStore}
    steps::Int
    epochs::Int

    function Workout(model::Layer, loss::Union{Loss,Function}, opt=Knet.SGD(); metrics...)
        new(model, loss, opt, collect(metrics), IdDict(), 0, 0)
    end
end

"""
Enable saving and loading of models by specialized KnetArray methods for Julia serialization
This will effectively move a GPU weight to the CPU before serialing it and move it back to
the GPU when deserializing.
"""
function Serialization.serialize(s::Serialization.AbstractSerializer, p::Knet.KnetArray)
    Serialization.serialize_type(s, typeof(p))
    Serialization.serialize(s, Array(p))
end

function Serialization.deserialize(s::Serialization.AbstractSerializer, t::Type{<:Knet.KnetArray})
    arr = Serialization.deserialize(s)
    return Knet.KnetArray(arr)
end

"""
Save a workout to a file. This will save all the state that is captured in the workout
and enables to continue at a later stage.
"""
function saveWorkout(workout::Workout, filename="workout_$(workout.steps).dat")::String
    # serialize
    Serialization.serialize(filename, workout)
    return filename
end

"""
Load a workout from file and return it.

# Usage

```julia
workout = loadWorkout("workout_1000.dat")
fit!(workout, mydata)
```
"""
function loadWorkout(filename)::Workout
    workout = Serialization.deserialize(filename)
    return workout
end

"""
Stop a training session. If this is called outside the scope of
a trianing session, just an error is thrown.
"""
function stop(workout::Workout, reason::String)
    throw(ErrorException(reason))
end

"""
Does the workout have any recorded values for a certain metric
"""
hasmetric(workout::Workout, metricname::Symbol)::Bool = haskey(workout.history, metricname)

"""
Function to generate the fully qualified metric name. It uses the metric name
and phase (train or valid) to come up with a unique name.
"""
function getmetricname(metric::Symbol, phase=:train)::Symbol
    metricname = phase == :train ? metric : Symbol("val_", metric)
end

"""
Invoke the configured metrics. The loss metric will always be logged and available.
Metrics are stored in the history attribute of the workout.
"""
function updatemetrics!(workout::Workout, loss::Number, y, y_pred, phase=:train)

    # First register the loss
    metricname = getmetricname(:loss, phase)
    e = get!(workout.history, metricname, SmartReducer())
    update!(e, workout.steps, loss)

    # And now run and register any additional metrics
    for (name,fn) in workout.metrics
        try
            metricname = getmetricname(name, phase)
            metricvalue = fn(y_pred, y)
            e = get!(workout.history, metricname, SmartReducer())
            update!(e, workout.steps, metricvalue)
        catch
            @warn "Failed executing metric." metricname maxlog=1
        end
    end
    return loss
end

"""
Get the metric value for a fully qualified metric name and a certain step. If
step is not provided the last step will be used. If no value is found the passed
function will not be invoked.

# Usage

```julia
getmetricvalue(workout, :val_loss) do value
    println("validation loss", value)
end
```
"""
function getmetricvalue(f::Function, workout::Workout, metricname::Symbol, step=workout.steps)
    if haskey(workout.history, metricname)
        m = workout.history[metricname]
        value = get(m.state, step, nothing)
        value !== nothing && f(value)
    end
end

"""
Utility function to calculate the gradients. Useful when checking for vanishing or
exploding gradients. The returned value is a Vector of (Param, Gradient).
"""
function gradients(workout::Workout, minibatch=first(workout.dl); convertor=autoConvertor)
    x, y = convertor(minibatch)
    J = Knet.@diff begin
        y_pred = workout.model(x)
        workout.loss(y_pred, y)
    end
    gradients = []
    for p in Knet.params(J)
        if p.opt === nothing; p.opt = Knet.clone(opt); end
        push!(gradients, (p, Knet.grad(J,p)))
    end
    return gradients
end

"""
Freeze a parameter so it no longer will be updated during training.
"""
function freeze!(p::Knet.Param)
    p.opt = :NOUPDATE
end

"""
Unfreeze a parameter so it will be updated again during training.
"""
function unfreeze!(p::Knet.Param)
    p.opt = nothing
end

"""
Perform the back propagation and update of weights in a single go.
"""
function back!(J::Knet.Tape, opt)
    for p in Knet.params(J)
        if p.opt === nothing; p.opt = Knet.clone(opt); end
        p.opt === :NOUPDATE && continue
        Knet.update!(p, Knet.grad(J,p))
    end
end

"""
Take a single step in updating the weights of a model. This function
will be invoked from fit! to do the actual learning.

For a minibatch (x,y) of data, the folowing sequence will be executed:

1. perform the forward pass
2. calculate the loss
3. update and remember the metrics, if any
4. do the backpropagation and update the weights
"""
function step!(workout::Workout, x, y; zerograd=true)
    workout.steps += 1
    J = Knet.@diff begin
        y_pred = workout.model(x)
        loss = workout.loss(y_pred, y)
        updatemetrics!(workout, Knet.value(loss), y, y_pred)
        loss
    end
    back!(J, workout.opt)
end

"""
Predict a sample, either a single value or a batch. Compared to invoking
the model directory with model(x), predit takes care of:

- Moving the data to the GPU if required.
- Making the data into a batch (controlled by makebatch parameter)

# Usage

```julia
x = randn(Float32, 224, 224, 3)
predict(model, x)
```
"""
function predict(model, x; makebatch=true, convertor=autoConvertor)
    x = makebatch ? addlast(x) : x
    y = model(convertor(x))
    makebatch ? droplast(y) : y
end

"""
Validate a minibatch and calculate the loss and metrics. Typically this function
is called from the fit! method. But if required can also be invoked directly.
"""
function validate(workout::Workout, x, y)
    y_pred = workout.model(x)
    loss = workout.loss(y_pred, y)
    updatemetrics!(workout, loss, y, y_pred, :valid)
end

"""
Train the model based on a supervised dataset and for a number
of epochs. fit! can be called multiple times and will continue
to train where is left of last time.

By default the fit! function will try to ensure the data is of the right
type (e.g. Float32) and on the right device (e.g. GPU) before feeding it to
the model.

# Usage

```julia
fit!(workout, traindata)
fit!(workout, traindata, testdata, epochs=50)
```
If you don't want any data conversion, just pass the identity funciton
as the convertor:

```julia
fit!(workout, traindata, convertor=identity)
```

"""
function fit!(workout::Workout, data, validation=nothing;
    epochs=1, convertor=autoConvertor, cb = ConsoleMeter())

    cb = runall(cb)

    for epoch in 1:epochs
        workout.epochs += 1
        d = data isa Function ? data() : data
        for minibatch in d
            step!(workout, convertor(minibatch)...)
            cb(workout, :train)
        end

        if validation !== nothing
            d = validation isa Function ? validation() : validation
            for minibatch in d
                validate(workout, convertor(minibatch)...)
            end
        end
        cb(workout, :valid)
    end
end

@debug "Loaded Training module"


text:
xform(d::EqualMix, _data) = xform(d.components[1], _data)

# StudentT(ν, μ = 0.0, σ = 1.0) = LocationScale(μ, σ, TDist(ν))

xform(d::Dists.Dirichlet, _data) = UnitSimplex(length(d.alpha))


text:
using CairoMakie
using GLMakie
using Colors

# Theme used
const BG = colorant"#191919"
const AX = colorant"#ababac"
const TX = colorant"#cdcdcf"
const GR = colorant"#717280"
const THEME = Theme(
    backgroundcolor = BG,
    Axis = (
        backgroundcolor = :transparent,
        bottomspinecolor = AX,
        leftspinecolor = AX,
        rightspinecolor = AX,
        topspinecolor = AX,
        xgridcolor = GR,
        xtickcolor = AX,
        xlabelcolor = TX,
        xticklabelcolor = TX,
        ygridcolor = GR,
        ytickcolor = AX,
        ylabelcolor = TX,
        yticklabelcolor = TX,
    ),
    Legend = (bgcolor = :transparent, framecolor = AX, labelcolor = TX),
    Label = (color = TX,),
)
set_theme!(THEME)
# Full names of parameters
const PARAM_MAP =
    Dict("repr" => "Reproduction Threshold", "speed" => "Speed", "sens" => "Sensory Radius")

# Plot the food as a heatmap, with a slider for interactivity
function plot_food(foodlog::String)
    fig = Figure(resolution = (600, 600))
    food = h5open(foodlog, "r")
    maxtime = size(food["log"], 1)
    time = Slider(fig[1, 1], range = 1:maxtime, startvalue = 1)
    cfood = @lift(food["log"][$(time.value), :, :])
    ax = fig[0, 1] = Axis(fig)
    heatmap!(ax, cfood, colorrange = (0, food["cap"][]))

    fig, food
end

# Record the food map playing through steps
function record_food(
    foodlog::String,
    filename::String;
    duration::Real = 10,
    nsteps::Int = -1,
)
    CairoMakie.activate!()
    fig = Figure(resolution = (600, 600), figure_padding = 0)
    food = h5open(foodlog, "r")
    maxtime = size(food["log"], 1)
    nsteps = nsteps <= 0 ? maxtime : min(nsteps, maxtime)
    ax = fig[1, 1] = Axis(fig; aspect = AxisAspect(1))
    hidedecorations!(ax)
    time = Observable(1)
    cfood = @lift(food["log"][$(time), :, :])
    titletext = @lift("Step: $($(time))")
    fig[0, 1] = Label(fig, titletext; tellwidth = false)
    heatmap!(ax, cfood, colorrange = (0, food["cap"][]))
    record(fig, filename, 1:nsteps; framerate = floor(Int, nsteps / duration)) do i
        time[] = i
    end
    close(food)
    GLMakie.activate!()
end

# Plot data of the simulation
function plot_agent_scatter(
    dir::String;
    params = ["sens", "repr", "speed"],
    fig = Figure(resolution = (1600, 900)),
    plotbact = false,
)
    agg_df = CSV.File(joinpath(dir, "agg.csv")) |> DataFrame
    n_species = maximum(agg_df.species)
    stepcount = maximum(agg_df.step)
    print(agg_df.step)
    print(stepcount)
    print(typeof(stepcount))
    agg_df = groupby(agg_df, :species; sort = true)

    plotline(ax, steps, data, color) = lines!(ax, steps, data; color, linewidth = 2)

    if plotbact
        df =
            CSV.File(joinpath(dir, "bact.csv"); select = ["step", "species", params...]) |>
            DataFrame
        df = groupby(df, :species; sort = true)
        scatter_cols = distinguishable_colors(n_species, HSV(0, 0.6, 1))
    end

    sp_cols = distinguishable_colors(n_species, HSV(0, 0.6, 1))

    ax =
        fig[1, :] = Axis(
            fig;
            ylabel = "Population",
            xticks = 0:500:stepcount,
            xticksvisible = false,
            xticklabelsvisible = false,
        )
    xlims!(ax, (0, stepcount))
    lines = [plotline(ax, agg_df[i].step, agg_df[i].nbact, sp_cols[i]) for i = 1:n_species]

    for (ind, par) in enumerate(params)
        ax =
            fig[ind+1, :] = Axis(
                fig;
                ylabel = PARAM_MAP[par],
                xticks = 0:500:stepcount,
                xticksvisible = ind == length(params),
                xticklabelsvisible = ind == length(params),
            )
        for i = 1:n_species
            xlims!(ax, (0, stepcount))
            plotbact && scatter!(
                ax,
                df[i].step,
                df[i][!, par];
                color = (scatter_cols[i], 0.01),
                markersize = 4,
                strokewidth = 0,
            )
            plotline(ax, agg_df[i].step, agg_df[i][!, "μ_$par"], sp_cols[i])
            band!(
                ax,
                agg_df[i].step,
                agg_df[i][!, "μ_$par"] .- agg_df[i][!, "σ_$par"],
                agg_df[i][!, "μ_$par"] .+ agg_df[i][!, "σ_$par"],
                color = (sp_cols[i], 0.25),
            )
        end
    end

    leg =
        fig[length(params)+2, :] = Legend(
            fig,
            lines,
            ["Mean" for _ = 1:n_species];
            orientation = :horizontal,
            tellwidth = false,
            tellheight = true,
        )
    fig
end

# Multiple `plot_agent_scatter` side by side
function comparative_agent_scatter(dirs::Vector{String}; kwargs...)
    fig = Figure(resolution = (1200, 1200))

    for (i, dir) in enumerate(dirs)
        plot_agent_scatter(dir; fig = fig[1, i], kwargs...)
    end
    fig
end

max_stars_repo_name: jiaangou/SimpleEvolution

max_issues_repo_name: jiaangou/SimpleEvolution

max_forks_repo_name: jiaangou/SimpleEvolution


text:
# Contents: Functions for hashing passwords using PBKDF2.

"""
Returns the hashed password as Array{UInt8, 1}.
"""
function pbkdf2(salt::Array{UInt8, 1}, password::AbstractString, niter::Int64, dklen::Int64)
    hlen = 64                    # Length of output of hash function in bytes (SHA512 >> length = 512 bits = 64 bytes)
    assert(dklen % hlen == 0)    # Assert dklen is divisible by hlen
    assert(niter >= 1000)
    derived_key = zeros(UInt8, dklen)
    nblocks     = convert(Int, round(dklen / hlen))
    for i = 1:nblocks
        dk_block_i = pbkdf2_block(salt, password, niter, i)    # dk_block_i has length hlen

	# Copy i^th block to result
	offset = (i - 1) * hlen
	for j = 1:hlen
	    derived_key[offset + j] = dk_block_i[j]
	end
    end
    derived_key
end

"""
Compute the i^th block of the derived key.
Uses SHA512 as the pseudorandom function.
"""
function pbkdf2_block(salt, password, niter, i)
    key    = password
    U0     = vcat(salt, num_to_bytearray(hton(Int32(i))))    # Length = 20 bytes = vcat(16 byte salt, 4 byte bigendian integer)
    U1     = digest(MD_SHA, U0, key)                         # MD_SHA == SHA512
    U2     = digest(MD_SHA, U1, key)
    result = U1 $ U2
    Ujm1   = U2    # U_(j-1)
    for j = 3:niter
	Uj      = digest(MD_SHA, Ujm1, key)
	result $= Uj
	Ujm1    = Uj
    end
    result
end

### EOF


text:
using HDF5

#to read SIFT1B
include("./bvecs_read.jl")
include("./ivecs_read.jl")
include("./fvecs_read.jl")

function read_dataset(
  dname::AbstractString,  # name of the dataset to load
  nvectors::Union{Integer, UnitRange},      # number of vectors to read
  V::Bool=false )         # whether to print progress

  if V print("Loading $(dname)... "); end

  # ====  SIFT1M ===
  if dname == "SIFT1M"

    fname = "./data/sift/sift_learn.fvecs";
    X     = fvecs_read(nvectors, fname)

  elseif dname == "SIFT1M_query"

    fname = "./data/sift/sift_query.fvecs"
    X     = fvecs_read(nvectors, fname)

  elseif dname == "SIFT1M_groundtruth"

    fname = "./data/sift/sift_groundtruth.ivecs"
    X     = ivecs_read(nvectors, fname)

  elseif dname == "SIFT1M_base"

    fname = "./data/sift/sift_base.fvecs";
    X     = fvecs_read(nvectors, fname)

  else

    error("Dataset $(dname) unknown")

  end

  if V println("done."); end
  return X

end

max_issues_repo_name: xiaosu-zhu/local-search-quantization

max_forks_repo_name: xiaosu-zhu/local-search-quantization


text:
module day8

function readInput()
    result = NamedTuple{(:op, :x),Tuple{Symbol,Int64}}[]
    for line in eachline(joinpath(@__DIR__, "input.txt"))
        op, x = split(line)
        push!(result, (op = Symbol(op), x = parse(Int, x)))
    end
    result
end

function run(input)
    acc, ptr = 0, 1
    states = Set{Int}()
    terminate = false
    while true
        ptr in states && return acc, false
        ptr == length(input) && (terminate = true)
        !(1 <= ptr <= length(input)) && return typemin(Int), false
        push!(states, ptr)
        instruction = input[ptr]
        if instruction.op === :acc
            acc += instruction.x
            ptr += 1
        elseif instruction.op === :jmp
            ptr += instruction.x
        elseif instruction.op === :nop
            ptr += 1
        else
            error("unknown instruction $(instruction.op)")
        end
        terminate && break
    end
    return acc, true
end

function part1(input)
    acc, terminated = run(input)
    return acc
end

function part2(input)
    for i in findall(i -> i.op !== :acc, input)

        if input[i].op === :jmp
            input[i] = (op = :nop, x = input[i].x)
        else
            input[i] = (op = :jmp, x = input[i].x)
        end

        acc, success = run(input)
        success && return acc
        if input[i].op === :jmp
            input[i] = (op = :nop, x = input[i].x)
        else
            input[i] = (op = :jmp, x = input[i].x)
        end

    end
end

end


text:
# ---
# title: 1657. Determine if Two Strings Are Close
# id: problem1657
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Greedy
# link: <https://leetcode.com/problems/determine-if-two-strings-are-close/description/>
# hidden: true
# ---
#
# Two strings are considered **close** if you can attain one from the other
# using the following operations:
#
#   * Operation 1: Swap any two **existing** characters.
#     * For example, `a _b_ cd _e_ -> a _e_ cd _b_`
#   * Operation 2: Transform **every** occurrence of one **existing** character into another **existing** character, and do the same with the other character.
#     * For example, `_aa_ c _abb_ -> _bb_ c _baa_` (all `a`'s turn into `b`'s, and all `b`'s turn into `a`'s)
#
# You can use the operations on either string as many times as necessary.
#
# Given two strings, `word1` and `word2`, return `true` _if_`word1` _and_`word2`
# _are **close** , and _`false` _otherwise._
#
#
#
# **Example 1:**
#
#
#
#     Input: word1 = "abc", word2 = "bca"
#     Output: true
#     Explanation: You can attain word2 from word1 in 2 operations.
#     Apply Operation 1: "a _bc_ " -> "a _cb_ "
#     Apply Operation 1: " _a_ c _b_ " -> " _b_ c _a_ "
#
#
# **Example 2:**
#
#
#
#     Input: word1 = "a", word2 = "aa"
#     Output: false
#     Explanation: It is impossible to attain word2 from word1, or vice versa, in any number of operations.
#
#
# **Example 3:**
#
#
#
#     Input: word1 = "cabbba", word2 = "abbccc"
#     Output: true
#     Explanation: You can attain word2 from word1 in 3 operations.
#     Apply Operation 1: "ca _b_ bb _a_ " -> "ca _a_ bb _b_ "
#     Apply Operation 2: " _c_ aa _bbb_ " -> " _b_ aa _ccc_ "
#     Apply Operation 2: " _baa_ ccc" -> " _abb_ ccc"
#
#
# **Example 4:**
#
#
#
#     Input: word1 = "cabbba", word2 = "aabbss"
#     Output: false
#     Explanation: It is impossible to attain word2 from word1, or vice versa, in any amount of operations.
#
#
#
#
# **Constraints:**
#
#   * `1 <= word1.length, word2.length <= 105`
#   * `word1` and `word2` contain only lowercase English letters.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
function retro_substituicao(A, b)
    n = size(A)[2]
    x = zeros(n)
    x[n] = b[n]/A[n,n]
    for i = n-1:-1:1
        temp = 0
        for k = i+1:n
            temp += A[i,k]*x[k]
        end
        x[i] = (b[i]-temp)/A[i,i]
    end
    return x
end

function substituicao(A, b)
    n = size(A)[2]
    x = zeros(n)
    x[1] = b[1]/A[1,1]
    for i = 2:n
        temp = 0
        for k = 1:i-1
            temp += A[i,k]*x[k]
        end
        x[i] = (b[i]-temp)/A[i,i]
    end
    return x
end

function eliminacao(A)
    row, col = size(A)
    for j in 1:row-1
        for i in j+1:row
            m = A[i,j]/A[j,j]
            A[i,:] = A[i,:]-m*A[j,:]
        end
    end
    return A
end

function fatoracao_lu(A)
    row, col = size(A)
    L = eye(row)

    for j in 1:row-1
        for i in j+1:row
            L[i,j] = A[i,j]/A[j,j]
            A[i,:] = A[i,:]-L[i,j]*A[j,:]
        end
    end

    U = A[:, 1:end-1]
    return L, A
end

function resolve_fatlu(A)
    for i in 6:8
        L, U = fatoracao_lu(A[:, [1:5; i]])
        println("A fatoração LU estará coreta se L*U for igual a A ")
        println(L*U, A)

        y = substituicao(L, b)
        x = retro_substituicao(U, y)
        println("Valor do y é $y o valor do x é $x")
    end
end

max_stars_repo_path: 003 - Sistemas_Lineares.jl

max_stars_repo_name: LudwigFonseca/Julia-Language-Examples

max_issues_repo_path: 003 - Sistemas_Lineares.jl

max_issues_repo_name: LudwigFonseca/Julia-Language-Examples

max_forks_repo_path: 003 - Sistemas_Lineares.jl

max_forks_repo_name: LudwigFonseca/Julia-Language-Examples


text:
## functions that work for all types but return output depending on
## the input type: Timematr for Timematr, Timenum for Timenum, ...

##########
## hcat ##
##########

## hcat requires completely equal indices and equal types
## hcat preserves type
import Base.hcat
macro pres_hcat(myType)
    esc(quote
        function hcat(inst::$(myType), inst2::$(myType))
            ## check for equal indices
            if idx(inst) != idx(inst2)
                error("indices must coincide for hcat")
            end

            instNew = $(myType)(hcat(inst.vals, inst2.vals), idx(inst))
            return instNew
        end
    end)
end

for t = (:Timedata, :Timenum, :Timematr)
    eval(macroexpand(:(@pres_hcat($t))))
end

macro pres_hcat_varargs(myType)
    esc(quote
        function hcat(tds::$(myType)...)
            nArgs = length(tds)

            res = deepcopy(tds[1])
            for ii=2:nArgs
                res = hcat(res, tds[ii])
            end
            res
        end
    end)
end

for t = (:Timedata, :Timenum, :Timematr)
    eval(macroexpand(:(@pres_hcat_varargs($t))))
end

function hcat(inst1::AbstractTimedata, inst2::AbstractTimedata)
    ## not defined for different types
    error("objects must be of equal type for hcat")
end

function hcat(inst::AbstractTimedata...)
    ## not defined for different types
    error("objects must be of equal type for hcat")
end

##########
## vcat ##
##########

import Base.vcat
macro pres_vcat(myType)
    esc(quote

        function vcat(inst::$(myType), inst2::$(myType))
            ## check for equal names
            if names(inst) != names(inst2)
                error("variable names must coincide for vcat")
            end

            if idxtype(inst) != idxtype(inst2)
                error("time index must be of same type for vcat")
            end

            catVals = vcat(inst.vals, inst2.vals)
            return $(myType)(catVals, [idx(inst), idx(inst2)])
        end
    end)
end

for t = (:Timedata, :Timenum, :Timematr)
    eval(macroexpand(:(@pres_vcat($t))))
end

macro pres_vcat_varargs(myType)
    esc(quote
        function vcat(tds::$(myType)...)
            ## all components must have equal names
            colnames = names(tds[1])
            nTds = length(tds)

            res = deepcopy(tds[1])
            for ii=2:nTds
                res = vcat(res, tds[ii])
            end
            res
        end
    end)
end

for t = (:Timedata, :Timenum, :Timematr)
    eval(macroexpand(:(@pres_vcat_varargs($t))))
end

function vcat(inst1::AbstractTimedata, inst2::AbstractTimedata)
    ## not defined for different types
    error("objects must be of equal type for vcat")
end

function vcat(inst::AbstractTimedata...)
    ## not defined for different types
    error("objects must be of equal type for vcat")
end

############
## flipud ##
############

import Base.flipud

macro pres_flipud(myType)
    esc(quote
        function flipud(inst::$(myType))
            ## flip data upside down

            flippedVals = DataFrame()
            for colName in names(inst.vals)
                flippedVals[colName] = flipud(inst.vals[colName])
            end
            dats = flipud(idx(inst))
            return $(myType)(flippedVals, dats)
        end
    end)
end

for t = (:Timedata, :Timenum, :Timematr)
    eval(macroexpand(:(@pres_flipud($t))))
end

##########
## narm ##
##########

function narm(td::AbstractTimedata)
    ## return days without NA
    smallTd = td[complete_cases(td), :]
    return smallTd
end

function narm(tm::AbstractTimematr)
    ## no NAs allowed
    return deepcopy(tm)
end

########################
## asTd / asTm / asTn ##
########################

function isaRowVector(x::Array)
    isa(x, Matrix) & (size(x, 1) == 1 )
end

function isaRowVector(x::Any)
    error("isaRowVector only is defined for Arrays")
end

function asArrayOfEqualDimensions(arr::Array, td::AbstractTimedata)
    if !isa(arr, Vector) & !isaRowVector(arr)
        error("array must be either row or column vector for
    conversion")
    end

    (nObs, nVars) = size(td)
    nElem = length(arr)

    if isaRowVector(arr)
        # interpretation as row vector
        if nElem != nVars
            error("wrong dimensions: number of columns not equal")
        end
        resVals = repmat(arr, nObs, 1)
    else # column vector
        if nElem == nObs
            resVals = repmat(arr, 1, nVars)
        elseif nElem == nVars
            resVals = repmat(arr', nObs, 1)
        else
            error("wrong dimensions: array can not be converted")
        end
    end
    resVals
end

function asArrayOfEqualDimensions(val::Any, td::AbstractTimedata)
    (nObs, nVars) = size(td)

    resVals = repmat([val], nObs, nVars)
    return resVals
end

function asTd(arr::Array, td::Timedata)
    resVals = asArrayOfEqualDimensions(arr, td)
    td = Timedata(resVals, names(td), idx(td))
end

function asTn(arr::Array, td::Timenum)
    resVals = asArrayOfEqualDimensions(arr, td)
    td = Timenum(resVals, names(td), idx(td))
end

function asTm(arr::Array, td::Timematr)
    resVals = asArrayOfEqualDimensions(arr, td)
    td = Timematr(resVals, names(td), idx(td))
end


text:
using MathOptInterface, LinearAlgebra, Statistics, Plots, QuantEcon, Interpolations, NLsolve, Optim, Random, IterTools, JuMP, Ipopt, MathOptInterface;

u(x) = sum(log(x[i]) for i in 1:numSectors)
w(x) = sum(5*log(x[i]) for i in 1:numSectors)

β = 0.96
α = 0.4
f(x) = x^α

numSectors = 2

numPoints1D = 2

grid = Vector{Vector{Float64}}(undef,(numPoints1D)^numSectors)

gridMax = 3
gridMin = 1

iter=1
if numSectors == 1
	grid = LinRange(gridMin, gridMax, numPoints1D)
elseif numSectors == 2
for p in product(LinRange(gridMin,gridMax,numPoints1D),LinRange(gridMin,gridMax,numPoints1D))
    grid[iter] = collect(p)
    global iter += 1
end
else
	print("Currently we only allow one or two sectors")
end
wVal = w.(grid)

function interp(x,y)
    return interpolate((LinRange(gridMin,gridMaxnumPoints1D),LinRange(gridMin,gridMax,numPoints1D)), wVals, Gridded(Linear()))(x,y)
end

function T(wVal, grid, β, f ; compute_policy = false)
	if numSectors == 1
		wVals = zeros(numPoints1D)
		wVals = wVal
		wFunc(y) = interpolate(grid)
	elseif numSectors == 2
		wVals = zeros(numPoints1D,numPoints1D);
		for j in numPoints1D
		    for i in numPoints1D
                        wVals[i,j] = wVal[i+(i-1)*j]
                    end
                end
		wFunc(x,y) = log(x) + log(y)
		#wFunc(x, y) = interpolate((LinRange(gridMin,gridMax,numPoints1D),LinRange(gridMin,gridMax,numPoints1D)),			  wVals, Gridded(Linear()))(x,y)
		print("wFunc created")
		print(wFunc(1.5,1.5))
	else
		print("Can only handle upto 2 sectors")
	end
    print(wFunc(1.5, 1.5))
    global Tw = zeros(length(grid))
    global σ = similar(grid)
    for n in 1:length(grid)
        y = grid[n]
        modTrial = Model(Ipopt.Optimizer);
        @variable(modTrial,  c[1:numSectors] >= 0.0001)
        @variable(modTrial, k[1:numSectors] >= 0.0001, start=2)
        for i in 1:numSectors
            @constraint(modTrial, gridMin <= c[i] <= y[i])
	    @constraint(modTrial, k[i] == y[i] - c[i])
        end
        register(modTrial, :wFunc, numSectors, wFunc; autodiff = true)
        register(modTrial, :f, 1, f; autodiff = true)
	if numSectors == 1
		@NLobjective(modTrial, Max, sum(log(c[i]) for i in 1:numSectors) + β*wFunc(f(k[1])))
	elseif numSectors == 2
        @NLobjective(modTrial, Max, sum(log(c[i]) for i in 1:numSectors) + β*wFunc(f(k[1]),f(k[2])))
	else
	print("check the number of sectors")
	end
        optimize!(modTrial)
        Tw[n] = JuMP.objective_value(modTrial)
        if compute_policy
            σ[n] = value.(c)
        end
    end
    if compute_policy
        return Tw, σ
    end

    return Tw

end

wVal = T(wVal, grid, β, f; compute_policy = false)
#=
u(c) = log(c[1]^0.5*c[2]^0.5);

numSectors == 2;

w(x,y) = x^2+y^2;

valGrid = zeros(numPoints1D,numPoints1D)
for i in 1:numPoints1D
    for j in 1:numPoints1D
        valGrid[i,j] = w(grid[j+(i-1)*numPoints1D][1],grid[j+(i-1)*numPoints1D][2])
    end
end

valGrid = zeros(numPoints1D,numPoints1D)
for i in 1:numPoints1D
    for j in 1:numPoints1D
        valGrid[i,j] = w(grid[j+(i-1)*numPoints1D][1],grid[j+(i-1)*numPoints1D][2])
    end
end

α = 0.4
β = 0.96
μ = 0
s = 0.1
ϵ_D = 1
ϵ_X = 0.1
Γ = [0.7 0.3; 0.3 0.7]

#initial capital flows
x= [1 2; 3 4]
#Investment is a function from R^2J to R^J
X = zeros(numSectors,1)
XInter = zeros(numSectors,numSectors)
for J in numSectors;
    for I in numSectors;
        XInter[I,J] = (Γ[I,J]^(1/ϵ_X)*x[I,J]^((ϵ_X-1)/ϵ_X))
    end
    X[J] = sum(XInter[:,J])^(ϵ_X/(ϵ_X-1))
end

K = zeros(numSectors,1)
Q = zeros(numSectors,1)
C = zeros(numSectors,1)

#Kapital is function from R^J to R^J
K = X;

Q = K.^α

C = Q - sum(x[:,J] for J in 1:numSectors)

#=
c1 = log(1 - α * β) / (1 - β)
c2 = (μ + α * log(α * β)) / (1 - α)
c3 = 1 / (1 - β)
c4 = 1 / (1 - α * β)
=#

# Utility is a function from R^J to R
u(C) = log(sum(C[j]^((ϵ_D-1)/ϵ_D) for j in 1:numSectors)^((ϵ_D)/(ϵ_D*0.999-1)));

# Deterministic part of production function a function from R^J to R^J
f(k) = k.^α

#= True optimal policy
c_star(y) = (1 - α * β) * y

# True value function
v_star(y) = c1 + c2 * (c3 - c4) + c4 * log(y)
=#

Random.seed!(42) # For reproducible results.

#=
grid = Array{Float64}[]
for n in 1:numSectors
    x = 1:numSectors
    push!(grid, x)
end
=#

#=
grid[j] = 1:numSectors
y = 1:numSectors
x' .* ones(5)
ones(3)' .* y
g         # Largest grid point
grid_size = 200      # Number of grid points
shock_size = [250,250]     # Number of shock draws in Monte Carlo integral

grid_y = range.(1e-5,  gridMax, length = grid_size)
shocks = exp.(μ .+ s * randn(shock_size))
w = T(v_star.(grid_y), grid_y, β, (log∘sum), k -> k^α, shocks)

#=
plt = plot(ylim = (-35,-24))
plot!(plt, grid_y, w, linewidth = 2, alpha = 0.6, label = "T(v_star)")
plot!(plt, v_star, grid_y, linewidth = 2, alpha=0.6, label = "v_star")
plot!(plt, legend = :bottomright)
=#
w = 5 * log.(grid_y)  # An initial condition -- fairly arbitrary
n = 35
#=
plot(xlim = (extrema(grid_y)), ylim = (-50, 10))
lb = "initial condition"
plt = plot(grid_y, w, color = :black, linewidth = 2, alpha = 0.8, label = lb)
for i in 1:n
    global w = T(w, grid_y, β, log, k -> k^α, shocks)
    plot!(grid_y, w, color = RGBA(i/n, 0, 1 - i/n, 0.8), linewidth = 2, alpha = 0.6,
          label = "")
end

lb = "true value function"
plot!(plt, v_star, grid_y, color = :black, linewidth = 2, alpha = 0.8, label = lb)
plot!(plt, legend = :bottomright)
=#

function solve_optgrowth(initial_w; tol = 1e-6, max_iter = 500)
    fixedpoint(w -> T(w, grid_y, β, u, f, shocks), initial_w).zero # gets returned
end

#=
initial_w = 5 * log.(grid_y)
v_star_approx = solve_optgrowth(initial_w)

plt = plot(ylim = (-35, -24))
plot!(plt, grid_y, v_star_approx, linewidth = 2, alpha = 0.6,
      label = "approximate value function")
plot!(plt, v_star, grid_y, linewidth = 2, alpha = 0.6, label = "true value function")
plot!(plt, legend = :bottomright)

Tw, σ = T(v_star_approx, grid_y, β, log, k -> k^α, shocks;
                         compute_policy = true)
cstar = (1 - α * β) * grid_y

plt = plot(grid_y, σ, lw=2, alpha=0.6, label = "approximate policy function")
plot!(plt, grid_y, cstar, lw = 2, alpha = 0.6, label = "true policy function")
plot!(plt, legend = :bottomright)
=#
=#
=#

max_stars_repo_name: uq-aibe/aim

max_issues_repo_name: uq-aibe/aim

max_forks_repo_name: uq-aibe/aim


text:
# create a timer object, used throughout this package, see TimerOutputs.jl
global const to = TimerOutput();

# struct to store FD params that are used in @initialize
mutable struct FD_params
    ndims::Int64 # 2D or 3D
    order::Int64 # stencil order
    use_gpu::Bool #
    datatype::DataType
    npml::Int64
    npextend::Int64
end
# dummy finite-difference stencil order; we need to initialize anyway
const _fd = FD_params(2, 2, false, Float32, 20, 20+2)

"""
Return axis names of 1D, 2D or 3D fields
"""
function dim_names(N, prefix = "", suffix = ""; order = 1)
    if (N == 3)
        names = (order == 2) ? [:zz, :yy, :xx, :xz, :xy, :yz] : [:z, :y, :x]
    elseif (N == 2)
        names = (order == 2) ? [:zz, :xx, :xz] : [:z, :x]
    elseif (N == 1)
        names = [:z]
    else
        error("invalid dim num")
    end
    return [Symbol(prefix, string(m), suffix) for m in names]
end


text:
module NeuralynxNCS
using PiecewiseIncreasingRanges, Compat

immutable NCSContinuousChannel
    header::ByteString
    samples::Vector{Int16}
    times::PiecewiseIncreasingRange{Float64,StepRange{Int64,Int64},Int64}
end

immutable NCSRecordHeader
    qwTimeStamp::UInt64
    dwChannelNumber::UInt32
    dwSampleFreq::UInt32
    dwNumValidSamples::UInt32
end

Base.read(io::IO, ::Type{NCSRecordHeader}) =
    NCSRecordHeader(read(io, UInt64), read(io, UInt32), read(io, UInt32), read(io, UInt32))

function compute_times(rec::NCSRecordHeader, sample_multiplier::Int64, step::Int64)
    @compat first_sample = Int64(rec.qwTimeStamp)*sample_multiplier
    first_sample:step:first_sample+step*(rec.dwNumValidSamples-1)
end

function readncs(filename::String; kwargs...)
    io = open(filename, "r")
    try
        readncs(io, filesize(io); kwargs...)
    finally
        close(io)
    end
end

function readncs(io::IO, sz::Union(Integer,Nothing)=nothing; skip_nonmonotonic::Bool=false)
    header = rstrip(bytestring(read(io, UInt8, 16384)), '\0')
    eof(io) && return NCSContinuousChannel(header, Int16[], PiecewiseIncreasingRange(StepRange{Int,Int}[], 1))

    nrecs::Int = isa(sz, Integer) ? div(sz-16384, 1044) : -1
    sample_buffer = Array(Int16, 512)
    samples = Array(Int16, isa(sz, Integer) ? nrecs*512 : 512)
    nsamples = 0
    times = Array(StepRange{Int,Int}, isa(sz, Integer) ? nrecs : 1)

    # Read first header so we can save its channel number to make
    # sure we only have one
    rec1 = read(io, NCSRecordHeader)
    read!(io, sample_buffer)
    copy!(samples, nsamples+1, sample_buffer, 1, rec1.dwNumValidSamples)
    nsamples += rec1.dwNumValidSamples

    # Figure out range multipliers

    # This would be correct if the reported sample rate were:
    # @compat divisor = lcm(Int64(rec1.dwSampleFreq), 10^6)
    # sample_multiplier = div(divisor, 10^6)
    # @compat step = div(divisor, Int64(rec1.dwSampleFreq))

    # However, it seems that the reported sample rate is wrong. The
    # actual sample rate is rounded so that the delta between
    # timestamps is constant.
    sample_multiplier = 512
    divisor = 512*10^6
    step = div(divisor, rec1.dwSampleFreq)

    nrecs = 1
    reccount = 1
    times[1] = compute_times(rec1, sample_multiplier, step)

    while !eof(io)
        reccount += 1
        rec = read(io, NCSRecordHeader)
        rec.dwChannelNumber == rec1.dwChannelNumber || error("only one channel supported per file")
        rec.dwSampleFreq == rec1.dwSampleFreq || error("sample rate is non-constant")
        read!(io, sample_buffer)
        rectimes = compute_times(rec, sample_multiplier, step)
        if first(rectimes) < last(times[nrecs])
            if skip_nonmonotonic
                warn("""timestamp at record $(reccount) is not monotonically increasing
                        (previous end time = $(last(times[nrecs])/divisor), start time = $(first(rectimes)/divisor); skipping""")
                continue
            else
                error("""timestamp at record $(reccount) is not monotonically increasing.
                         Use readncs(..., skip_nonmotonoic=true) to skip nonmonotonic timestamps""")
            end
        end
        isa(sz, Nothing) && resize!(samples, nsamples+rec.dwNumValidSamples)
        copy!(samples, nsamples+1, sample_buffer, 1, rec.dwNumValidSamples)
        nsamples += rec.dwNumValidSamples
        nrecs += 1
        if isa(sz, Integer)
            times[nrecs] = rectimes
        else
            push!(times, rectimes)
        end
    end

    resize!(samples, nsamples)
    resize!(times, nrecs)

    return NCSContinuousChannel(header, samples, PiecewiseIncreasingRange(times, divisor))
end

export readncs
end


text:
using Flux
using Transformers.Basic

struct Decoder
    layers::Array{TransformerDecoder}
end

Flux.@functor Decoder
Flux.trainable(decoder::Decoder) = (decoder.layers,)

function (decoder::Decoder)(target_embedding, encoded_embedding)
    embedding = target_embedding
    for transformer ∈ decoder.layers
        embedding = transformer(embedding, encoded_embedding)
    end
    return embedding
end

function Decoder(
    size::Integer,
    head::Integer,
    hs::Integer,
    ps::Integer,
    layer::Integer;
    act=relu,
    pdrop=0.1
)::Decoder
    return Decoder([
        TransformerDecoder(size, head, hs, ps, act=act, pdrop=pdrop)
        for i ∈ 1:layer
    ])
end


text:
# 32-bit
# pays a heavy cost for the array indexing
@inline function load32u(a::Vector{UInt8}, i::Integer)
    local val::UInt32 = (0x00000000
        | UInt32(a[i])
        | UInt32(a[i+1]) << 8
        | UInt32(a[i+2]) << 16
        | UInt32(a[i+3]) << 24)
    return val
end

@inline function store32u!(a::Vector{UInt8}, i::Integer, u::UInt32)
    a[i] = u % UInt8
    a[i+1] = (u >> 8) % UInt8
    a[i+2] = (u >> 16) % UInt8
    a[i+3] = (u >> 24) % UInt8
end
@inline function unaligned_copy_64u!(dst::Vector{UInt8}, di::Integer, src::Vector{UInt8}, si::Integer)
    for i in 0:7
        dst[di+i] = src[si+i]
    end
end
@inline function unaligned_copy_128u!(dst::Vector{UInt8}, di::Integer, src::Vector{UInt8}, si::Integer)
    for i in 0:15
        dst[di+i] = src[si+i]
    end
end


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/FEMBasis.jl/blob/master/LICENSE

function subs(p::Number, ::Any)
    return p
end

function subs(p::Symbol, data::Pair{Symbol, T}) where T
    k, v = data
    if p == k
        return v
    end
    return p
end

function subs(p::Symbol, data::NTuple{N,Pair{Symbol, T}}) where {N, T}
    for (k, v) in data
        if p == k
            return v
        end
    end
    return p
end

function subs(p::Expr, d::Pair)
    v = copy(p)
    for j in 2:length(p.args)
        v.args[j] = subs(v.args[j], d)
    end
    return v
end

"""
    subs(expression, data)

Given expression and pair(s) of `symbol => value` data, substitute to expression.

# Examples

Let us have polynomial `1 + u + v + u*v^2`, and substitute u=1 and v=2:
```julia
expression = :(1 + u + v + u*v^2)
data = (:u => 1.0, :v => 2.0)
subs(expression, data)
8.0
```

"""
function subs(p::Expr, data::NTuple{N,Pair{Symbol, T}}) where {N, T}
    for di in data
        p = subs(p, di)
    end
    return Calculus.simplify(p)
end


text:
function simguidedlm_llikelihoodθ!(::LeftRule,  Xᵒ, x0, W, Q::GuidedProposal!, k; skip = 0, ll0 = true)
    Pnt = eltype(x0)
    tt =  Xᵒ.tt
    Xᵒ.yy[1] .= deepvalue(x0)

    typeθel = typeof(Q.target.a)
    som::typeθel  = 0.  # here is the only change!!!

    dump(typeθel)
    # initialise objects to write into
    # srout and strout are vectors of Points
    dwiener = dimwiener(Q.target)
    srout = zeros(Pnt, dwiener)
    strout = zeros(Pnt, dwiener)

    #x = copy(x0)  ##### ADJUST: this should also be of dual type !!!
    n = Q.target.n
    x = NState(zeros(Point{typeθel},n),zeros(Point{typeθel},n))

    xarb = NState(zeros(Point{typeθel},n),zeros(Point{typeθel},n))
    rout = copy(xarb)
    bout = copy(xarb)
    btout = copy(xarb)
    wout = copy(xarb)

    if !constdiff(Q)
        At = Bridge.a((1,0), x0, auxiliary(Q,k))  # auxtimehomogeneous switch
        A = zeros(Unc{deepeltype(x0)}, 2Q.target.n,2Q.target.n)
    end

    for i in 1:length(tt)-1
        dt = tt[i+1]-tt[i]
        #println(typeof(bout)==typeof(x))
        b!(tt[i], x, bout, target(Q)) # b(t,x)
        _r!((i,tt[i]), x, rout, Q,k) # tilder(t,x)
        σt!(tt[i], x, rout, srout, target(Q))      #  σ(t,x)' * tilder(t,x) for target(Q)
        Bridge.σ!(tt[i], x, srout*dt + W.yy[i+1] - W.yy[i], wout, target(Q)) # σ(t,x) (σ(t,x)' * tilder(t,x) + dW(t))
        # likelihood terms
        if i<=length(tt)-1-skip
            _b!((i,tt[i]), x, btout, auxiliary(Q,k))
            som += dot(bout-btout, rout) * dt
            if !constdiff(Q)
                σt!(tt[i], x, rout, strout, auxiliary(Q,k))  #  tildeσ(t,x)' * tilder(t,x) for auxiliary(Q)
                som += 0.5*Bridge.inner(srout) * dt    # |σ(t,x)' * tilder(t,x)|^2
                som -= 0.5*Bridge.inner(strout) * dt   # |tildeσ(t,x)' * tilder(t,x)|^2
                Bridge.a!((i,tt[i]), x, A, target(Q))
                som += 0.5*(dot(At,Q.guidrec.Ht[i]) - dot(A,Q.guidrec.Ht[i])) * dt
            end
        end
        x .= x + dt * bout + wout
        Xᵒ.yy[i+1] .= deepvalue(x)
    end
    if ll0
        logρ0 = lρtilde(x0,Q,k)
    else
        logρ0 = 0.0 # don't compute
    end
    copyto!(Xᵒ.yy[end], Bridge.endpoint(Xᵒ.yy[end],Q))
    som + logρ0
end

function simguidedlm_llikelihoodθ!(::LeftRule,  X, x0, W, Q::GuidedProposal!; skip = 0, ll0 = true) # rather would like to dispatch on type and remove '_mv' from function name
    soms  = zeros(typeof(Q.target.a), Q.nshapes)
    for k in 1:Q.nshapes
        soms[k] = simguidedlm_llikelihoodθ!(LeftRule(), X[k],x0,W[k],Q,k ;skip=skip,ll0=ll0)
    end
    soms
end

# attempt to include automatic differentiation on parameter updates

function slogρ!(x0deepv, θ, Q, W, X, obs_info, llout)
    #(a, c, γ) = θ
    a = θ[1]; c= θ[2]; γ = θ[3]
    println("θ in slogρ! ")
    dump(a)
    P = Q.target
    if isa(P,MarslandShardlow)
        #Ptemp = MarslandShardlow(a, c, γ, P.λ, P.n)
        Ptemp = MarslandShardlow(a, c, γ, 0.0a, P.n) # be careful, only works if λ=0 is assumed
    elseif isa(P,Landmarks)
        nfs = construct_nfs(P.db, P.nfstd, γ)
        Ptemp = Landmarks(a, c, P.n, P.db, P.nfstd, nfs)
    end
    #Paux = [auxiliary(Ptemp,Q.xobsT[k]) for k in 1:Q.nshapes]

    Ptempaux = [auxiliary(Ptemp,State(Q.xobsT[k],Q.mT)) for k in 1:Q.nshapes] # auxiliary process for each shape
    Qtemp = GuidedProposal!(Ptemp,Ptempaux,Q.tt,Q.guidrec,Q.xobs0,Q.xobsT,Q.nshapes,Q.mT)
    update_guidrec!(Qtemp, obs_info)   # compute backwards recursion

    x0 = deepvec2state(x0deepv)
    lltemp = simguidedlm_llikelihoodθ!(LeftRule(), X, x0, W, Qtemp; skip=sk)  #overwrites X

    llout .= ForwardDiff.value.(lltemp)
    sum(lltemp)
end
slogρ!(x0deepv, Q, W, X, obs_info, llout) = (θ) -> slogρ!(x0deepv,θ, Q, W,X,obs_info,llout)

llout = copy(ll)
θ = [0.2, 0.4, 1.2]
slogρ!(deepvec(x0), θ, Q, W, X, obs_info, llout)

∇θ = zeros(3)
ForwardDiff.gradient!(∇θ, slogρ!(deepvec(x0), Q, W, X,obs_info,llout),θ) # X gets overwritten but does not change

ForwardDiff.gradient(slogρ!(deepvec(x0), Q, W, X,obs_info,llout),θ) # X gets overwritten but does not change


text: print_range(r) = ( for i in r; print(i, " "); end; println())

max_stars_repo_name: JuliaTagBot/JuliaBook-Samples

max_issues_repo_name: JuliaTagBot/JuliaBook-Samples

max_forks_repo_name: JuliaTagBot/JuliaBook-Samples


text:
##### Beginning of file

Test.@test( isdir(PredictMDFull.package_directory()) )

Test.@test( isdir(PredictMDFull.package_directory("ci")) )

Test.@test( isdir(PredictMDFull.package_directory("ci", "travis")) )

Test.@test( isdir(PredictMDFull.package_directory(PredictMDFullTestModuleA)) )

Test.@test( isdir(PredictMDFull.package_directory(PredictMDFullTestModuleB)) )

Test.@test( isdir( PredictMDFull.package_directory(
            PredictMDFullTestModuleB, "directory2",
            ) ) )

Test.@test( isdir( PredictMDFull.package_directory(
            PredictMDFullTestModuleB, "directory2", "directory3",
            ) ) )

Test.@test_throws(
    ErrorException,
    PredictMDFull.package_directory(PredictMDFullTestModuleC),
    )

##### End of file