text:
using Optim

""" Create a new model of the same type as `model`, but with the given `params`.
This is meant to be used with Optim.jl. Inspired from sklearn's `set_params`.

This is an internal function."""
function set_params(model::Model, params::AbstractVector, names=fieldnames(typeof(model)))
    # Not efficient, but doesn't really have to be for significant input length.
    i = 1
    upd = Dict()
    for name in names
        v = getfield(model, name)
        nvals = length(v)
        upd[name] = params[v isa Number ? i : (i:i+nvals-1)]
        i += nvals
    end
    kwargs = map(fieldnames(typeof(model))) do f
        f=>get(upd, f) do
            getfield(model, f)
        end
    end
    return roottypeof(model)(; kwargs...)
end

""" Do one round of flattening, eg.

```
julia> struct MyModel <: MiniKalman.Model  # of course, should use @with_kw
           a
           b
       end

julia> MiniKalman.get_params(MyModel(2, [3,4]))
3-element Array{Int64,1}:
 2
 3
 4
```
"""
get_params(model::Model, names=fieldnames(typeof(model))) =
    [x for v in names for x in getfield(model, v)]

################################################################################

split_units(vec::Vector) = ustrip.(vec), unit.(vec)

# get_inputs is to handle `inputs` that are Functions. It's useful when the inputs
# includes model-parameter-dependent quantities (eg. for speed). It's not a very
# nice separation of concerns, though. There might be a better way to arrange things.
get_inputs(model, inputs::Function) = inputs(model)
get_inputs(model, inputs) = inputs

""" Finds a set of model parameters that attempts to maximize the log-likelihood
on the given dataset. Returns `(best_model, optim_output_object)`. """
function Optim.optimize(model0::Model, inputs,
                        observations::Union{Nothing, AbstractVector}=nothing;
                        initial_state=MiniKalman.initial_state(model0),
                        min=0.0, # 0.0 is arbitrary... see below
                        parameters_to_optimize=parameters(model0),
                        method=LBFGS(linesearch=Optim.LineSearches.BackTracking()),
                        kwargs...)
    # It would be nice not to need split_units
    initial_x, units = split_units(get_params(model0, parameters_to_optimize))
    function objective(params)
        model = set_params(model0, params .* units, parameters_to_optimize)
        return -log_likelihood(model, get_inputs(model, inputs), observations;
                               initial_state=initial_state)
    end
    td = OnceDifferentiable(objective, initial_x; autodiff=:forward)
    # Tell Optim that no parameter can be below `mins`. This should be made optional, TODO
    mins = min isa AbstractVector ? min : fill(min, length(initial_x))
    maxes = fill(Inf, length(initial_x))
    o = optimize(td, mins, maxes, initial_x, Fminbox(method), Optim.Options(; kwargs...))
    best_model = set_params(model0, Optim.minimizer(o) .* units, parameters_to_optimize)
    return (best_model, o)
end


text:
using ProgressMeter
using SQLite
using Tables
using Dates

SQLITE_FILE = "socat.sqlite"

DJF = 0
MAM = 1
JJA = 2
SON = 3

function main()
  db = initdb()

  # Open the file and skip the header
  infile = open("SOCATv2021.tsv", "r")
  prog = ProgressUnknown("Lines read:")

  # Skip headers
  colheadersfound = false
  while !colheadersfound
    line = readline(infile)
    ProgressMeter.next!(prog)
    colheadersfound = startswith(line, "Expocode\tversion\tSource_DOI")
  end

  # Start a transaction
  SQLite.transaction(db)

  eof = false
  while !eof

    # Read the next line
    line = readline(infile)
    if line == ""
      # If there's no next line, we've finished. Write the last series
      eof = true
    else

      fields = split(line, "\t")

      lonindex = floor(parse(Float64, fields[11])) + 1
      if lonindex == 361
        lonindex = 360
      end

      latindex = floor(parse(Float64, fields[12]) + 90) + 1
      if latindex == 181
        latindex = 180
      end

      expocode = fields[1]
      year = parse(UInt16, fields[5])
      month = parse(UInt8, fields[6])
      day = parse(UInt8, fields[7])
      hour = parse(UInt8, fields[8])
      minute = parse(UInt8, fields[9])
      second = floor(parse(Float64, fields[10]))
      if second == 60
        second = 59
      end

      linetime = Dates.DateTime(
        year,
        month,
        day,
        hour,
        minute,
        second,
      )

      doy = dayofyear(linetime)

      if isleapyear(linetime) && ((month == 2 && day == 29) || month > 2)
        doy -= 1
      end

      season = getseason(doy)

      minofday = minuteofday(linetime)

      sst = parse(Float64, fields[15])
      sss = parse(Float64, fields[14])
      fco2 = parse(Float64, fields[30])

      SQLite.execute(db,
      "INSERT INTO socat VALUES " *
      "($lonindex, $latindex, $year, $doy, $season, $minofday, '$expocode', " *
      "$(isnan(sst) ? "NULL" : sst), $(isnan(sss) ? "NULL" : sss), " *
      "$(isnan(fco2) ? "NULL" : fco2))"
    )

    end

    ProgressMeter.next!(prog)
  end

  # End the transaction and tidy up
  SQLite.commit(db)
  ProgressMeter.finish!(prog)
  close(infile)

  println("Creating index...")
  SQLite.execute(db,
    "CREATE INDEX cellidx ON socat(lonindex, latindex)"
  )

  close(db)

  return nothing
end

function initdb()::SQLite.DB
  if ispath(SQLITE_FILE)
    rm(SQLITE_FILE)
  end

  db = SQLite.DB(SQLITE_FILE)

  timeseriestable = Tables.Schema(
    (:lonindex, :latindex, :year, :dayofyear, :season, :minuteofday,
     :expocode, :sst, :sss, :fco2),
    Tuple{Int64,Int64,Int64,Int64,Int64,Int64,
      String,Union{Missing,Float64},Union{Missing,Float64},Union{Missing,Float64}}
  )

  SQLite.createtable!(db, "socat", timeseriestable)

  return db
end

function minuteofday(date::DateTime)::Int64
  return hour(date) * 60 + minute(date)
end

function getseason(dayofyear)::Int64
  if dayofyear ≥ 335 || dayofyear < 60
    return DJF
  elseif dayofyear ≥ 60 && dayofyear < 152
    return MAM
  elseif dayofyear ≥ 152 && dayofyear < 244
    return JJA
  else
    return SON
  end

  return result
end

main()

ext: jl


text:
"""
        marcos(decisionMat, weights, fns)

Apply MARCOS (Measurement Alternatives and Ranking according to COmpromise Solution) for a given matrix and weights.

# Arguments:
 - `decisionMat::DataFrame`: n × m matrix of objective values for n alterntives and m criteria
 - `weights::Array{Float64, 1}`: m-vector of weights that sum up to 1.0. If the sum of weights is not 1.0, it is automatically normalized.
 - `fns::Array{Function, 1}`: m-vector of functions to be applied on the columns.

# Description
marcos() applies the MARCOS method to rank n alterntives subject to m criteria which are supposed to be
either maximized or minimized.

# Output
- `::MARCOSResult`: MARCOSResult object that holds multiple outputs including scores, rankings, and best index.

# Examples
```julia-repl
julia> decmat = [8.675 8.433 8.000 7.800 8.025 8.043;
       8.825 8.600 7.420 7.463 7.825 8.229;
       8.325 7.600 8.040 7.700 7.925 7.600;
       8.525 8.667 7.180 7.375 7.750 8.071]
4×6 Array{Float64,2}:
 8.675  8.433  8.0   7.8    8.025  8.043
 8.825  8.6    7.42  7.463  7.825  8.229
 8.325  7.6    8.04  7.7    7.925  7.6
 8.525  8.667  7.18  7.375  7.75   8.071

julia> df = makeDecisionMatrix(decmat)

4×6 DataFrame
 Row │ Crt1     Crt2     Crt3     Crt4     Crt5     Crt6
     │ Float64  Float64  Float64  Float64  Float64  Float64
─────┼──────────────────────────────────────────────────────
   1 │   8.675    8.433     8.0     7.8      8.025    8.043
   2 │   8.825    8.6       7.42    7.463    7.825    8.229
   3 │   8.325    7.6       8.04    7.7      7.925    7.6
   4 │   8.525    8.667     7.18    7.375    7.75     8.071

julia> weights = [0.1901901901901902 , 0.15915915915915918 , 0.19819819819819823, 0.1901901901901902, 0.15115115115115116, 0.11111111111111112];

julia> fns = [maximum, maximum, maximum, maximum, maximum, maximum];

julia> Fns = convert(Array{Function, 1} , fns)

julia> result = marcos(df, weights, Fns);

julia> result.scores
4-element Array{Float64,1}:
 0.6848657890705123
 0.6727670074308345
 0.6625969531206817
 0.6611030275027843

julia> result.bestIndex
1
```

# References
Stević, Z., Pamučar, D., Puška, A., Chatterjee, P., Sustainable supplier selection in healthcare industries using a new MCDM method: Measurement Alternatives and Ranking according to COmpromise Solution (MARCOS), Computers & Industrial Engineering (2019), doi: https://doi.org/10.1016/j.cie.2019.106231

Puška, A., Stojanović, I., Maksimović, A., & Osmanović, N. (2020). Evaluation software of project management used measurement of alternatives and ranking according to compromise solution (MARCOS) method. Operational Research in Engineering Sciences: Theory and Applications, 3(1), 89-102.
"""
function marcos(decisionMat::DataFrame, weights::Array{Float64,1}, fns::Array{Function,1})::MARCOSResult

    # df = convert(Matrix, decisionMat)
    df = Matrix(decisionMat)

    row, col = size(df)

    w = unitize(weights)

    AAI = zeros(Float64, col)
    AI  = zeros(Float64, col)

    temp = [df; AI'; AAI']

    normalizedDecisionMat = similar(temp)

    @inbounds for i in 1:col
        if fns[i] == maximum
            AI[i] = maximum(df[:, i])
            temp[row + 1, i] = AI[i]
            AAI[i] = minimum(df[:, i])
            temp[row + 2, i] = AAI[i]
            normalizedDecisionMat[:, i] = temp[:, i] ./ AI[i]
        elseif fns[i] == minimum
            AI[i] = minimum(df[:, i])
            temp[row + 1, i] = AI[i]
            AAI[i] = maximum(df[:, i])
            temp[row + 2, i] = AAI[i]
            normalizedDecisionMat[:, i] = AI[i] ./ temp[:, i]
        end
    end

    S  = zeros(Float64, col)

    for i in 1:row + 2
        S[i] = w .* normalizedDecisionMat[i, :] |> sum
    end

    KPlus  = S[1:row] ./ S[row + 1]
    KMinus  = S[1:row] ./ S[row + 2]

    fKPlus = KPlus ./ (KPlus .+ KMinus)
    fKMinus = KMinus ./ (KPlus .+ KMinus)

    scores = zeros(Float64, row)
    for i in 1:row
        scores[i] = (KPlus[i] + KMinus[i]) / ((1 + (1 - fKPlus[i]) / fKPlus[i]) + ((1 - fKMinus[i]) / fKMinus[i]))
    end
    rankings = sortperm(scores)

    bestIndex = rankings |> last

    result = MARCOSResult(
        decisionMat,
        w,
        scores,
        rankings,
        bestIndex
    )

    return result
end

max_stars_repo_name: bahadirfyildirim/JMcDM

max_issues_repo_name: bahadirfyildirim/JMcDM

max_forks_repo_name: bahadirfyildirim/JMcDM


text:
#!/usr/bin/julia

using Pkg
#Pkg.add("StatsModels")
using ExperimentalDesign, StatsModels, GLM, DataFrames, Distributions, Random
include("gcc.jl")

function generate_file_distribution()
	prefix = "processed"

	files = readdir(prefix)
	maximum = 0

	for file in files
		splitted = split(file, ".")
		base = splitted[1]
		ext = splitted[2]

		base_split = split(base, "-")
		num_func_str = base_split[1]
		num_funcs = parse(Int64, num_func_str)

		maximum = max(maximum, num_funcs)
	end

	names = fill(String[], maximum)

	for i in 1:maximum
		names[i] = String[]
	end

	for file in files
		splitted = split(file, ".")
		base = splitted[1]
		ext = splitted[2]

		base_split = split(base, "-")
		num_func_str = base_split[1]
		size_str = base_split[2]
		ret_str = base_split[3]

		num_funcs = parse(Int64, num_func_str)

		if num_funcs == 0
			continue
		end

		push!(names[num_funcs], file)
	end

	return names
end

function compute_fill_rate(names)
	len = length(names)
	acc = 0

	for i in 1:len
		if names[i] != []
			acc = acc + 1
		end
	end

	print(acc/len)
end

function y(x)
    return run_gcc(x)
end

# Used for DiscreteNonParametric(fp, prob), but use DiscreteUniform for now.
# fp   = [20 + 14*n for n in 0:14]
# prob = [1/15 for n in 0:14]

#random_design_generator = RandomDesign((
#          Distributions.DiscreteUniform(1, 300),     # File size, in number of functions
#          Distributions.DiscreteUniform(1, 10000),   # Min lto partition
#          Distributions.Bernoulli(0.5)))             # Run in parallel?
#
#println(random_design_generator)
#
#random_design = rand(random_design_generator, 8)
#random_design[!, :response] = y.(eachrow(random_design[:, :]))
#
#println(random_design)

names = generate_file_distribution()
compute_fill_rate(names)

max_stars_repo_name: giulianobelinassi/LTO-Timing-Analysis

max_issues_repo_name: giulianobelinassi/LTO-Timing-Analysis

max_forks_repo_name: giulianobelinassi/LTO-Timing-Analysis


text:
makekernel(f,r::UnitRange{Int}) = Kernel{(r,)}(w->f(Tuple(w)))
makekernel(f,t::NTuple{2,Int}) = Kernel{(first(t):last(t),)}(w->f(Tuple(w)))
makekernel(f,A::NTuple{<:Any,UnitRange{Int}}) = Kernel{A}(w->f(Tuple(w)))

function applyrolling(f,data,window,filler)
    k=makekernel(f,window)
    return map(k,extend(data,StaticKernels.ExtensionConstant(filler)))
end

function applyrolling(k::Kernel,data,filler)
    return map(k,extend(data,StaticKernels.ExtensionConstant(filler)))
end

function applyrolling(k::Kernel,data)
    return map(k,extend(data,StaticKernels.ExtensionNothing()))
end

logdiff(x) = log(last(x)) - log(first(x))

# below also works but needs to be cleaned up
# for w in 1:12
#     for f in (:sum,:mean,:std,:var,:minimum,:maximum,:median,:middle,:logdiff)
#         nm = Symbol("RollingKernel_",w,"_",f)
#         @eval const $nm = Kernel{(-$w:0,)}(x -> $f(Tuple(x)))
#     end
# end


text:
# As of this date (May 2016) nl2sno (finite difference calc of jacobian)
# is still encountering random memory corruption issues.  Use at you own
# risk...
@enum OptoAlgo  NL2SOL NL2SNO

mutable struct Regression{T<:Number}
    sigma::T
    rms::T
    coef_determ::T
    covmx::Matrix{T}
    cormx::Matrix{T}
    std_param::Vector{T}
    t_ratio::Vector{T}
    std_wresid::Vector{T}
    results::Any
    rank::Int
end

function Regression(T, q, n, m)
    Regression(
        zero(T),              # sigma
        zero(T),              # rms
        zero(T),              # coef_determ
        zeros(T, q+n, q+n),   # covmx
        zeros(T, q+n, q+n),   # cormx
        zeros(T, m),
        zeros(T, q+n),
        zeros(T, m),
        0,
        n
    )
end

# We intend T to be either Float64 or Complex128
mutable struct FitContext{T<:Number}
    # Required
    y::Vector{T}         # the m 'measurements'. Can be real or complex
    t::Vector{Float64}   # Independent variable (at m points) and is always real
    w::Matrix{Float64}   # the weight matrix (input initially as vector)
    alpha::Vector{T}     # Initial value of nonlinear variables. Can be real or complex
    n::Int               # number of linear parameters
    ind::Matrix{Int}     # integer matrix maps non zero entries of dphi
    ada::Function        # function to calculate phi
    gada::Function       # function to calculate dphi

    # Optional (only to change defaults)
    n1::Int              # Default to n
    debug::Bool          # turn on debugging
    neglect::Bool        # this turns on the Kaufman approximation to the jacobian
    opto::OptoAlgo       # choices are NL2SOL, NL2SNO (fd-jacobian), or LEVENBERG
    verbose::Bool        # More chatty if true
    mxfcal::Int          # Max number of function call for NL2sol
    mxiter::Int          # Max iterations for NL2sol

    # "Internal" values.  Some are for memory preallocation
    q::Int               # length(alpha)
    m::Int               # length(y)
    c::Vector{T}         # the linear variables
    y_est::Vector{T}     # current estimate of y
    rank::Int            # Rank of the phi matrix
    phi::Matrix{T}       # Current estimate of phi
    dphi::Matrix{T}      # current estimate of derivative of phi
    wresid::Vector{T}    # current weighted residual
    jac1::Matrix{T}      #
    jac2::Matrix{T}      #
    jac::Matrix{T}       #
    iscomplex::Bool      # are y and alpha complex?

    U::Matrix{T}         # SVD decomp of
    s::Vector{T}
    V::Matrix{T}
    alpha_real::Vector{Float64}
    wresid_real::Vector{Float64}
    jac_real::Matrix{Float64}
    skip_stats::Bool          # skip all the stats stuff
end

function FitContext(y::Vector{T}, t, w, alpha::Vector{T}, n, ind, ada, gada) where {T}
    if length(size(y)) != 1 || length(size(t)) != 1 || length(size(w)) != 1 ||
       length(size(alpha)) != 1
        error("y, t, w, and alpha must all be vectors")
    end
    q = length(alpha)
    m = length(y)
    if length(t) != m || length(w) != m
        error("t (independent variable) and w (weight vector) must be same " *
              "size as y (sample vector)")
    end

    a, b = size(ind)
    if a != 2
        error("ind must by a 2 by mm array, where mm is the number of nonzero" *
              " columns of dphi")
    end

    if b > n*q
        error("The max size of the second dimension of ind is $n * $q")
    end
    n1 = n

    ctx = FitContext(
        # required
        y,                       # sample values to be fit
        t,                       # independent variable
        diagm(0 => w),                # weight matrix  (TODO: matlab code used spdiagm(w, 0, m, m))
        alpha,                   # current estimate on non-linear parameters
        n,                       # n
        ind,                     # ind::Matrix{Int}  constant integer matrix
        ada,                     # phi and residual calculation
        gada,                    # dphi and jacobian calculation

        # Defaults
        n1,                      # n1::Int
        false,                   # debug flag.
        false,                   # neglect jac2 in jacobian calculation
        NL2SOL,                  # default to nl2sol
        true,                    # default verbose on
        200,                     # default max func calls for NL2sol
        500,                     # default max iterations for NL2sol

        # "Internals"
        q,                       # number of nonlinear parameters
        m,                       # number of samples
        zeros(T, n),             # c::Vector{T} current estimate of linear coefficients
        zeros(T, m),             # y_est::Vector{T}: current estimate of modeled y
        n,                       # rank::Int initialized to n

        zeros(T, m, n),          # phi::Matrix{T}    Current estimate of phi
        zeros(T, m, b),          # dphi::Matrix{T}   Current estimate of derivative of phi
        zeros(T, m),             # wresd::Vector{T}  Current weighted residual
        zeros(T, m, q),          # jac1::Matrix{T}
        zeros(T, m, q),          # jac2::Matrix{T}
        zeros(T, m, q),          # jac::Matrix{T}
        T <: Complex,            # are y and t complex?

        zeros(T, m, n),          # U, a factor from the svd decomp
        zeros(T, n),             # s, a factor (as vector) from svd decomp
        zeros(T, n, n),          # V, a factor from the svd decomp

        # Need double size if y and alpha are complex
        zeros(Float64, T <: Complex ? 2*n : n),                         # alpha_real
        zeros(Float64, T <: Complex ? 2*m : m),                         # wresid_real
        zeros(Float64, T <: Complex ? 2*m : m, T <: Complex ? 2*size(ind,2) :
              size(ind,2)),  # jac_real
        false
    )
end


text:
using EvolutionaryModelingTools.Scalar

x = scalar(1)
y = scalar(1im)
@test repr(x) == "scalar(1)"
@test repr(y) == "scalar($(1im))"
@test repr("text/plain", x) == "1"
@test repr("text/plain", y) == repr("text/plain", 1im)
@test x[] == 1
@test x[1] == 1
@test x[1, 1] == 1
@test x[CartesianIndex()] == 1
@test x[CartesianIndex(1)] == 1
@test x[CartesianIndex(1, 1)] == 1
@test (+x)::Int == 1
@test (-x)::Int == -1
@test (x + 1)::Int == 2
@test (x - 0x1)::Int == 0
@test (x * 1.0)::Float64 == 1.0
@test (x / 1.0f0)::Float32 == 1.0
@test convert(typeof(y), x) isa typeof(y)
@test x + x == 2
@test x - y == 1 - im
@test x * y == im
@test x / y == -im
@test x \ y == im
@test x ^ y == 1
@test real(x)::Int == 1
@test imag(x)::Int == 0
@test real(y) == 0
@test imag(y) == 1
@test float(x)::Float64 == 1.0
@test float(y)::ComplexF64 == im
@test exp(x) == exp(1)
@test ℯ^y == exp(im)
x[] = 0x2
@test x[] == 2
x[1] += 1
@test x[1] == 3
x[1, 1] += 1
@test x[1, 1] == 4
y[CartesianIndex()] *= im
@test y[] == -1
y[CartesianIndex(1)] *= im
@test y[1] == -im
y[CartesianIndex(1, 1)] *= im
@test y[1, 1] == 1
@test exp(x) == exp(4)
@test ℯ^y == exp(1)
@test promote_type(typeof(x), Int) == Int
@test promote_type(typeof(x), Float64) == Float64
@test promote_type(typeof(x), typeof(y)) == Complex{Int}


text:
GL_COLOR_ATTACHMENT(i) = GLuint(GL_COLOR_ATTACHMENT0 + i)
#TODO talk about Contexts!
"""
Holds the id, format and attachment of an OpenGL RenderBuffer.
RenderBuffers cannot be read by Shaders.
"""
mutable struct RenderBuffer{T}
    id        ::GLuint
    format    ::GLenum
    attachment::GLenum
    context   ::Context
    size      ::Tuple{Int, Int}
    function RenderBuffer{T}(format::GLenum, attachment::GLenum, dimensions) where {T}
        @assert length(dimensions) == 2
        id = glGenRenderbuffers(format, attachment, dimensions)
        obj = new(id, format, attachment, current_context(), dimensions)
        finalizer(free!, obj)
        obj
    end
end

"Creates a `RenderBuffer` with purpose for the depth component `T` of a `FrameBuffer`."
function RenderBuffer(::Type{T}, dimensions) where  {T <: DepthFormat}
    RenderBuffer{T}(gl_internal_format(T), gl_attachment(T), dimensions)
end

free!(rb::RenderBuffer) =
    context_command(rb.context, glDeleteRenderbuffers(1, [rb.id]))

bind(b::RenderBuffer) = glBindRenderbuffer(GL_RENDERBUFFER, b.id)
unbind(b::RenderBuffer) = glBindRenderbuffer(GL_RENDERBUFFER, b.id)

Base.size(b::RenderBuffer) = b.size
Base.eltype(r::RenderBuffer{T}) where {T} = T
width(b::RenderBuffer)  = size(b, 1)
height(b::RenderBuffer) = size(b, 2)
depth(b::RenderBuffer)  = size(b, 3) #possible need to assert

function resize_nocopy!(b::RenderBuffer, dimensions)
    bind(b)
    b.size = dimensions
    glRenderbufferStorage(GL_RENDERBUFFER, b.format, dimensions...)
end

"""
A FrameBuffer holds all the data related to the usual OpenGL FrameBufferObjects.
The `attachments` field gets mapped to the different possible GL_COLOR_ATTACHMENTs, which is bound by GL_MAX_COLOR_ATTACHMENTS,
and to one of either a GL_DEPTH_ATTACHMENT or GL_DEPTH_STENCIL_ATTACHMENT.
"""
mutable struct FrameBuffer{ElementTypes, T}
    id         ::GLuint
    attachments::T
    context    ::Context
    function FrameBuffer(fb_size::NTuple{2, Int}, attachments::Union{Texture, RenderBuffer}...; kwargs...)
        framebuffer = glGenFramebuffers()
        glBindFramebuffer(GL_FRAMEBUFFER, framebuffer)
        max_ca = glGetIntegerv(GL_MAX_COLOR_ATTACHMENTS)

        depth_attachments = Union{Texture, RenderBuffer}[]
        for (i, a) in enumerate(attachments)
            if eltype(a) <: DepthFormat
                push!(depth_attachments, a)
            else
                attach2framebuffer(a, GL_COLOR_ATTACHMENT(i-1))
            end
        end

        if length(depth_attachments) > 1
            error("The amount of DepthFormat types in texture types exceeds the maximum of 1.")
        end

        if length(attachments) > max_ca
            error("The length of texture types exceeds the maximum amount of framebuffer color attachments! Found: $N, allowed: $max_ca")
        end

        !isempty(depth_attachments) && attach2framebuffer(depth_attachments[1])
        #this is done so it's a tuple. Not entirely sure why thats better than an
        #array?
        _attachments = (attachments...,depth_attachments...,)

        @assert glCheckFramebufferStatus(GL_FRAMEBUFFER) == GL_FRAMEBUFFER_COMPLETE "FrameBuffer (id $framebuffer) with attachments $attachment_types failed to be created."
        glBindFramebuffer(GL_FRAMEBUFFER, 0)
        obj = new{Tuple{eltype.(_attachments)...}, typeof(_attachments)}(framebuffer, _attachments, current_context())
        finalizer(free!, obj)
        return obj
    end
    function FrameBuffer(::Val{0})
        obj = new{Val{0}, Vector{Nothing}}(0, [nothing])
        return obj
    end
end

# Constructor that takes care of creating the FBO with the specified types and then filling in the data
# Might be a bit stupid not to just put this into the main constructor
function FrameBuffer(fb_size::Tuple{<: Integer, <: Integer}, texture_types, texture_data::Vector{<:Matrix}; kwargs...)
	fbo = FrameBuffer(fb_size, texture_types; kwargs...)
	for (data, attachment) in zip(texture_data, fbo.attachments)
		xrange = 1:size(data)[1]
		yrange = 1:size(data)[2]
		gpu_setindex!(attachment, data, xrange, yrange)
	end
	return fbo
end

context_framebuffer() = FrameBuffer(Val(0))

function attach2framebuffer(t::Texture{T, 2}, attachment) where T
    glFramebufferTexture2D(GL_FRAMEBUFFER, attachment, GL_TEXTURE_2D, t.id, 0)
end
function attach2framebuffer(t::Texture{T, 2}) where {T <: DepthFormat}
    glFramebufferTexture2D(GL_FRAMEBUFFER, gl_attachment(T), GL_TEXTURE_2D, t.id, 0)
end
function attach2framebuffer(x::RenderBuffer)
    glFramebufferRenderbuffer(GL_FRAMEBUFFER, x.attachment, GL_RENDERBUFFER, x.id)
end

function free!(fb::FrameBuffer)
    for attachment in fb.attachments
        free!(attachment)
    end
    context_command(fb.context, () -> glDeleteFramebuffers(1, [fb.id]))
    return
end

Base.size(fb::FrameBuffer) = size(fb.attachments[1])
width(fb::FrameBuffer)     = size(fb, 1)
height(fb::FrameBuffer)    = size(fb, 2)
depth(fb::FrameBuffer)     = size(fb, 3)

eltypes(fb::FrameBuffer{elt}) where elt                = elt.parameters
attachtypes(fb::FrameBuffer{elt, int}) where {elt,int} = int.parameters

Base.resize!(fb::FrameBuffer{Val{0}, Vector{Nothing}}, dimensions) = glViewport(0, 0, dimensions...)

function Base.resize!(fb::FrameBuffer, dimensions)
    ws = dimensions[1], dimensions[2]
    if ws != size(fb) && all(x -> x > 0, dimensions)
        dimensions = tuple(dimensions...)
        for attachment in fb.attachments
            resize_nocopy!(attachment, dimensions)
        end
    end
    nothing
end

bind(fb::FrameBuffer, target=GL_FRAMEBUFFER) = glBindFramebuffer(target, fb.id)
unbind(fb::FrameBuffer) = glBindFramebuffer(GL_FRAMEBUFFER, 0)

#I think a lot could be optimized when knowing for sure whether an FBO has depth
#or not, and where it is located. Could even just be a different type.
function draw(fb::FrameBuffer)
    ntexts = length(color_attachments(fb))
    glDrawBuffers(GLuint(ntexts), GL_COLOR_ATTACHMENT.(0:ntexts-1))
end
draw(fb::FrameBuffer, i::Int) = glDrawBuffer(GL_COLOR_ATTACHMENT.(i-1))
function draw(fb::FrameBuffer, i::AbstractUnitRange)
    ntexts = length(color_attachments(fb)[i])
    glDrawBuffers(GLuint(ntexts), GL_COLOR_ATTACHMENT.(i.-1))
end
draw(fb::FrameBuffer{Val{0}, Nothing}) = nothing

#All this is not very focussed on performance yet
function clear!(fb::FrameBuffer, color)
    glClearColor(GLfloat(color[1]), GLfloat(color[2]), GLfloat(color[3]), GLfloat(color[4]))
    draw(fb)
    glClear(GL_COLOR_BUFFER_BIT)
    fm = depthformat(fb)
    if fm <: DepthStencil
        glClear(GL_STENCIL_BUFFER_BIT)
        glClear(GL_DEPTH_BUFFER_BIT)
    elseif fm <: Depth
        glClear(GL_DEPTH_BUFFER_BIT)
    end
end
clear!(fb::FrameBuffer) = clear!(fb, (0.0, 0.0, 0.0, 0.0))
clear!(fb::FrameBuffer{Val{0}, Vector{Nothing}}) = nothing
color_attachments(fb::FrameBuffer) =
    fb.attachments[findall(x -> !(x <: DepthFormat), eltypes(fb))]

color_attachment(fb::FrameBuffer, i) =
    fb.attachments[findall(x -> !(x <: DepthFormat), eltypes(fb))[i]]

depth_attachment(fb::FrameBuffer) =
    fb.attachments[findfirst(x -> x <: DepthFormat, eltypes(fb))]

function depthformat(fb::FrameBuffer)
    id = findfirst(x -> x <: DepthFormat, eltypes(fb))
    if id != 0
        return eltypes(fb)[id]
    else
        return Void
    end
end

gpu_data(fb::FrameBuffer, i) = gpu_data(fb.attachments[i])


text:
mutable struct AlphaEdge
    vec::Vector{Float64}
    id::UInt64
end
Base.hash(a::AlphaEdge, h::UInt64=zero(UInt64)) = hash(a.vec, hash(a.id, h))

mutable struct MCVINode
    id::UInt64
    act::Any
    states::Any
    alpha_edges::Vector{AlphaEdge}
    MCVINode() = new()
    MCVINode(id, act, states, aes) = new(id, act, states, aes)
end

hasnext(n::MCVINode) = !isempty(n.alpha_edges)

copy(n::MCVINode) = MCVINode(n.id, n.act, n.states, n.alpha_edges)

"""
Returns the next node in the policygraph
"""
mutable struct MCVIUpdater{P<:POMDP} <: POMDPs.Updater
    problem::P
    root
    root_belief
    nodes::Dict{UInt64, MCVINode}
    nodes_queue::Vector{MCVINode}
end

MCVIUpdater(problem) = MCVIUpdater(problem, nothing, nothing, Dict{UInt64, MCVINode}(), Vector{MCVINode}())

Base.length(ps::MCVIUpdater) = length(ps.nodes)

hasroot(ps::MCVIUpdater) = ps.root != nothing

hasnode(ps::MCVIUpdater, n::MCVINode) = haskey(ps.nodes, n.id)

create_belief(ps::MCVIUpdater) = MCVINode()

function initialize_belief(up::MCVIUpdater, b::Any)
    if up.root_belief != b
        if @implemented initialstate_distribution(::typeof(up.problem))
            is = initialstate_distribution(up.problem)
        else
            is = "<initialstate_distribution(::$(typeof(up.problem))) not implemented>"
        end
        @warn("""
             The belief used to start MCVI policy execution was (potentially) different from the initial belief used in the MCVI solution.

             updater root belief: $(up.root_belief)
             b0 for policy execution: $b
             initialstate_distribution(pomdp): $is
             """)
    end
    return up.root
end

function create_node(ps::MCVIUpdater, a::A, states::Any, alpha_edges::Vector{AlphaEdge}) where {A}
    if states == nothing
        st_hash = Base.hash(states)
    else
        st_hash = hash(states) # TODO: recursive hash of array?
    end
    id = hash(a, hash(alpha_edges, st_hash))

    return MCVINode(id, a, states, alpha_edges)
end

mutable struct MCVIPolicy <: POMDPs.Policy
    problem::POMDPs.POMDP
    updater::Union{Nothing, POMDPs.Updater}
    MCVIPolicy() = new()
    MCVIPolicy(p) = new(p, nothing)
    MCVIPolicy(p, up) = new(p, up)
    MCVIPolicy(p, up, root) = new(p, up)
end

function init_node(ps::MCVIUpdater, problem::POMDPs.POMDP)
    return create_node(ps, init_lower_action(problem), nothing, Vector{AlphaEdge}())
end

function init_nodes(policy::MCVIPolicy, ps::MCVIUpdater)
    ns = Vector{MCVINode}()
    for a in actions(policy.problem)
        if isterminal(policy.problem, a) # FIXME: Can use init_lower_action instead?
            push!(ns, create_node(ps, a, nothing, Vector{AlphaEdge}()))
        end
    end
    return ns
end

function addnode!(ps::MCVIUpdater, n::MCVINode)
    ps.nodes[n.id] = n
    push!(ps.nodes_queue, n)
end

updater(policy::MCVIPolicy) = policy.updater

"""
Initialize and return the MCVIUpdater
"""
function initialize_updater!(policy::MCVIPolicy)
    ps = MCVIUpdater(policy.problem)
    ns = init_nodes(policy, ps)
    for n in ns
        addnode!(ps, n)
        ps.root = n
    end
    policy.updater = ps
end

"""
Move to the next policy state given observation
"""
function update(ps::MCVIUpdater, n::MCVINode, act::A, obs::O, np::MCVINode=create_belief(ps)) where {A,O}
    @assert hasnext(n) "No next policy state exists"

    local nid::UInt64
    if length(n.alpha_edges) == 1
        nid = n.alpha_edges[1].id
    else
        @assert n.states != nothing "No states in the current node"
        obs_wt = zeros(length(n.states))
        for (i,s) in enumerate(n.states)
            obs_wt[i] = obs_weight(ps.problem, act, s, obs)
        end
        maxv = -Inf
        for ae in n.alpha_edges
            v = dot(obs_wt, ae.vec)
            if v > maxv
                maxv = v
                nid = ae.id
            end
        end
    end
    @assert ps.nodes[nid] != nothing
    np = ps.nodes[nid]
    return np
end

"""
Return all next policy states for the given node
"""
function next_nodes(ps::MCVIUpdater, n::MCVINode)
    nodes = Vector{MCVINode}(length(n.alpha_edges))
    for i in 1:length(nodes)
        nid = n.alpha_edges[i].id
        nodes[i] = ps.nodes[nid]
    end
    return nodes
end

function to_dict(n::MCVINode)
    return Dict("id"          => n.id,
                "act"         => n.act,
                "states"      => n.states,
                "alpha_edges" => n.alpha_edges)
end

"""
Dump policy graph for inspection
"""
function dump_json(policy::MCVIPolicy, filename::AbstractString)
    root = policy.updater.root
    nodes = policy.updater.nodes
    open(filename, "w") do f
        write(f, JSON.json(Dict(root.id => to_dict(root))))
        for k in keys(nodes)
            write(f, JSON.json(Dict(k => to_dict(nodes[k]))))
        end
    end
end


text:
search("Julia","uli") #=> 2:4
typeof(1:1000) #=> UnitRange{Int64}

a = split("A,B,C,D",",")
typeof(a) #=> Array{SubString{ASCIIString},1}
show(a) #=> SubString{ASCIIString}["A","B","C","D"]

arr = [100, 25, 37]
arr[1] #=> 100
arr[end] #=> 37
# arr[6] #=> BoundsError
println(eltype(arr)) #=> Int64
println(length(arr)) #=> 3
println(ndims(arr)) #=> 1
println(size(arr,1)) #=> 3
println(size(arr)) #=> (3,)

arr2 = Array(Int, 5)
show arr2 #=> [0,0,0,0,0]
sizehint(arr2, 10^5)

a = Float64[] #=> 0-element Array{Float64,1}
b = [1:7] #=> 7-element Array{Int64,1}: [1,2,3,4,5,6,7]
join(b, ", ") #=> "1, 2, 3, 4, 5, 6, 7"
matrix = [1 2; 3 4] # => 2x2 Int64 Array: [1 2; 3 4]

b[1:3] # => [1, 2, 3]
b[2:end] # => [2, 3, 4, 5, 6, 7]  , b[2:] is deprecated
c = [100,200,300]
append!(b,c) # Now b is [1, 2, 3, 4, 5, 6, 7, 100, 200, 300]
pop!(b) #=> 300, b is now [1, 2, 3, 4, 5, 6, 7, 100, 200]
push!(b, 42) # b is now [1, 2, 3, 4, 5, 6, 7, 100, 200, 42]
shift!(b) #=> 1, b is now [2, 3, 4, 5, 6, 7, 100, 200, 42]
unshift!(b, 42) # b is now [42, 2, 3, 4, 5, 6, 7, 100, 200, 42]
splice!(b,8) #=> 100, b is now [42, 2, 3, 4, 5, 6, 7, 200, 42]
in(42, b) #=> true
in(43, b) #=> false
sort(b) #=> [2,3,4,5,6,7,42,42,200], b is not changed
sort!(b) #=> b is now changed to [2,3,4,5,6,7,42,42,200]

arr = [1, 5, 3]
# looping
for e in arr
    print("$e ")
end # prints 1 5 3
arr + 2 #=> [3, 7, 5]
arr * 2 #=> [2, 10, 6]

ma1 = rand(3, 5)
ndims(ma1) #=> 2
size(ma1) #=> (3, 5)
size(ma1,1) #=> 3
size(ma1,2) #=> 5
length(ma1) #=> 15

idm = eye(3)
zeros(5)
ones(4)
ones(3, 2)
eqa = linspace(0, 10, 5)
show(eqa) #=> [0.0,2.5,5.0,7.5,10.0]
fill!(arr, 42) #=> [42,42,42]

A = [1,2,4,6]
A0 = A
show(A0)
A[4] = 0
show(A) #=> [1,2,4,0]
show(A0) #=> [1,2,4,0]

max_stars_repo_path: Module 1/Chapter02/Old/arrays.jl

max_issues_repo_path: Module 1/Chapter02/Old/arrays.jl

max_forks_repo_path: Module 1/Chapter02/Old/arrays.jl


text:
#error analysis of fault prediction network
cd(@__DIR__)
using Pkg
Pkg.activate(".")

using Flux
using Plots
using BSON: @load
using Random
using Statistics

#load network to test
@load "Models/Fault_Model.bson" Fault_Model

#function for parsing data
function get_data(filename)
    f1 = open(filename)
    tempdata = readlines(f1)
    tempdata[1] = tempdata[1][4:end]  #cutting of beggining bc some weird chars
    close(f1)
    data = []
    regvec = zeros(length(tempdata),
                   length(split(tempdata[1], ",")))
    #parsing datavector to float array... not regularized
    for i in 1:length(tempdata)
        regvec[i,:] =  parse.(Float64,split(tempdata[i], ","))
    end
    return regvec
end

#normalizing data
function norm_data(data, scalingmatrix=nothing)
    #getting regularization values and normalizing data
    regvec = data
    normvec = zeros(length(regvec[:,1]),
                    length(regvec[1,:]))
    #checking if optional scalingmatrix was passed, if not will calc mean and std
    if scalingmatrix == nothing
        scalingmatrix = zeros(length(regvec[1,:]),2)
        for i in 1:length(regvec[1,:])
            scalingmatrix[i,1] = Statistics.mean(regvec[:,i]) #first index is mean
            scalingmatrix[i,2] = Statistics.std(regvec[:,i])  #second index is std
        end
        for i in 1:length(regvec[1,:])
            #normalizing regvec
            normvec[:,i] = (regvec[:,i] .- scalingmatrix[i,1]) ./ scalingmatrix[i,2]
        end
    else
        #for case of normalizing test data by training std and mean
        for i in 1:length(regvec[1,:])
            #normalizing regvec
            normvec[:,i] = (regvec[:,i] .- scalingmatrix[i,1]) ./ scalingmatrix[i,2]
        end
    end
    return normvec, scalingmatrix
end

#function to convert vector of vector to matrix
function vecvec_to_matrix(vecvec)
    dim1 = length(vecvec)
    dim2 = length(vecvec[1])
    my_array = zeros(dim1, dim2)
    for i in 1:dim1
        for j in 1:dim2
            my_array[i,j] = vecvec[i][j]
        end
    end
    return my_array
end

#getting training data
fname = "Dataset/FINAL_DATA.txt"
data_tmp = get_data(fname)
#seperating test and train data
train_data_tmp = []; test_data_tmp = [];
for i in 1:length(data_tmp[:,1])
    if data_tmp[i,2] in [0.3,0.38,0.46,0.53,0.61,0.65] #training case
        push!(train_data_tmp, data_tmp[i,:])
    else #testing case
        push!(test_data_tmp, data_tmp[i,:])
    end
end
train_data = vecvec_to_matrix(train_data_tmp)
test_data = vecvec_to_matrix(test_data_tmp)
#normalizing data
train_norm, scalingmatrix = norm_data(train_data)
test_norm, scalingmatrix = norm_data(test_data,scalingmatrix)
data_norm, scalingmatrix = norm_data(data_tmp,scalingmatrix)
#converting data to format which works with network
data = []
for i in 1:length(train_norm[:,1])
    #adjust x and y for fault detection here
    push!(data,(train_norm[i,[1,2,3]],train_norm[i,4]))
end

#Normalized Root Mean Squared Error
function RMSE(data_set,model)
    #calculate root mean squared error
    N = length(data_set[:,1])
    print(N,'\n')
    MSE_sum = 0.0
    for i in 1:N
        MSE_sum += 1/N * (data_set[i,3] - model(data_set[i,[1,2,4]])[1])^2.0
    end
    return sqrt(MSE_sum)  # no norm needed bc data is already normalized
end

print("Error Anlysis")
print("\nOverall root mean squared error: ")
print(RMSE(data_norm,Fault_Model))
print("\nTest root mean squared error: ")
print(RMSE(test_norm,Fault_Model))
print("\nTrain root mean squared error: ")
print(RMSE(train_norm,Fault_Model))
#TODO: represent interval of prediction (normalied goes -2 to 0.5 ish)

#plotting
alt_mean = scalingmatrix[1,1]
alt_std = scalingmatrix[1,2]
mach_SCAT = [];pressure_SCAT = [];fault_SCAT = [];model_SCAT = [];
specified_alt = (12000 - alt_mean)/alt_std
for i in 1:length(test_norm[:,1])
    if test_norm[i,1] == specified_alt
        push!(mach_SCAT,test_norm[i,2])
        push!(fault_SCAT,test_norm[i,3])
        push!(pressure_SCAT,test_norm[i,4])
        push!(model_SCAT,Fault_Model(test_norm[i,[1,2,4]])[1])
    end
end
#creating 3d scatter for showing network results
scatter(pressure_SCAT,mach_SCAT,fault_SCAT,label="Actual",
        title="Predictions with Testing Data [All]",
        xlabel="Pressure",
        ylabel="Mach",
        zlabel="Fault Parameter")
scatter!(pressure_SCAT,mach_SCAT,model_SCAT,label="Predicted")

max_stars_repo_name: NoahLopezGit/PitotTubeNeuralNetwork

max_issues_repo_name: NoahLopezGit/PitotTubeNeuralNetwork

max_forks_repo_name: NoahLopezGit/PitotTubeNeuralNetwork


text:
using Test, GeoDatasets

include("land_sea_mask.jl")

@testset "gshhg" begin
    res = 'c'
    @test length(GeoDatasets.gshhg(res,[1,5])) > 300
    @test length(GeoDatasets.WDBII_river(res,[1,2])) > 1000
    @test length(GeoDatasets.WDBII_border(res,[1])) >  100

    # resolution 'x' does not exists
    @test_throws ErrorException GeoDatasets.gshhg('x',[1,5])
end

max_stars_repo_name: JuliaGeo/GeoDatasets

max_issues_repo_name: JuliaGeo/GeoDatasets

max_forks_repo_name: JuliaGeo/GeoDatasets


text:
"""
    eom_var1(x, p, n) -> Function

Equations of motion for a vector autoregressive system where X₁ → X₂ → X₃.
"""
function eom_var1(x, p, n)
    σ₁, σ₂, σ₃ = p[1], p[2], p[3]
    x₁, x₂, x₃ = x[1], x[2], x[3]
    θ = rand(Normal(0, σ₁))
    η = rand(Normal(0, σ₂))
    ϵ = rand(Normal(0, σ₃))

    dx₁ = θ
    dx₂ = x₁ * η
    dx₃ = 0.5*x₃ * x₂ + ϵ
    return SVector{3}(dx₁, dx₂, dx₃)
end

function var1(uᵢ, σ₁, σ₂, σ₃)
    p = [σ₁, σ₂, σ₃]
    DiscreteDynamicalSystem(eom_var1, uᵢ, p)
end

"""
    var1(x, p, n) -> DiscreteDynamicalSystem

Initialise a discrete vector autoregressive system where X₁ → X₂ → X₃.
"""
var1(;uᵢ = rand(3), σ₁ = 1.0, σ₂ = 0.2, σ₃ = 0.3) = var1(uᵢ, σ₁, σ₂, σ₃)


text:
using Plots
# Model
include_model("particle")
model = Particle{Discrete, FixedTime}(n, m, d,
				 mass, 0.0, μ,
				 nq,
				 nu,
				 nc,
				 nf,
				 nb,
				 ns,
				 idx_u,
				 idx_λ,
				 idx_b,
				 idx_ψ,
				 idx_η,
				 idx_s)
# Horizon
T = 101

# Time step
tf = 1.0
h = tf / (T-1)

# Path
px = range(-h, stop = tf, length = T+1)
py = 0.05 * sin.(2 * π * px)

plot(px, py, aspect_ratio = :equal)
q_ref = [[px[t]; py[t]; 1.0] for t = 1:length(px)]
x_ref = configuration_to_state(q_ref)

# Bounds
_uu = Inf * ones(model.m)
_ul = zeros(model.m)
_ul[model.idx_u] .= -Inf

ul, uu = control_bounds(model, T, _ul, _uu)

# Initial and final states
q1 = q_ref[1]#[0.0, 0.0, 1.0]
# v1 = [3.0, 5.0, 0.0]
# v2 = v1 - C_func(model, q1, v1) * h
# q2 = q1 + 0.5 * h * (v1 + v2)
#
# x1 = [q1; q2]

xl, xu = state_bounds(model, T, x1 = [q_ref[1]; q_ref[2]], xT = [q_ref[end-1]; q_ref[end]])

# Objective
obj_tracking = quadratic_tracking_objective(
    [Diagonal(1000.0 * ones(model.n)) for t = 1:T],
    [Diagonal([1.0e-3, 1.0e-3, 1.0e-3,
		ones(model.nc)..., ones(model.nb)...,
		zeros(model.m - model.nu - model.nc - model.nb)...])
		for t = 1:T-1],
   	x_ref,
    [zeros(model.m) for t = 1:T])
obj_penalty = PenaltyObjective(1.0e5, model.m)

obj = MultiObjective([obj_tracking, obj_penalty])

# Constraints
include_constraints(["contact", "loop"])
con_contact = contact_constraints(model, T)
# con_loop = loop_constraints(model, collect([(2:model.nq)..., (nq .+ (2:model.nq))...]), 1, T)
con = multiple_constraints([con_contact])#, con_loop])

# Problem
prob = trajectory_optimization_problem(model,
               obj,
               T,
               h = h,
               xl = xl,
               xu = xu,
               ul = ul,
               uu = uu,
               con = con)

# Trajectory initialization
x0 = deepcopy(x_ref) #linear_interpolation(x1, x1, T) # linear interpolation on state
u0 = [0.001 * rand(model.m) for t = 1:T-1] # random controls

# Pack trajectories into vector
z0 = pack(x0, u0, prob)

#NOTE: may need to run examples multiple times to get good trajectories
# Solve nominal problem
@time z̄, info = solve(prob, copy(z0), tol = 1.0e-5, c_tol = 1.0e-5)

check_slack(z̄, prob)
x̄, ū = unpack(z̄, prob)
q̄ = state_to_configuration(x̄)

q = state_to_configuration(x̄)
u = [u[model.idx_u] for u in ū]
γ = [u[model.idx_λ] for u in ū]
b = [u[model.idx_b] for u in ū]
h̄ = h
@save joinpath(@__DIR__, "particle_sinusoidal.jld2") z̄ x̄ ū h̄ q u γ b

include(joinpath(pwd(), "models/visualize.jl"))
vis = Visualizer()
render(vis)
visualize!(vis, model,
    q̄,
    Δt = h)

plot(hcat(q_ref...)', color = :black, width = 2.0)
plot!(hcat(q̄...)', color = :red, width = 1.0)

vT = (q̄[end] - q̄[end-1]) / h
v1 = (q̄[2] - q̄[1]) / h
@show norm(v1 - vT)
# open(vis)
#
# vT = (q_ref[end] - q_ref[end-1]) / h
# v1 = (q_ref[2] - q_ref[1]) / h
# @show norm(v1 - vT)

max_issues_repo_name: thowell/DirectMotionPlanning


text:
#default implementation
@impl begin
    struct ScoreGetLogScore end
    function get_log_score(sf::Score{I}, i::I)::AbstractFloat where {I}
        return log(get_score(sf, i))
    end
end


text:
using Documenter, MolecularConformation

makedocs(sitename= "MolecularConformation.jl",
				   pages = ["Home"=> "index.md",
							"Quick Start" => "quick.md",
							"API" => "api.md",
							"Contributing"=> "contributing.md",
							"Development Status"=> "devstatus.md"],
				   format = Documenter.HTML(prettyurls = false)

		)

deploydocs(
		   	repo = "github.com/evcastelani/MolecularConformation.jl.git"
				)


text:
using GLVisualize, GeometryTypes, GLAbstraction, Colors, FileIO
if !isdefined(:runtests)
	window = glscreen()
end

# a few helper functions to generate images
typealias NColor{N, T} Colorant{T, N}
fillcolor{T <: NColor{4}}(::Type{T}) = T(0,1,0,1)
fillcolor{T <: NColor{3}}(::Type{T}) = T(0,1,0)

# create different images with different color types (not an exhaustive list of supported types)
arrays = map((RGBA{U8}, RGBA{Float32}, RGB{U8}, RGB{Float32}, BGRA{U8}, BGR{Float32})) do C
     C[fillcolor(C) for x=1:45,y=1:45]
 end
# load a few images from the asset folder with FileIO.load (that's what loadasset calls)
loaded_imgs = map(x->loadasset("test_images", x), readdir(assetpath("test_images")))

# combine them all into one array and add an animated gif and a few other images
x = Any[
    arrays..., loaded_imgs...,
    loadasset("kittens-look.gif"),
    loadasset("mario", "stand", "right.png"),
    loadasset("mario", "jump", "left.gif"),
]

# visualize all images and convert the array to be a vector of element type context
# This shouldn't be necessary, but it seems map is not able to infer the type alone
images = convert(Vector{Context}, map(visualize, x))
# make it a grid
images = reshape(images, (4,4))
# GLVisualize offers a few helpers to visualize arrays of render objects
# spaced out as the underlying array. So this will create a grid whereas every
# item is 128x128x128 pixels big
img_vis = visualize(images, scale=Vec3f0(128))
_view(img_vis, window)

if !isdefined(:runtests)
	renderloop(window)
end


text:
@testset "ch08_fastslam02345" begin
    dt = 0.1
    # environment
    xlim = [-5.0, 5.0]
    ylim = [-5.0, 5.0]
    # id of landmark must start from 0 with 1 step
    landmarks =
        [Landmark([2.0, -3.0], 0), Landmark([3.0, 3.0], 1), Landmark([-4.0, 2.0], 2)]
    envmap = Map()
    push!(envmap, landmarks)
    world = World(xlim, ylim)
    push!(world, envmap)
    # robot side
    initial_pose = [0.0, 0.0, 0.0]
    estimator = FastSlam1(initial_pose, 100, length(landmarks))
    circling_agent = EstimatorAgent(0.2, 10.0 * pi / 180, dt, estimator)
    robot = RealRobot(initial_pose, circling_agent, RealCamera(landmarks); color = "red")
    push!(world, robot)
    anim = @animate for i = 1:400
        t = dt * i
        annota = "t = $(round(t, sigdigits=3))[s]"
        z = observations(
            robot.sensor_,
            robot.pose_;
            noise = true,
            bias = true,
            phantom = true,
        )
        p = draw(world, annota)
        v, ω = decision(circling_agent, z, envmap; resample = true)
        state_transition(
            robot,
            v,
            ω,
            dt;
            move_noise = true,
            vel_bias_noise = true,
            kidnap = false,
        )
    end
    if GUI
        gif(anim, "ch08_fastslam02345.gif", fps = 20)
    end
end

@testset "ch08_fastslam067" begin
    dt = 0.1
    # environment
    xlim = [-5.0, 5.0]
    ylim = [-5.0, 5.0]
    # id of landmark must start from 0 with 1 step
    landmarks =
        [Landmark([2.0, -3.0], 0), Landmark([3.0, 3.0], 1), Landmark([-4.0, 2.0], 2)]
    envmap = Map()
    push!(envmap, landmarks)
    world = World(xlim, ylim)
    push!(world, envmap)
    # robot side
    initial_pose = [0.0, 0.0, 0.0]
    estimator = FastSlam2(initial_pose, 100, length(landmarks))
    circling_agent = EstimatorAgent(0.2, 10.0 * pi / 180, dt, estimator)
    robot = RealRobot(initial_pose, circling_agent, RealCamera(landmarks); color = "red")
    push!(world, robot)
    anim = @animate for i = 1:400
        t = dt * i
        annota = "t = $(round(t, sigdigits=3))[s]"
        z = observations(
            robot.sensor_,
            robot.pose_;
            noise = true,
            bias = true,
            phantom = true,
        )
        p = draw(world, annota)
        v, ω = decision(circling_agent, z, envmap; resample = true)
        state_transition(
            robot,
            v,
            ω,
            dt;
            move_noise = true,
            vel_bias_noise = true,
            kidnap = false,
        )
    end
    if GUI
        gif(anim, "ch08_fastslam067.gif", fps = 20)
    end
end

max_stars_repo_name: soblin/JuliaProbo

max_issues_repo_name: soblin/JuliaProbo

max_forks_repo_name: soblin/JuliaProbo


text:
# kuramoto_7-25-19.jl
#
# Continuing with our research. Trying to mix some cool new literate
# programming techniques in today, too.

#------Cell boundary: Use Alt+⬆ and Alt+⬇ to skip quickly in Atom.-----#

using Plots, DifferentialEquations, BenchmarkTools

#------Cell boundary: Use Alt+⬆ and Alt+⬇ to skip quickly in Atom.-----#

# It takes me a while to get going sometimes, because I use a lot of the
# first part of my day for things like reviewing flashcards of common
# programs I like to use. Once I have all that in place, though, watch
# out -- I'll be a serious force to be reckoned with.

#------Cell boundary: Use Alt+⬆ and Alt+⬇ to skip quickly in Atom.-----#

# Today, I"m going to take a bit and actually read through the Julia
# documentation. I already have enough programming knowledge to do a
# good job working in it, but it's probably not a bad idea to give it
# some of its own time of day as well. These things can be... finicky.

#------Cell boundary: Use Alt+⬆ and Alt+⬇ to skip quickly in Atom.-----#

# Oh, interesting. UTF-8 is allowed in the names of things.

const ϵ = 0.000002
const δ = 0.000001

function f()
    println(ϵ)
    println(δ)
end

f()

#------Cell boundary: Use Alt+⬆ and Alt+⬇ to skip quickly in Atom.-----#

# I did basically no research today, but I *did* learn a lot about the Julia language. I'll pick up again tomorrow starting from
#
# https://docs.julialang.org/en/v1/manual/strings/#Unicode-and-UTF-8-1
#
# 🙂

#------Cell boundary: Use Alt+⬆ and Alt+⬇ to skip quickly in Atom.-----#


text:
struct PolynomialStep
    η :: Real
    c :: Real
end
function get(step::PolynomialStep, k::Real)
    step.c * (k + 1.) ^ (-step.η)
end

struct AdaptiveState
    swap_target_ar :: Real
    scale          :: Base.RefValue{<:Real}
    step           :: PolynomialStep
end
function AdaptiveState(swap_target::Real, scale::Real, step::PolynomialStep)
    AdaptiveState(swap_target, Ref(log(scale)), step)
end

function init_adaptation(
    Δ::Vector{<:Real},
    swap_target::Real,
    scale::Real,
    γ::Real
)
    Nt = length(Δ)
    step = PolynomialStep(γ, Nt - 1)
    Ρ = [AdaptiveState(swap_target, scale, step) for _ in 1:(Nt - 1)]
    return Ρ
end

function rhos_to_ladder(Ρ, Δ)
    β′ = Δ[1]
    for i in 1:length(Ρ)
        β′ += exp(Ρ[i].scale[])
        Δ[i + 1] = Δ[1] / β′
    end
    return Δ
end

function adapt_rho(ρ::AdaptiveState, swap_ar, n)
    swap_diff = swap_ar - ρ.swap_target_ar
    γ = get(ρ.step, n)
    return γ * swap_diff
end

function adapt_ladder(Ρ, Δ, k, swap_ar, n)
    Ρ[k].scale[] += adapt_rho(Ρ[k], swap_ar, n)
    return Ρ, rhos_to_ladder(Ρ, Δ)
end


text:
using JSON

function dump_ir(script, ir_funcs = (code_llvm,))
    @eval Main include($script)
    results = Dict()
    for f in ir_funcs
        results[Symbol(f)] = sub = Dict()
        for (key, args) in Main.IR_DUMP
            sub[key] = sprint(f, args...)
        end
    end
    return results
end

function _json_ir(script, ir_funcs)
    JSON.print(STDOUT, dump_ir(script, ir_funcs))
end

function assert_git_is_clean()
    status = readstring(`git status --short`)
    if strip(status) == ""
        return
    end
    error("""Git repository is not clean:
          $status""")
end

function _subproc_dump_ir(package_name, rev, script)
    ir_funcs = (code_llvm,)
    juliacmd = `julia -O3`
    return cd(Pkg.dir(package_name)) do
        assert_git_is_clean()
        current = strip(readstring(`git rev-parse $rev`))
        try
            run(`git checkout $rev`)
            code = """
            using PkgBenchUtils: _json_ir
            _json_ir($(repr(script)), $ir_funcs)
            """
            return open(JSON.parse, `$juliacmd -e $code`)
        finally
            run(`git checkout $current`)  # rollback
        end
    end
end

struct IRComparison
    package_name
    script
    target
    baseline
    target_results
    baseline_results
end

function compare_ir(script;
                    package_name = guess_package(),
                    target = "HEAD",
                    baseline = "HEAD^")
    return IRComparison(
        package_name,
        script,
        target,
        baseline,
        _subproc_dump_ir(package_name, target, script),
        _subproc_dump_ir(package_name, baseline, script),
    )
end


text:
using CUDA

export NodePairIds,
       BitCircuitPair,
       SumNodePairIds,
       ProdNodePairIds,
       BitCircuitPair,
       ParamBitCircuitPair

const NodeId = Int32

"Integer identifier for a circuit pair node"
const NodePairId = Int32

"The BitCircuitPair ids associated with a node"
abstract type NodePairIds end

mutable struct SumNodePairIds <: NodePairIds
    layer_id::NodePairId
    node_id::NodePairId
end

mutable struct ProdNodePairIds <: NodePairIds
    layer_id::NodePairId
    # Assuming both Product nodes have two children
    left_left_id::NodePairId
    right_right_id::NodePairId

    ProdNodePairIds(ll::SumNodePairIds, rr::SumNodePairIds) = begin
        l = max(ll.layer_id, rr.layer_id)
        new(l, ll.node_id, rr.node_id)
    end
end

const NODES_LENGTH = 6
const ELEMENTS_LENGTH = 3

"""
A bit circuit pair is a low-level representation of pairs of nodes traversed for a pair of probabilistic circuit logical circuit structure.
For example, this is used for taking expectation of a Logistic/Regression circuit w.r.t. a probabilistic circuit.

The wiring of the circuit is captured by two matrices: nodes and elements.
  * Nodes are either leafs (input) or pair of Sum nodes
  * Elements are Pair of Product nodes
  * In addition, there is a vector of layers, where each layer is a list of node ids.
    Layer 1 is the leaf/input layer. Layer end is the circuit root.
  * And there is a vector of parents, pointing to element id parents of decision nodes.

Nodes are represented as a 6xN matrix where
  * nodes[1,:] is the first child pair id
  * nodes[2,:] is the last child pair id
  * nodes[3,:] is the first parent index
  * nodes[4,:] is the last parent index
  * nodes[5,:] is the id of corresponding node from the PC (first circuit)
  * nodes[6,:] is the id of corresponding node from the LC (first circuit)

Elements are represented by a 3xE matrix, where
  * elements[1,:] is the Product pair node id,
  * elements[2,:] is the (left,left) child node id
  * elements[3,:] is the (right right) child node id

  Elements belonging to node pair `i` are `elements[:, nodes[1,i]:nodes[2,i]]`

Parents are represented as a one dimentional array
    Parents belonging to node pair `i` are `parents[nodes[3,i]:nodes[4,i]]`

"""
struct BitCircuitPair{V,M}
    layers::Vector{V}
    nodes::M
    elements::M
    parents::V
end

struct ParamBitCircuitPair{V,M, WPC, WLC}
    pc_bit::BitCircuit{V,M}
    lc_bit::BitCircuit{V,M}
    pair_bit::BitCircuitPair{V,M}
    pc_params::WPC
    lc_params::WLC
end

function ParamBitCircuitPair(pc::ProbCircuit, lc::LogisticCircuit; Float=Float32)
    pc_thetas::Vector{Float} = Vector{Float}()
    lc_thetas::Vector{Vector{Float}} = Vector{Vector{Float}}()

    sizehint!(pc_thetas, num_edges(pc))

    pc_cache = Dict{DAG, NodeId}() # only for sum nodes
    lc_cache = Dict{DAG, NodeId}() # only for sum nodes

    lc_num_classes = num_classes(lc);
    pc_on_decision(n, cs, layer_id, decision_id, first_element, last_element) = begin
        if isnothing(n) # this decision node is not part of the PC
            push!(pc_thetas, zero(Float))
        else
            pc_cache[n] = decision_id
            append!(pc_thetas, n.log_probs)
        end
    end

    lc_on_decision(m, cs, layer_id, decision_id, first_element, last_element) = begin
        if isnothing(m)
            throw("here, some node is not part of the logistic circuit")
            # push!(lc_thetas, zeros(Float32, lc_num_classes))
        else
            lc_cache[m] = decision_id
            for theta in eachrow(m.thetas)
                push!(lc_thetas, theta)
            end
        end
    end

    pbc_callback(n, m, results, layer_id, last_dec_id, first_el_id, last_el_id) = begin
        nothing
    end

    pc_bit = BitCircuit(pc, num_variables(pc); on_decision=pc_on_decision)
    lc_bit = BitCircuit(lc, num_variables(pc); on_decision=lc_on_decision)
    bcp = BitCircuitPair(pc, lc; on_sum_callback = pbc_callback, pc_cache, lc_cache)

    lc_thetas_reshaped = permutedims(hcat(lc_thetas...), (2, 1))
    ParamBitCircuitPair(pc_bit, lc_bit, bcp, pc_thetas, lc_thetas_reshaped)
end

function BitCircuitPair(pc::ProbCircuit, lc::LogisticCircuit;
    on_sum_callback=noop,
    pc_cache=nothing, lc_cache=nothing, property_check=false)

    if property_check
        @assert num_variables(pc) == num_variables(lc)
        vtree1::Vtree = infer_vtree(pc);
        @assert respects_vtree(lc, vtree1);
    end

    num_features = num_variables(pc)
    num_leafs = 4*num_features
    layers::Vector{Vector{NodePairId}} = Vector{NodePairId}[collect(1:num_leafs)]
    nodes::Vector{NodePairId} = zeros(NodePairId, NODES_LENGTH*num_leafs)
    elements::Vector{NodePairId} = NodePairId[]
    parents::Vector{Vector{NodePairId}} = Vector{NodePairId}[NodePairId[] for i = 1:num_leafs]
    last_dec_id::NodePairId = num_leafs
    last_el_id::NodePairId = zero(NodePairId)

    cache = Dict{Pair{DAG, DAG}, NodePairIds}()

    func(n,m) = begin
        throw("Unsupported pair of nodes!! $n, $m")
    end

    function func(n::Union{PlainProbLiteralNode, StructProbLiteralNode},
                    m::LogisticLiteralNode)::NodePairIds
        get!(cache, Pair(n, m)) do
            @assert variable(n) == variable(m)
            if ispositive(n) && ispositive(m)
                SumNodePairIds(one(NodePairId), NodePairId(variable(n)))
            elseif !ispositive(n) && ispositive(m)
                SumNodePairIds(one(NodePairId), NodePairId(variable(n) + num_features))
            elseif ispositive(n) && !ispositive(m)
                SumNodePairIds(one(NodePairId), NodePairId(variable(n) + 2*num_features))
            elseif !ispositive(n) && !ispositive(m)
                SumNodePairIds(one(NodePairId), NodePairId(variable(n) + 3*num_features))
            end
        end
    end

    function func(n::Union{PlainProbLiteralNode, StructProbLiteralNode},
                    m::Logistic⋁Node)::NodePairIds

        @assert num_children(m) == 1
        func(n, children(m)[1])
    end

    function func(n::Union{PlainMulNode, StructMulNode},
            m::Logistic⋀Node)::NodePairIds

        get!(cache, Pair(n, m)) do
            @assert num_children(n) == 2
            @assert num_children(m) == 2

            ll = func(children(n)[1], children(m)[1])
            rr = func(children(n)[2], children(m)[2])
            ProdNodePairIds(ll, rr)
        end
    end

    function func(n::Union{PlainSumNode, StructSumNode},
                    m::Logistic⋁Node)::NodePairIds

        get!(cache, Pair(n, m)) do
            # First process the children
            results::Vector{Vector{NodePairIds}} = NodePairIds[]
            for i in 1:num_children(n)
                push!(results, NodePairIds[])
                for j in 1:num_children(m)
                    push!(results[i], func(children(n)[i], children(m)[j]))
                end
            end

            # Do not move this, otherwise el_ids would be wrong.
            first_el_id::NodePairId = last_el_id + one(NodePairId)
            layer_id::NodePairId = zero(NodePairId)
            push!(parents, NodePairId[])

            # Add processed children to bitcircuitpair
            for i in 1:num_children(n)
                for j in 1:num_children(m)
                    cur = results[i][j];
                    layer_id = max(layer_id, cur.layer_id)
                    last_el_id += one(NodePairId)

                    if typeof(cur) == ProdNodePairIds
                        push!(elements, last_dec_id, cur.left_left_id, cur.right_right_id);
                        @inbounds push!(parents[cur.left_left_id], last_el_id)
                        @inbounds push!(parents[cur.right_right_id], last_el_id)
                    else
                        !property_check &&  @assert children(m)[1] isa LogisticLeafNode
                        push!(elements, last_dec_id, cur.node_id, cur.node_id);
                        @inbounds push!(parents[cur.node_id], last_el_id)
                        @inbounds push!(parents[cur.node_id], last_el_id)
                    end
                end
            end
            layer_id += one(NodePairId)
            length(layers) < layer_id && push!(layers, NodePairId[])

            last_dec_id::NodePairId += one(NodePairId)

            PC_ID = pc_cache !== nothing ? pc_cache[n] : zero(NodePairId)
            LC_ID = lc_cache !== nothing ? lc_cache[m] : zero(NodePairId)

            push!(nodes, first_el_id, last_el_id, zero(NodePairId), zero(NodePairId), PC_ID, LC_ID)
            push!(layers[layer_id], last_dec_id)
            on_sum_callback(n, m, results, layer_id, last_dec_id, first_el_id, last_el_id)
            SumNodePairIds(layer_id, last_dec_id)
        end
    end

    # Call on roots
    func(pc, lc)

    nodes_m = reshape(nodes, NODES_LENGTH, :)
    elements_m = reshape(elements, ELEMENTS_LENGTH, :)
    parents_m = Vector{NodePairId}(undef, size(elements_m,2)*2)

    last_parent = zero(NodePairId)[]
    @assert last_dec_id == size(nodes_m,2) == size(parents,1)
    @assert sum(length, parents) == length(parents_m)
    for i in 1:last_dec_id-1
        if !isempty(parents[i])
            nodes_m[3,i] = last_parent + one(NodePairId)
            parents_m[last_parent + one(NodePairId):last_parent + length(parents[i])] .= parents[i]
            last_parent += length(parents[i])
            nodes_m[4,i] = last_parent
        else
            !property_check && @assert i <= num_leafs "Only root and leaf nodes can have no parents: $i <= $num_leafs"
        end
    end
    return BitCircuitPair(layers, nodes_m, elements_m, parents_m)
end

######################
## Helper Functions ##
######################

import LogicCircuits: num_nodes, num_elements, num_features, num_leafs, nodes, elements
import LogicCircuits: to_gpu, to_cpu, isgpu #extend

# most of this are identical with bitcircuit maybe make the BitCircuitPair a subtype of BitCircuit?
nodes(c::BitCircuitPair) = c.nodes
elements(c::BitCircuitPair) = c.elements

num_nodes(c::BitCircuitPair) = size(c.nodes, 2)

num_elements(c::BitCircuitPair) = size(c.elements, 2)

to_gpu(c::BitCircuitPair) =
    BitCircuitPair(map(to_gpu, c.layers), to_gpu(c.nodes), to_gpu(c.elements), to_gpu(c.parents))

to_cpu(c::BitCircuitPair) =
    BitCircuitPair(map(to_cpu, c.layers), to_cpu(c.nodes), to_cpu(c.elements), to_cpu(c.parents))

isgpu(c::BitCircuitPair{<:CuArray,<:CuArray}) = true
isgpu(c::BitCircuitPair{<:Array,<:Array}) = false

#############################
## Param Helper functions ###
#############################

pc_params(c::ParamBitCircuitPair) = c.pc_params
lc_params(c::ParamBitCircuitPair) = c.lc_params

num_nodes(c::ParamBitCircuitPair) = num_nodes(c.pair_bit)
num_elements(c::ParamBitCircuitPair) = num_elements(c.pair_bit)
num_features(c::ParamBitCircuitPair) = num_features(c.pair_bit)
num_leafs(c::ParamBitCircuitPair) = num_leafs(c.pair_bit)

nodes(c::ParamBitCircuitPair) = nodes(c.pair_bit)
elements(c::ParamBitCircuitPair) = elements(c.pair_bit)

to_gpu(c::ParamBitCircuitPair) =
    ParamBitCircuitPair(to_gpu(c.pc_bit), to_gpu(c.lc_bit),
        to_gpu(c.pair_bit), to_gpu(c.pc_params), to_gpu(c.lc_params))

to_cpu(c::ParamBitCircuitPair) =
    ParamBitCircuitPair(to_cpu(c.pc_bit), to_cpu(c.lc_bit), to_cpu(c.pair_bit), to_cpu(c.pc_params), to_cpu(c.lc_params))

isgpu(c::ParamBitCircuitPair) =
    isgpu(c.pc_bit) && isgpu(c.lc_bit) && isgpu(c.pair_bit) && isgpu(c.pc_params) && isgpu(c.lc_params)

max_stars_repo_name: Juice-jl/DiscriminativeCircuits

max_issues_repo_name: Juice-jl/DiscriminativeCircuits

max_forks_repo_name: Juice-jl/DiscriminativeCircuits


text:
module Parameters

export FreeParameters, lognormal, ScaledLogitNormal

using Oceananigans.Architectures: CPU, arch_array, architecture
using Oceananigans.Utils: prettysummary
using Oceananigans.TurbulenceClosures: AbstractTurbulenceClosure
using Oceananigans.TurbulenceClosures: AbstractTimeDiscretization, ExplicitTimeDiscretization

using Printf
using Distributions
using DocStringExtensions
using LinearAlgebra

using SpecialFunctions: erfinv
using Distributions: AbstractRNG, ContinuousUnivariateDistribution

#####
##### Priors
#####

"""
    lognormal(; mean, std)

Return `Lognormal` distribution parameterized by
the distribution `mean` and standard deviation `std`.

Notes
=====

A variate `X` is `LogNormal` distributed if

```math
\\log(X) ∼ 𝒩(μ, σ²) ,
```

where ``𝒩(μ, σ²)`` is the `Normal` distribution with mean ``μ``
and variance ``σ²``.

The `mean` and variance ``s²`` (where ``s`` is the standard
deviation or `std`) are related to the parameters ``μ``
and ``σ²`` via

```math
 m = \\exp(μ + σ² / 2),
```
```math
s² = [\\exp(σ²) - 1] m².
```

These formula allow us to calculate ``μ`` and ``σ`` given
``m`` and ``s²``, since rearranging the formula for ``s²``
gives

```math
\\exp(σ²) = m² / s² + 1
```

which then yields

```math
σ = \\sqrt{\\log(m² / s² + 1)}.
```

We then find that

```math
μ = \\log(m) - σ² / 2 .
```

See also
[wikipedia](https://en.wikipedia.org/wiki/Log-normal_distribution#Generation_and_parameters).
"""
function lognormal(; mean, std)
    k = std^2 / mean^2 + 1 # intermediate variable
    σ = sqrt(log(k))
    μ = log(mean) - σ^2 / 2
    return LogNormal(μ, σ)
end

struct ScaledLogitNormal{T} <: ContinuousUnivariateDistribution
    μ :: T
    σ :: T
    lower_bound :: T
    upper_bound :: T

    ScaledLogitNormal{T}(μ, σ, L, U) where T = new{T}(T(μ), T(σ), T(L), T(U))
end

"""Return a logit-normally distributed variate given the normally-distributed variate `X`."""
normal_to_scaled_logit_normal(L, U, X) = L + (U - L) / (1 + exp(X))

"""Return a normally-distributed variate given the logit-normally distributed variate `Y`."""
scaled_logit_normal_to_normal(L, U, Y) = log((U - Y) / (Y - L))

Base.rand(rng::AbstractRNG, d::ScaledLogitNormal) =
    normal_to_scaled_logit_normal(d.lower_bound, d.upper_bound, rand(rng, Normal(d.μ, d.σ)))

unit_normal_std(mass) = 1 / (2 * √2 * erfinv(mass))

"""
    ScaledLogitNormal([FT=Float64;] bounds=(0, 1), mass=0.5, interval=nothing, μ=nothing, σ=nothing)

Return a `ScaledLogitNormal` distribution with compact support within `bounds`.

`interval` is an optional 2-element tuple or Array. When specified,
the parameters `μ` and `σ` of the underlying `Normal` distribution
are calculated so that `mass` fraction of the probability density
lies within `interval`.

If `interval` is not specified, then `μ=0` and `σ=1` by default.

Notes
=====

`ScaledLogitNormal` is a four-parameter distribution
generated by the transformation

```math
Y = L + (U - L) / [1 + \\exp(X)],
```

of the normally-distributed variate ``X ∼ 𝒩(μ, σ)``. The four parameters
governing the distribution of ``Y`` are thus

- ``L``:  lower bound (0 for the `LogitNormal` distribution)
- ``U``:  upper bound (1 for the `LogitNormal` distribution)
- ``μ``:  mean of the underlying `Normal` distribution
- ``σ²``: variance of the underlying `Normal` distribution
"""
function ScaledLogitNormal(FT=Float64; bounds=(0, 1), mass=0.5, interval=nothing, μ=nothing, σ=nothing)
    L, U = bounds

    if isnothing(interval) # use default μ=0 and σ=1 if not set
        isnothing(μ) && (μ = 0)
        isnothing(σ) && (σ = 1)

    elseif !isnothing(interval) # try to compute μ and σ

        Li, Ui = interval

        # User friendliness
        (!isnothing(μ) || !isnothing(σ)) && @warn "Using interval and mass to determine μ and σ."
        0 < mass < 1 || throw(ArgumentError("Mass must lie between 0 and 1."))
        Li > L && Ui < U || throw(ArgumentError("Interval limits must lie between `bounds`."))

        # Compute lower and upper limits of midspread in unconstrained space
        #
        # Note that the _lower_ bound in unconstrained space is associated with the
        # _upper_ bound in constrained space, and vice versa.
        L̃i = scaled_logit_normal_to_normal(L, U, Ui)
        Ũi = scaled_logit_normal_to_normal(L, U, Li)

        μ = (Ũi + L̃i) / 2

        # Note that the mass beneath a half-width `δ` of the
        # standard Normal distribution is
        #
        # mass = 2 / √(2π) ∫₀ᵟ exp(-x^2 / 2) dx
        #      = erf(δ / √2)
        #
        # For an `interval = (Ũi, L̃i)` of the normal distribution,
        # the non-dimensional half-width is
        #
        # δ = (Ũi - L̃i) / 2σ
        #
        # where σ is the distribution's standard deviation.
        # We then find
        #
        # erfinv(mass) = (Ũi - L̃i) / (2 * √2 * σ) ,
        #
        # and rearranging to solve for σ yields
        #
        σ = (Ũi - L̃i) / (2 * √2 * erfinv(mass))
    end

    return ScaledLogitNormal{FT}(μ, σ, L, U)
end

# Calculate the prior in unconstrained space given a prior in constrained space
unconstrained_prior(Π::LogNormal)         = Normal(Π.μ / abs(Π.μ), Π.σ / abs(Π.μ))
unconstrained_prior(Π::Normal)            = Normal(0, 1)
unconstrained_prior(Π::ScaledLogitNormal) = Normal(Π.μ, Π.σ)

"""
    transform_to_unconstrained(Π, Y)

Transform the "constrained" (physical) variate `Y` into it's
unconstrained (normally-distributed) counterpart `X` through the
forward map associated with `Π`.

If some mapping between ``Y`` and the normally-distributed ``X`` is
defined via

```math
Y = g(X).
```

Then `transform_to_unconstrained` is the inverse ``X = g^{-1}(Y)``.
The change of variables ``g(X)`` determines the distribution `Π` of `Y`.

Example
=======

The logarithm of a `LogNormal(μ, σ)` distributed variate is normally-distributed,
such that the forward trasform ``f ≡ \\exp``,

```math
Y = \\exp(X),
```

and the inverse trasnform is the natural logarithm ``f^{-1} ≡ \\log``,

```math
\\log(Y) = X ∼ 𝒩(μ, σ).
```
"""
transform_to_unconstrained(Π::Normal,    Y) = (Y - Π.μ) / Π.σ
transform_to_unconstrained(Π::LogNormal, Y) = log(Y^(1 / abs(Π.μ))) # log(Y) / abs(Π.μ)

transform_to_unconstrained(Π::ScaledLogitNormal, Y) =
    scaled_logit_normal_to_normal(Π.lower_bound, Π.upper_bound, Y)

"""
    transform_to_constrained(Π, X)

Transform an "unconstrained", normally-distributed variate `X`
to "constrained" (physical) space via the map associated with
the distribution `Π` of `Y`.
"""
transform_to_constrained(Π::Normal, X)    = X * Π.σ + Π.μ
transform_to_constrained(Π::LogNormal, X) = exp(X * abs(Π.μ))
transform_to_constrained(Π::ScaledLogitNormal, X) =
    normal_to_scaled_logit_normal(Π.lower_bound, Π.upper_bound, X)

# Convenience vectorized version
transform_to_constrained(priors::NamedTuple, X::AbstractVector) =
    NamedTuple(name => transform_to_constrained(priors[name], X[i])
               for (i, name) in enumerate(keys(priors)))

# Convenience matrixized version assuming particles vary on 2nd dimension
transform_to_constrained(priors::NamedTuple, X::AbstractMatrix) =
    [transform_to_constrained(priors, X[:, k]) for k = 1:size(X, 2)]

function inverse_covariance_transform(Π, X, covariance)
    diag = [covariance_transform_diagonal(Π[i], X[i]) for i=1:length(Π)]
    dT = Diagonal(diag)
    return dT * covariance * dT'
end

covariance_transform_diagonal(::LogNormal, X) = exp(X)
covariance_transform_diagonal(::Normal, X)    = 1
covariance_transform_diagonal(Π::ScaledLogitNormal, X) = - (Π.upper_bound - Π.lower_bound) * exp(X) / (1 + exp(X))^2

#####
##### Free parameters
#####

"""
    struct FreeParameters{N, P, D}

A container for free parameters that includes the parameter names and their
corresponding prior distributions.

$(FIELDS)
"""
struct FreeParameters{N, P, D}
    "free parameters"
    names :: N
    "prior distributions for free parameters"
    priors :: P
    "dependent parameters"
    dependent_parameters :: D
end

"""
    FreeParameters(priors; names = Symbol.(keys(priors)), dependent_parameters=NamedTuple())

Return named `FreeParameters` with priors. Free parameter `names` are inferred from
the keys of `priors` if not provided. Optionally, `dependent_parameters` are prescribed
as a `NamedTuple` whose keys are the names of "additional" parameters, and whose values
are functions that return those parameters given a vector of free parameters in `names`.

Example
=======

```jldoctest
julia> using Distributions, ParameterEstimocean

julia> priors = (ν = Normal(1e-4, 1e-5), κ = Normal(1e-3, 1e-5))
(ν = Normal{Float64}(μ=0.0001, σ=1.0e-5), κ = Normal{Float64}(μ=0.001, σ=1.0e-5))

julia> free_parameters = FreeParameters(priors)
FreeParameters with 2 parameters
├── names: (:ν, :κ)
├── priors: Dict{Symbol, Any}
│   ├── ν => Normal{Float64}(μ=0.0001, σ=1.0e-5)
│   └── κ => Normal{Float64}(μ=0.001, σ=1.0e-5)
└── dependent parameters: Dict{Symbol, Any}

julia> c(p) = p.ν + p.κ # compute a third dependent parameter `c` as a function of `ν` and `κ`
c (generic function with 1 method)

julia> free_parameters_with_a_dependent = FreeParameters(priors, dependent_parameters=(; c))
FreeParameters with 2 parameters and 1 dependent parameter
├── names: (:ν, :κ)
├── priors: Dict{Symbol, Any}
│   ├── ν => Normal{Float64}(μ=0.0001, σ=1.0e-5)
│   └── κ => Normal{Float64}(μ=0.001, σ=1.0e-5)
└── dependent parameters: Dict{Symbol, Any}
    └── c => c (generic function with 1 method)
```
"""
function FreeParameters(priors; names = Symbol.(keys(priors)), dependent_parameters=NamedTuple())
    priors = NamedTuple(name => priors[name] for name in names)
    return FreeParameters(Tuple(names), priors, dependent_parameters)
end

Base.summary(fp::FreeParameters) = "$(fp.names)"

function prior_show(io, priors, name, prefix, width)
    print(io, @sprintf("%s %s => ", prefix, lpad(name, width, " ")))
    show(io, priors[name])
    return nothing
end

function dependent_parameter_show(io, dependent_parameters, name, prefix, width)
    print(io, @sprintf("%s %s => ", prefix, lpad(name, width, " ")))
    print(io, prettysummary(dependent_parameters[Symbol(name)]))
    return nothing
end

parameter_str(N) = N>1 ? "parameters" : "parameter"

function Base.show(io::IO, p::FreeParameters)
    Np, Nd = length(p), length(p.dependent_parameters)

    free_parameters_summary = "FreeParameters with $Np " * parameter_str(Np)

    title = Nd > 0 ?
            free_parameters_summary * " and $Nd dependent " * parameter_str(Nd) :
            free_parameters_summary

    print(io, title, '\n',
              "├── names: $(p.names)", '\n',
              "├── priors: Dict{Symbol, Any}")

    maximum_name_length = maximum([length(string(name)) for name in p.names])

    for (i, name) in enumerate(p.names)
        prefix = i == length(p.names) ? "│   └──" : "│   ├──"
        print(io, '\n')
        prior_show(io, p.priors, name, prefix, maximum_name_length)
    end

    print(io, '\n')

    print(io, "└── dependent parameters: Dict{Symbol, Any}")

    if !isempty(p.dependent_parameters)
        maximum_name_length = maximum([length(string(name)) for name in p.dependent_parameters])

        for (i, name) in enumerate(p.dependent_parameters)
            prefix = i == length(p.dependent_parameters) ? "    └──" : "    ├──"
            print(io, '\n')
            dependent_parameter_show(io, p.dependent_parameters, name, prefix, maximum_name_length)
        end
    end

    return nothing
end

Base.length(p::FreeParameters) = length(p.names)

function build_parameters_named_tuple(p::FreeParameters, free_θ)
    if free_θ isa Dict # convert to NamedTuple with
        free_θ = NamedTuple(name => free_θ[name] for name in p.names)
    elseif !(free_θ isa NamedTuple) # mostly likely a Vector: convert to NamedTuple with
        free_θ = NamedTuple{p.names}(Tuple(free_θ))
    end

    # Compute dependent parameters
    dependent_names = keys(p.dependent_parameters)
    maps = values(p.dependent_parameters)
    dependent_θ = NamedTuple(name => maps[name](free_θ) for name in dependent_names)

    return merge(dependent_θ, free_θ) # prioritize free_θ
end

#####
##### Setting parameters
#####

const ParameterValue = Union{Number, AbstractArray}

"""
    construct_object(specification_dict, parameters; name=nothing, type_parameter=nothing)

    construct_object(d::ParameterValue, parameters; name=nothing)

Return a composite type object whose properties are prescribed by the `specification_dict`
dictionary. All parameter values are given the values in `specification_dict` *unless* they
are included as a parameter name-value pair in the named tuple `parameters`, in which case
the value in `parameters` is asigned.

The `construct_object` is recursively called upon every property that is included in `specification_dict`
until a property with a numerical value is reached. The object's constructor name must be
included in `specification_dict` under key `:type`.

Example
=======

```jldoctest; filter = [r".*Dict{Symbol.*", r".*:type       => Closure.*", r".*:c          => 3.*", r".*:subclosure => Dict{Symbol.*"]
julia> using ParameterEstimocean.Parameters: construct_object, dict_properties, closure_with_parameters

julia> struct Closure; subclosure; c end

julia> struct ClosureSubModel; a; b end

julia> sub_closure = ClosureSubModel(1, 2)
ClosureSubModel(1, 2)

julia> closure = Closure(sub_closure, 3)
Closure(ClosureSubModel(1, 2), 3)

julia> specification_dict = dict_properties(closure)
Dict{Symbol, Any} with 3 entries:
  :type       => Closure
  :c          => 3
  :subclosure => Dict{Symbol, Any}(:a=>1, :b=>2, :type=>ClosureSubModel)

julia> new_closure = construct_object(specification_dict, (a=2.1,))
Closure(ClosureSubModel(2.1, 2), 3)

julia> another_new_closure = construct_object(specification_dict, (b=π, c=2π))
Closure(ClosureSubModel(1, π), 6.283185307179586)
```
"""
construct_object(d::ParameterValue, parameters; name=nothing) =
    name ∈ keys(parameters) ? getproperty(parameters, name) : d

function construct_object(specification_dict, parameters;
                          name=nothing, type_parameter=nothing)

    type = Constructor = specification_dict[:type]
    kwargs_vector = [construct_object(specification_dict[name], parameters; name)
                        for name in fieldnames(type) if name != :type]

    return isnothing(type_parameter) ? Constructor(kwargs_vector...) :
                                       Constructor{type_parameter}(kwargs_vector...)
end

"""
    dict_properties(object)

Return a dictionary with all properties of an `object` and their values, including the
`object`'s type name. If any of the `object`'s properties is not a numerical value but
instead a composite type, then `dict_properties` is called recursively on that `object`'s
property returning a dictionary with all properties of that composite type. Recursion
ends when properties of type `ParameterValue` are found.
"""
function dict_properties(object)
    p = Dict{Symbol, Any}(n => dict_properties(getproperty(object, n)) for n in propertynames(object))
    p[:type] = typeof(object).name.wrapper

    return p
end

dict_properties(object::ParameterValue) = object

"""
    closure_with_parameters(closure, parameters)

Return a new object where for each (`parameter_name`, `parameter_value`) pair
in `parameters`, the value corresponding to the key in `closure` that matches
`parameter_name` is replaced with `parameter_value`.

Example
=======

Create a placeholder `Closure` type that includes a parameter `c` and a sub-closure
with two parameters: `a` and `b`. Then construct a closure with values `a, b, c = 1, 2, 3`.

```jldoctest closure_with_parameters
julia> struct Closure; subclosure; c end

julia> struct ClosureSubModel; a; b end

julia> sub_closure = ClosureSubModel(1, 2)
ClosureSubModel(1, 2)

julia> closure = Closure(sub_closure, 3)
Closure(ClosureSubModel(1, 2), 3)
```

Providing `closure_with_parameters` with a named tuple of parameter names and values,
and a recursive search in all types and subtypes within `closure` is done and whenever
a parameter is found whose name exists in the named tuple we provided, its value is
then replaced with the value provided.

```jldoctest closure_with_parameters
julia> new_parameters = (a = 12, d = 7)
(a = 12, d = 7)

julia> using ParameterEstimocean.Parameters: closure_with_parameters

julia> closure_with_parameters(closure, new_parameters)
Closure(ClosureSubModel(12, 2), 3)
```
"""
closure_with_parameters(closure, parameters) = construct_object(dict_properties(closure), parameters)

closure_with_parameters(closure::AbstractTurbulenceClosure{ExplicitTimeDiscretization}, parameters) =
    construct_object(dict_properties(closure), parameters, type_parameter=nothing)

closure_with_parameters(closure::AbstractTurbulenceClosure{TD}, parameters) where {TD <: AbstractTimeDiscretization} =
    construct_object(dict_properties(closure), parameters; type_parameter=TD)

closure_with_parameters(closures::Tuple, parameters) =
    Tuple(closure_with_parameters(closure, parameters) for closure in closures)

"""
    update_closure_ensemble_member!(closures, p_ensemble, parameters)

Use `parameters` to update the `p_ensemble`-th closure from and array of `closures`.
The `p_ensemble`-th closure corresponds to ensemble member `p_ensemble`.
"""
update_closure_ensemble_member!(closure, p_ensemble, parameters) = nothing

update_closure_ensemble_member!(closures::AbstractVector, p_ensemble, parameters) =
    closures[p_ensemble] = closure_with_parameters(closures[p_ensemble], parameters)

function update_closure_ensemble_member!(closures::AbstractMatrix, p_ensemble, parameters)
    for j in 1:size(closures, 2) # Assume that ensemble varies along first dimension
        closures[p_ensemble, j] = closure_with_parameters(closures[p_ensemble, j], parameters)
    end

    return nothing
end

function update_closure_ensemble_member!(closure_tuple::Tuple, p_ensemble, parameters)
    for closure in closure_tuple
        update_closure_ensemble_member!(closure, p_ensemble, parameters)
    end
    return nothing
end

"""
    new_closure_ensemble(closures, θ, arch=CPU())

Return a new set of `closures` in which all closures that have free parameters are updated.
Closures with free parameters are expected as `AbstractArray` of `TurbulenceClosures`, and
this allows `new_closure_ensemble` to go through all closures in `closures` and only update
the parameters for the any closure that is of type `AbstractArray`. The `arch`itecture
(`CPU()` or `GPU()`) defines whethere `Array` or `CuArray` is returned.
"""
function new_closure_ensemble(closures::AbstractArray, θ, arch)
    cpu_closures = arch_array(CPU(), closures)

    for (p, θp) in enumerate(θ)
        update_closure_ensemble_member!(cpu_closures, p, θp)
    end

    return arch_array(arch, cpu_closures)
end
new_closure_ensemble(closures::Tuple, θ, arch) =
    Tuple(new_closure_ensemble(closure, θ, arch) for closure in closures)

new_closure_ensemble(closure, θ, arch) = closure

end # module


text:
mutable struct RandomSampleGrid <: SVDD.InitializationStrategyC
    solver
    C_range
    pools
    selection_metric
end

function SVDD.calculate_C(model, strategy::RandomSampleGrid, kernel)
    m = deepcopy(model)

    init_strategy = FixedParameterInitialization(kernel, 1.0)
    SVDD.initialize!(m, init_strategy)
    model_range = typeof(model)[]

    for C in strategy.C_range
        @debug "Training with C=$c"
        set_param!(m, Dict(:C => C))
        set_adjust_K!(m, true)
        SVDD.fit!(m, strategy.solver)
        push!(model_range, m)
        m = deepcopy(m)
    end

    idx_pools = (strategy.pools .!= :U)
    labels = ifelse.(strategy.pools[idx_pools] .== :Lin, :inlier, :outlier)

    class_range = map(x -> SVDD.classify.(SVDD.predict(x, m.data[:, idx_pools])), model_range)
    cms_pools = map(x -> ConfusionMatrix(x, labels), class_range)
    qm = strategy.selection_metric.(cms_pools)
    qm_best_idx = findmax(qm)[2]
    C = strategy.C_range[qm_best_idx]
    return C
end

max_stars_repo_name: holtri/lama

max_issues_repo_name: holtri/lama

max_forks_repo_name: holtri/lama


text:
info("Test read header info")

f = HTSLIB.sam_open(joinpath(Pkg.dir("HTSLIB"),"test/data/100.sam"),"r")
phdr = f.phdr
#@show HTSLIB.strptr(phdr)
HTSLIB.close(f)

info("succefully read header info")


text:
using DrWatson
quickactivate(@__DIR__)
@show projectdir()

## ============================================================================
@info "load packages"
using CovidSurvey
using Random, Dates, Turing
using Base.Threads
using StatsPlots
using PrettyTables
using DataFrames
using BSON

using Memoization
using ReverseDiff
setadbackend(:reversediff)
Turing.setrdcache(true)

plotlyjs()
@info "number of threads available: $(nthreads())"
## ============================================================================

ps = ( warmup = 100,
    model  = "intdeaths",
    steps  = 100,
    seed   = nothing,
    observ = nothing,
    chains = 1,
    preds = nothing,#"CF,CC,CR,CS",#"CS,MG,MA"
    semipara=true,
    rwstep = 7,
    include= true,
    std = false,
)
## ============================================================================
# load data
@info "load data"

data_params = (
      observationsend  = ps.observ
    , predictors        = ps.preds|> CovidSurvey.parse_predictors
    # , cases_start       = ps.cases
    # , hospitmodel       = Regional.HospitInit1(obs_stop="2020-07-01")
    , rwstep            = ps.rwstep
    , epidemicstart     = Date("2020-08-01")
    , numimpute         = 20
    , include_early_dynamic = ps.include
    , link              = KLogistic(4.5)
    , invlink           = KLogit(4.5)
    , predictive        = false
    , covariates_kwargs = Dict(
      :semiparametric => ps.semipara,
      :fname          => [projectdir("data/contacts/int/smoothed_contact_rates_region=$(r).csv" ) for r in International.regions],#[projectdir("data", "mean_contact_rates_region=$r.csv" ) for r in 1:Regional.nregions],
      :shift          => -1,
      :conditions     => :date => x->x>Date("2020-11-10"), # >= 1000 cases / day
      # :datecol        => "date",
      :aggregation    => ps.rwstep,
      :standardize    => ps.std,
      :normalize      => true,
      :mobility       => [projectdir("data/mobility/mobility_country=$(r).csv") for r in International.regions]
      )
    )

data = International.load_data(; data_params... )
d = turing_data = data.turing_data;
## ============================================================================
@info "prepare model"

model = Regional.model_deaths
model = Regional.model_international_deaths
#!isnothing(ps.seed) && Random.seed!(ps.seed);
@info ps
m = model(turing_data; vectorize=true)
m();
## ============================================================================
@info "start sampling"

@time chain = let #460s; 366s
    thinning = 1
    if ps.chains > 1
        sample(m, NUTS(ps.warmup, 0.9; max_depth=6), MCMCThreads(), ps.steps + ps.warmup, ps.chains; progress=true, thinning)
    else
        sample(m, NUTS(ps.warmup, 0.9; max_depth=6), ps.steps + ps.warmup; progress=true, thinning) #; max_depth=15
    end
end

divervences = sum( chain[ps.warmup+1:end, :numerical_error, :] )
@show divervences

## ==========================================================================
@info "save"
fdir = projectdir("reports/tmp", "$(today())")
fname = Regional.save_results(fdir, ps, data, chain)
## ==========================================================================
@info "post processing"
Regional.postprocessing(fname; plot_results=false)
## ==========================================================================
## ==========================================================================
@info "post processing"
fdir, ps, ignores = Regional.parse_fname(fname)
p = Regional.PostProcessing2(fdir, ps, ignores, [], fname)
## ==========================================================================
@info "plot chain"
plot_results = true
Regional.plot_chains(p; plot_results)

## ==========================================================================
@info "meanplot"
Regional.plot_means(p; plot_results)
## ==========================================================================
# perform diagnostics
p = Regional.skip_warmup(p)
Regional.diagnostics(p)
data.cases
## ==========================================================================
@info "make predictions"
gp = Regional.generate_posterior(p);
## ==========================================================================
@info "plot regions"
Regional.plot_regions(p, gp; plot_results=true)
## ==========================================================================
@info "plot predictors"
# pgfplotsx()
# default(titlefontsize = 20, legendfontsize = 18, labelfontsize = 18, guidefontsize = 18, tickfontsize = 12, framestyle = :zerolines, yminorgrid = true)
Regional.plot_effects(p, gp; plot_results , grouped = false, effect_on_Rt = 0.)
Regional.plot_effects(p, gp; plot_results, grouped = false, effect_on_Rt = 2.)
Regional.plot_effects(p, gp; plot_results, grouped = false, effect_on_Rt = -0.5)

Regional.plot_effects(p, gp; plot_results, grouped = true, effect_on_Rt = 0.)
Regional.plot_effects(p, gp; plot_results, grouped = true, effect_on_Rt = 2.)
Regional.plot_effects(p, gp; plot_results, grouped = true, effect_on_Rt = -0.5)

## ===================================================================

TV = Vector{Float64}
V = Float64
vectorize = true
θ = turing_data
@unpack num_observations, num_total_days, num_regions, num_rt_steps, deathmodel, num_covariates, rwscale, include_early_dynamic,init_infected = θ
  # If we don't want to predict the future, we only need to compute up-to time-step `num_obs_countries[m]`
num_time_steps = num_observations

############# 2.) time varying reproduction number

R0_std     = truncated(Normal(.5, .2), 0.01, Inf) |> rand
R0s        = filldist(truncated(Normal(0.8, R0_std), 0., 4.5), num_regions) |> rand
σ_rt       = truncated(Normal(0.3*rwscale, .02*rwscale), 0, .5*rwscale) |> rand

n = first( num_rt_steps )
latent_Rts_z = filldist( CovidSurvey.RandomWalk( n ), num_regions) |> rand
latent_Rts = CovidSurvey.rescale( latent_Rts_z, σ_rt, d.invlink.(R0s) )

Rts = TV[TV(undef, num_time_steps[m]) for m in 1:num_regions]
if num_covariates > 0
    grouped_effect = filldist( Laplace(0,0.2), num_covariates) |> rand
    effect_std = filldist( GammaMeanCv(0.1, 0.5), num_covariates) |> rand
    effects_z = filldist( MvNormal( num_covariates, 1.), num_regions) |> rand
    effects = [ effects_z[:,m] .* effect_std .+ grouped_effect for m in 1:num_regions]

    # grouped_effect = filldist( Exponential(0.2), num_covariates)
    # effect_std = filldist( GammaMeanCv(0.1, 0.5), num_covariates)
    # pooled_effects_distr_raw = Normal.( grouped_effect, effect_std )
    # pooled_effects_distr = arraydist( truncated.( pooled_effects_distr_raw, zero(V), V(Inf) ) )
    # effects = filldist( pooled_effects_distr, num_regions )

    if θ.semiparametric
        CovidSurvey.semiparametric!(Rts, θ, false, latent_Rts, R0s, σ_rt, effects)
    else
        CovidSurvey.mixed!(Rts, θ, false, latent_Rts, R0s, σ_rt, effects)
    end
else
    CovidSurvey.random_walks!(Rts, θ, false, latent_Rts, R0s, σ_rt)
end
plot(latent_Rts_z[:,3], ticks=:native)
plot(θ.link.(latent_Rts[:,1]), ticks=:native)
plot(d.dates[1], Rts[1], ticks=:native)
############ 3.) infection dynamics
ys = arraydist(Exponential.(init_infected)) |> rand

newly_infecteds      = TV[TV(undef, num_time_steps[m]) for m in 1:num_regions]
CovidSurvey.infections!(newly_infecteds, θ, ys, Rts)
plot(d.dates[1], newly_infecteds[1], ticks=:native)

# infections!(newly_infecteds, cumulative_infecteds, effective_Rts, θ, τ, ys, Rts)
########### 4.) derive observables
μ_i2d = truncated(Normal(21.82, 1.01), 18, 24) |> rand
σ_i2d = truncated(Normal(4.26, 5.18), 1, 14) |> rand
idr   = truncated(Normal(8/1000, 1/1000), 1/1000, 15/1000) |> rand
ϕ_d   = truncated(Normal(50, 10), 20, Inf) |> rand

weekdayeffect_simplex = filldist(Dirichlet([7,7,7,7,7,3,1]), num_regions) |> rand
weekdayeffect = TV[TV(undef, 7) for m in 1:num_regions]
CovidSurvey.weekdayeffect!(weekdayeffect, weekdayeffect_simplex)

### holiday effect
holidayeffect = filldist(Beta(1,1), num_regions) |> rand

expected_daily_deaths    = TV[TV(undef, num_time_steps[m]) for m in 1:num_regions]
deaths_observation_model = CovidSurvey.WeekdayHolidayObsModel2(
    deathmodel, μ_i2d, σ_i2d, idr, ϕ_d, expected_daily_deaths, weekdayeffect, holidayeffect, vectorize
)

# deaths_observation_model = CovidSurvey.SimpleObsModel2(deathmodel, μ_i2d, σ_i2d, idr, ϕ_d, expected_daily_deaths, vectorize)

CovidSurvey.expected!(deaths_observation_model, newly_infecteds)
plot(d.dates[1], expected_daily_deaths[1], ticks=:native)

########### 4.) compare model to observations
## 4.1) observe hospitalizations

ℓ  = zero(V)
ℓ += CovidSurvey.logpdf(deaths_observation_model, θ.deaths)

max_stars_repo_name: andreaskoher/CovidSurvey

max_issues_repo_name: andreaskoher/CovidSurvey

max_forks_repo_name: andreaskoher/CovidSurvey