text:
# Estimation problem
# A "problem" is the combination of data and model parameters.
# This naming is used in DynamicHMC. In Stan, this is called the data block

"
$TYPEDEF

SDE estimation problem containing both the data and model parameters needed for parameter
estimation.
"
abstract type Problem{E<:SdeEngine} end

const widest_parint = 1e-3..1e3

const widest_freqint = 1e-3..1e2

const default_N_min = 100

function minval_with_data(x,N)
  xm = median(x)
  if length(x)>N
    xs = sort(x)
    xm = min(xm,xs[N])
  end
  return xm
end

"
$SIGNATURES

Parameter interval for given data. Not tested!
"
function parint(data::Vector{Obs})
  @warn "parint from data untested"
  t = diff([x.t for x in data])
  Tmin = minval_with_data(diff(t),default_N_min)
  Tmax = maximum(t)-minimum(t)
  idata = (1/Tmax)..(1/Tmin)
  i = idata ∩ widest_parint
  return i
end

"
Likelihood problem for Stochastic Differential Equation

$FIELDS
"
struct ProblemLike{E}  <: Problem{E}
  "vector of (time,observation) tuples"
  data::Vector{Obs}
  "Model order"
  p::Int
  "Parameter interval"
  parint::Interval{Float64}
  "Frequency range, used for imaginary part of roots (SdeRoots only)"
  freqint::Interval{Float64}
  "Only allow sorted roots (SdeRoots only)"
  sorted_roots::Bool
  "Standard deviation of `y`. Is pre-computed for use in the likelihood computation"
  σy::Float64
end

function ProblemLike(data;p=1,
    parint=widest_parint,freqint=widest_freqint,sorted_roots=true,engine=default_engine,
    σy = √mean(x->x.y.^2,data))
  return ProblemLike{engine}(data,p,parint,freqint,sorted_roots,σy)
end

function Base.show(io::IO,p::ProblemLike{E}) where {E}
  print(io,"SDE($(p.p))")
  print(io," | a ∈ ",p.parint)
  if sdetype(E)===SdeRoots
    print(io," f ∈ ",p.freqint, " sorted roots: ", p.sorted_roots ? "yes" : "no")
  end
  print(io," | engine : $E")
  Tav = diff(map(x->x.t,p.data)) |> mean
  print(io," || N = ",length(p.data),@sprintf(" Tav = %0.2g",Tav))
end

function bringin(i::AbstractInterval,x,f)
  if !(i>0) error("Currently only works for intervals>0") end
  if x≪i  i.first/f  elseif x≫i  f*i.last else x end
end

"
$SIGNATURES

Is process `proc` in the problem range?
"
function Base.in(proc::SdeProcess{<:Any,<:Sde},problem::Problem)
  return order(proc)==problem.p && all(proc.eq.a.∈problem.parint)
end

function bringin(problem,eq::Sde)
  ain = bringin.(problem.parint,eq.a,0.9)
  return Sde(ain,eq.σϵ)
end

function Base.in(proc::SdeProcess{<:Any,<:SdeRoots},problem::Problem)
  rr,ri = splitroots(proc.eq.r)
  isin = order(proc)==problem.p &&
    all(-rr.∈problem.parint) && all(-real(ri).∈problem.parint) &&
    all(-imag(ri).<2π*problem.freqint.last)
  return isin
end

function bringin(problem,eq::SdeRoots)
  rr,ri = splitroots(eq.r)
  rrin = -bringin.(problem.parint,-rr,0.9)
  riin = map(ri) do ric
    a = -bringin(problem.parint,-real(ric),0.9)
    b = -bringin(problem.freqint,-imag(ric),0.9)
    a+b*im
  end
  rin = combineroots(rrin,riin)
  return SdeRoots(rin)
end

"
$SIGNATURES

Bring SDE process `proc` into the range of `problem` to guarantee
`(bringin(problem,proc) ∈ problem) == true`
"
function bringin(problem,proc::SdeProcess)
  eq = sde(problem,proc.eq)
  eqin = bringin(problem,eq)
  return SdeProcess(eqin)
end

"
$SIGNATURES

`SdeProblemLike` for lower model order `pc`.
"
function lowerorder(p::ProblemLike{E},pc) where {E}
  return ProblemLike{E}(p.data,pc,p.parint,p.freqint,p.sorted_roots,p.σy)
end

## Probabilistic probem ——————————————————————————————————

abstract type Prior end

"Flat prior in the SDE coefficients `a`"
struct FlatCoefPrior<:Prior end

Base.show(io::IO,::FlatCoefPrior) = print(io,"Flat a")

(p::FlatCoefPrior)(::NamedTuple{(:a,)}) = 0

"
Reference Prior derived from regularly sampled SDE(1) process at sampling interval `T`

$FIELDS
"
struct ReferencePrior<:Prior
  "Sampling time interval"
  T::Float64
end

Distributions.logpdf(p::ReferencePrior,θ::NamedTuple) = -1/2*sum(@. log(1-exp(2θ.r_re*p.T)))

# Note: A nested call is required for writing the custom adjoint below
(p::ReferencePrior)(θ::NamedTuple) = logpdf(p,θ)

# Ensure that the adjoint for prior(θ) is real-valued
Zygote.@adjoint (p::ReferencePrior)(θ::NamedTuple{(:r_re,:r_im)}) = begin
  lp,Bc = Zygote.pullback(logpdf,p,θ)
  if isinf(lp) return lp,Bc end
  B = function(l̄)
      p̄,θ̄c = Bc(l̄)
      # Because the prior does not depend on r_im, we know the derivative w.r.t. it is `nothing`
      # By hard coding this, type stability is maintained, hopefully resulting in faster execution.
      θ̄ = (r_re=real.(θ̄c.r_re),r_im=nothing)
      return p̄,θ̄
    end
  return lp,B
end

Zygote.@adjoint (p::ReferencePrior)(θ::NamedTuple{(:r_re,)}) = begin
  lp,Bc = Zygote.pullback(logpdf,p,θ)
  if isinf(lp) return lp,Bc end
  B = function(l̄)
      p̄,θ̄c = Bc(l̄)
      θ̄ = (r_re=real.(θ̄c.r_re),)
      return p̄,θ̄
    end
  return lp,B
end

"""
Probabilistic Stochastic Differential Equation problem.

$FIELDS
"""
struct ProblemProb{E,P<:Prior} <: Problem{E}
  "Likelihood"
  like::ProblemLike{E}
  "Prior"
  prior::P
end

ProblemProb(data;prior=FlatCoefPrior(),kwargs...) = ProblemProb(ProblemLike(data;kwargs...),prior)

Base.show(io::IO,p::ProblemProb) = print(io,p.like," || prior : ",p.prior)

engine(::ProblemLike{E},eq::AbstractSde) where {E} = E(eq)

function Distributions.logpdf(problem::ProblemLike,proc::SdeProcess)
  proc′ = SdeProcess(sde(problem,proc.eq))
  e = engine(problem,proc′.eq)
  lp,_ = logpdf_state(e,proc′.state,problem.data)
  return lp
end

"
$SIGNATURES

Log probability for named parameter tuple `θ` with maximum likelihood
estimate for `σϵ`.
"
function logpdf_ml_σϵ(problem::ProblemLike,θ)
  eq1 = sde(problem,θ,1)
  state1 = statstate(eq1)
  e = engine(problem,eq1)
  lp,σϵ = logpdf_ml_σϵ(e,state1,problem.data)
  return lp,σϵ
end

function (problem::ProblemLike{E})(θ) where{E}
  if sdetype(E)===SdeRoots && problem.sorted_roots
    if !issorted_rootstuple(θ)
      @debug "Rejecting point because roots are not sorted" θ
      return -Inf
    end
  end
  lp = try
    logpdf_ml_σϵ(problem,θ)[1]
  catch e
    @debug "Rejecting point because logpdf_ml_σϵ(problem,$θ) throws $e"
    -Inf
  end
  return lp
end

function safe_problem_adjoint(problem,θ)
  lp,B0 = try
    Zygote.pullback((p,t)->logpdf_ml_σϵ(p,t)[1],problem,θ)
  catch e
    @debug "∂logpdf_ml_σϵ(problem,$θ) throws $e"
    return -Inf,(_) -> nothing
  end
end

# Ensure that the adjoint for problem(θ) catches errors
Zygote.@adjoint (problem::ProblemLike)(θ) = safe_problem_adjoint(problem,θ)

# Ensure that the adjoint for problem(θ) catches errors and is real-valued
Zygote.@adjoint (problem::ProblemLike{SdeRootsEngine})(θ) = begin
  lp,Bc = safe_problem_adjoint(problem,θ)
  if isinf(lp) return lp,Bc end
  B = function(l̄)
      p̄,θ̄c = Bc(l̄)
      θ̄ = map(x->real.(x),θ̄c)
      return p̄,θ̄
    end
  return lp,B
end

transform_ss(problem::Problem,proc) = SdeProcess(sde(problem,proc.eq))

(problem::ProblemProb)(θ) = (problem.like)(θ)+(problem.prior)(θ)

"
Time interval for which sufficient data is available to estimate the autocorrelation.
"
function minimal_time(problem::ProblemLike,N=default_N_min)
  Δt = diff([x.t for x in problem.data])
  Tmin = minval_with_data(Δt,N)
  return Tmin
end

minimal_time(problem::ProblemProb) = minimal_time(problem.like)

function Statistics.mean(procs::Vector{<:SdeProcess},::Problem)
  eq = mean(getproperty.(procs,:eq))
  return SdeProcess(eq)
end


text:
#=
    Contents: Functions for server-side sessions that are common to all databases.
=#

"Read the session_id from the specified cookie."
function read_sessionid(req::Request, cookiename::AbstractString)
    session_id = ""
    if using_secure_cookies()
	session_id = get_securecookie_data(req, cookiename)
    else
	session_id = get_cookie_value(req, cookiename)
    end

    # Check whether session has timed out
    if !session_is_valid(session_id)
	session_id = ""
    end

    # Update lastreq
    if session_id != ""
	update_lastreq!(session_id)
    end
    session_id
end

read_sessionid(req::Request) = read_sessionid(req, config[:session][:cookiename])

"""
Returns true if:
    - session_id exists in the server-side data store, and
    - the session hasn't timed out

NOTE: If the session_id is not valid,  it is deleted from the database.
"""
function session_is_valid(session_id::AbstractString)
    result = session_is_valid(config[:session][:datastore], session_id)
    if result == false
	session_delete!(session_id)    # Delete session from database
    end
    result
end

function session_is_valid(cp::ConnectionPool, session_id::AbstractString)
    con    = get_connection!(cp)
    result = session_is_valid(con, session_id)
    release!(cp, con)
    result
end

"Updates session.lastreq."
function update_lastreq!(session_id::AbstractString)
    update_lastreq!(config[:session][:datastore], session_id)
end

function update_lastreq!(cp::ConnectionPool, session_id::AbstractString)
    con = get_connection!(cp)
    update_lastreq!(con, session_id)
    release!(cp, con)
end

"""
Init session_id => session on the database and set the specified cookie to session_id.
Return: session_id
"""
function session_create!(res::Response, username::AbstractString)
    cookiename = config[:session][:cookiename]
    session_create!(config[:session][:datastore], username, res, cookiename)
end

function session_create!(cp::ConnectionPool, username::AbstractString, res::Response, cookiename::AbstractString)
    con = get_connection!(cp)
    session_id = session_create!(con, username, res, cookiename)
    release!(cp, con)
    session_id
end

"Delete session from database and set the specified cookie to an invalid state."
function session_delete!(res::Response, session_id::AbstractString)
    session_delete!(res)           # Remove session ID cookie
    session_delete!(session_id)    # Delete session from database
end

"Delete session_id from the database."
function session_delete!(session_id::AbstractString)
    session_delete!(config[:session][:datastore], session_id)
end

function session_delete!(cp::ConnectionPool, session_id::AbstractString)
    con = get_connection!(cp)
    session_delete!(con, session_id)
    release!(cp, con)
end

function session_get(session_id, keys...)
    session_get(config[:session][:datastore], session_id, keys...)
end

function session_get(cp::ConnectionPool, session_id::AbstractString, keys...)
    con = get_connection!(cp)
    v   = session_get(con, session_id, keys...)
    release!(cp, con)
    v
end

function session_set!(session_id, keys_value...)
    k1 = keys_value[1]
    (k1 == "lastreq" || k1 == "username") && return    # Forbid changing username or lastreq
    session_set!(config[:session][:datastore], session_id, keys_value...)
end

function session_set!(cp::ConnectionPool, session_id::AbstractString, keys_value...)
    con = get_connection!(cp)
    session_set!(con, session_id, keys_value...)
    release!(cp, con)
end

# EOF

ext: jl


text:
#!/usr/bin/julia

# Copyright (C) 2021, 2022 Matheus Fernandes Bigolin <mfrdrbigolin@disroot.org>
# SPDX-License-Identifier: MIT

"Day Eight, Seven Segment Search."

#= This day involved some heavy logical deductions, so I attempted to provide
 = explananations in the program. Thus, I hope the reasoning I outlined
 = throught this code is clear to the reader.
 =#

include("Utils.jl")

using .Utils

"Count how many times the digits 1, 4, 7, or 8 appear in the `output` values."
solve1(output) = count.(d -> length(d) in [2,4,3,7], output) |> sum

"Segments used by each digit."
const DIGITS = ([
    [1,2,3,5,6,7], [3,6], [1,3,4,5,7], [1,3,4,6,7], [2,3,4,6],
    [1,2,4,6,7], [1,2,4:7...], [1,3,6], [1:7...], [1:4...,6,7]
])

"Decode the wires mapping in the `pattern` and return the configuration."
function decode_wires(pattern)
    conf = zeros(Int, 7)

    let_count = [sum(count.(in(letter), pattern)) for letter in 'a':'g']
    pat_length = length.(pattern)

    # Find the signal wires which have unique length.
    conf[2], conf[5], conf[6] = [findfirst(==(i), let_count) for i in [6,4,9]]

    # Segments which are useful for decoding the signal wires.
    one = Set(collect(pattern[findfirst(==(2), pat_length)]))
    seven = Set(collect(pattern[findfirst(==(3), pat_length)]))
    four = collect(pattern[findfirst(==(4), pat_length)]) .- ('a' - 1)

    # Which signal wire is in seven but not in one.
    conf[1] = first(setdiff(seven, one)) - 'a' + 1
    # The signal wire with length eight besides the top one.
    conf[3] = filter(!in(conf), findall(==(8), let_count))[1]
    # Signal wire in four not yet found.
    conf[4] = first(filter(!in(conf), four))
    # Find the only signal wire not yet marked.
    conf[7] = 28 - sum(conf)

    return conf
end

"Sum all the output values, given the `patterns` and the `digits`."
function solve2(patterns, digits)
    sum_output = 0

    for (i, pat) in enumerate(patterns)
        conf = decode_wires(pat)

        output = zeros(Int, 4)

        for (j, letter_digit) in enumerate(digits[i])
            digit = collect(letter_digit) .- ('a' - 1)
            conf_digit = sort([findfirst(==(d), conf) for d in digit])
            output[j] = findfirst(==(conf_digit), DIGITS) - 1
        end

        sum_output += array2num(output)
    end

    return sum_output
end

function main()
    usage_and_exit(length(ARGS) != 1)

    data = readlines(first(ARGS))

    entries = @. filter(!=("|"), split(data, " "))

    patterns = [getindex(entry, 1:10) for entry in entries]
    digits = [getindex(entry, 11:length(entries[1])) for entry in entries]

    println(solve1(digits))
    println(solve2(patterns, digits))
end

if abspath(PROGRAM_FILE) == @__FILE__
    main()
end


text:
module General

using ..TypeSynonyms

struct Contour
    center    ::RationalComplex
    radius    ::Float
    max_power ::Int
end
export Contour

abstract type ProblemType end
export ProblemType

abstract type Fiber{P <: ProblemType} end
export Fiber

abstract type Cohesive{P <: ProblemType} end
export Cohesive

abstract type SubMatrixType{P <: ProblemType} end
export SubMatrixType

struct SubMatrix_Zero{P <: ProblemType} <: SubMatrixType{P}
    contour ::Contour
end
export SubMatrix_Zero

struct SubMatrix_Fiber{P <: ProblemType} <: SubMatrixType{P}
    contour ::Contour
    fiber   ::Fiber{P}
end
export SubMatrix_Fiber

struct SubMatrix_Cohesive{P <: ProblemType} <: SubMatrixType{P}
    contour  ::Contour
    cohesive ::Cohesive{P}
end
export SubMatrix_Cohesive

abstract type SubMatrix_Averaging{P <: ProblemType} <: SubMatrixType{P} end
export SubMatrix_Averaging

struct ReferenceMatrix{P <: ProblemType}
    matrix ::Matrix{ SubMatrixType{P} }
end
export ReferenceMatrix

end # module General


text:
# adjusters.jl

module AdjustersTest

using Test
using Dates
using QDates

#trunc
@testset "trunc" begin
    qdt = QDates.QDate(2012,12,21)
    @test trunc(qdt, Dates.Year) == QDates.QDate(2012)
    @test trunc(qdt, Dates.Month) == QDates.QDate(2012,12)
    @test trunc(qdt, Dates.Day) == QDates.QDate(2012,12,21)
end

# Date functions
@testset "Date functions" begin
    # Name functions
    mtk = QDates.QDate(2017,1,1)
    ksrg = QDates.QDate(2017,2,2)
    yyi = QDates.QDate(2017,3,3)
    udk = QDates.QDate(2017,4,4)
    stk = QDates.QDate(2017,5,5)
    lstk = QDates.QDate(2017,5,true,6)
    mndk = QDates.QDate(2017,6,7)
    fdk = QDates.QDate(2017,7,8)
    hdk = QDates.QDate(2017,8,9)
    ngtk = QDates.QDate(2017,9,10)
    kmndk = QDates.QDate(2017,10,11)
    smtk = QDates.QDate(2017,11,12)
    sws = QDates.QDate(2017,12,13)

    @testset "firstdayofyear" begin
        @test QDates.firstdayofyear(sws) == QDates.QDate(2017,1,1)
        @test Dates.firstdayofyear(sws) == QDates.QDate(2017,1,1)
    end
    @testset "lastdayofyear" begin
        @test QDates.lastdayofyear(mtk) == QDates.QDate(2017,12,30)
        @test Dates.lastdayofyear(mtk) == QDates.QDate(2017,12,30)
    end
    @testset "lastdayofmonth" begin
        @test QDates.lastdayofmonth(mtk) == Dates.lastdayofmonth(mtk) == QDates.QDate(2017,1,29)
        @test QDates.lastdayofmonth(ksrg) == QDates.QDate(2017,2,30)
        @test QDates.lastdayofmonth(yyi) == QDates.QDate(2017,3,29)
        @test QDates.lastdayofmonth(udk) == QDates.QDate(2017,4,30)
        @test QDates.lastdayofmonth(stk) == QDates.QDate(2017,5,29)
        @test QDates.lastdayofmonth(lstk) == QDates.QDate(2017,5,true,29)
        @test QDates.lastdayofmonth(mndk) == QDates.QDate(2017,6,30)
        @test QDates.lastdayofmonth(fdk) == QDates.QDate(2017,7,29)
        @test QDates.lastdayofmonth(hdk) == QDates.QDate(2017,8,30)
        @test QDates.lastdayofmonth(ngtk) == QDates.QDate(2017,9,29)
        @test QDates.lastdayofmonth(kmndk) == QDates.QDate(2017,10,30)
        @test QDates.lastdayofmonth(smtk) == QDates.QDate(2017,11,30)
        @test QDates.lastdayofmonth(sws) == QDates.lastdayofyear(mtk) == QDates.QDate(2017,12,30)
    end
    @testset "firstdayofmonth" begin
        @test QDates.firstdayofmonth(mtk) == QDates.firstdayofyear(sws) == QDates.QDate(2017,1,1)
        @test QDates.firstdayofmonth(ksrg) == QDates.QDate(2017,2,1)
        @test QDates.firstdayofmonth(yyi) == QDates.QDate(2017,3,1)
        @test QDates.firstdayofmonth(udk) == QDates.QDate(2017,4,1)
        @test QDates.firstdayofmonth(stk) == QDates.QDate(2017,5,1)
        @test QDates.firstdayofmonth(lstk) == QDates.QDate(2017,5,true,1)
        @test QDates.firstdayofmonth(mndk) == QDates.QDate(2017,6,1)
        @test QDates.firstdayofmonth(fdk) == QDates.QDate(2017,7,1)
        @test QDates.firstdayofmonth(hdk) == QDates.QDate(2017,8,1)
        @test QDates.firstdayofmonth(ngtk) == QDates.QDate(2017,9,1)
        @test QDates.firstdayofmonth(kmndk) == QDates.QDate(2017,10,1)
        @test QDates.firstdayofmonth(smtk) == QDates.QDate(2017,11,1)
        @test QDates.firstdayofmonth(sws) == Dates.firstdayofmonth(sws) == QDates.QDate(2017,12,1)
    end
end

# # Adjusters
# # Adjuster Constructors # Not implemented
# @test QDates.QDate(Dates.isXXXX,2014) == QDates.QDate(2014,XX,XX)
# @test QDates.QDate(Dates.isYYYY,2014,5) == QDates.QDate(2014,5,YY)

# tonext, toprev, tofirst, tolast
@testset "tonext, toprev, tofirst, tolast" begin
    qdt = QDates.QDate(2017,5,3)
    @test Dates.tonext(qdt, QDates.先勝) == QDates.QDate(2017,5,9)
    @test Dates.tonext(qdt, QDates.先勝; same=true) == qdt
    @test Dates.tonext(qdt, QDates.友引) == QDates.QDate(2017,5,4)
    @test Dates.tonext(qdt, QDates.先負) == QDates.QDate(2017,5,5)
    @test Dates.tonext(qdt, QDates.仏滅) == QDates.QDate(2017,5,6)
    @test Dates.tonext(qdt, QDates.大安) == QDates.QDate(2017,5,7)
    @test Dates.tonext(qdt, QDates.赤口) == QDates.QDate(2017,5,8)
    # No dayofweek function for out of range values
    @test_throws KeyError Dates.tonext(qdt, 7)

    #test func, diff steps, negate, same
    @test Dates.tonext(QDates.is先勝, qdt) == QDates.QDate(2017,5,9)
    @test Dates.tonext(QDates.is先勝, qdt; same=true) == qdt
    @test Dates.tonext(QDates.is友引, qdt) == QDates.QDate(2017,5,4)
    @test Dates.tonext(QDates.is先負, qdt) == QDates.QDate(2017,5,5)
    @test Dates.tonext(QDates.is仏滅, qdt) == QDates.QDate(2017,5,6)
    @test Dates.tonext(QDates.is大安, qdt) == QDates.QDate(2017,5,7)
    @test Dates.tonext(QDates.is赤口, qdt) == QDates.QDate(2017,5,8)
    @test Dates.tonext(QDates.isleapmonth, qdt) == QDates.QDate(2017,5,true,1)
    @test Dates.tonext(QDates.isleapmonth, qdt; step=Dates.Month(1)) == QDates.QDate(2017,5,true,3)

    # @test Dates.tonext(x->!QDates.is大安(x), qdt; negate=true) == QDates.QDate(2017,5,7)
    @test Dates.tonext(!(x->!QDates.is大安(x)), qdt) == QDates.QDate(2017,5,7)
    # Reach adjust limit
    @test_throws ArgumentError Dates.tonext(QDates.is大安, qdt; limit=3)

    #toprev
    @test Dates.toprev(qdt, QDates.先勝) == QDates.QDate(2017,4,28)
    @test Dates.toprev(qdt, QDates.先勝; same=true) == qdt
    @test Dates.toprev(qdt, QDates.友引) == QDates.QDate(2017,4,29)
    @test Dates.toprev(qdt, QDates.先負) == QDates.QDate(2017,4,30)
    @test Dates.toprev(qdt, QDates.仏滅) == QDates.QDate(2017,4,25)
    @test Dates.toprev(qdt, QDates.大安) == QDates.QDate(2017,5,1)
    @test Dates.toprev(qdt, QDates.赤口) == QDates.QDate(2017,5,2)
    # No dayofweek function for out of range values
    @test_throws KeyError Dates.toprev(qdt, 8)

    #tofirst
    @test Dates.tofirst(qdt, QDates.先勝) == QDates.QDate(2017,5,3)
    @test Dates.tofirst(qdt, QDates.友引) == QDates.QDate(2017,5,4)
    @test Dates.tofirst(qdt, QDates.先負) == QDates.QDate(2017,5,5)
    @test Dates.tofirst(qdt, QDates.仏滅) == QDates.QDate(2017,5,6)
    @test Dates.tofirst(qdt, QDates.大安) == QDates.QDate(2017,5,1)
    @test Dates.tofirst(qdt, QDates.赤口) == QDates.QDate(2017,5,2)

    @test Dates.tofirst(qdt, QDates.先勝, of=Dates.Year) == QDates.QDate(2017,1,1)
    @test Dates.tofirst(qdt, QDates.友引, of=Dates.Year) == QDates.QDate(2017,1,2)
    @test Dates.tofirst(qdt, QDates.先負, of=Dates.Year) == QDates.QDate(2017,1,3)
    @test Dates.tofirst(qdt, QDates.仏滅, of=Dates.Year) == QDates.QDate(2017,1,4)
    @test Dates.tofirst(qdt, QDates.大安, of=Dates.Year) == QDates.QDate(2017,1,5)
    @test Dates.tofirst(qdt, QDates.赤口, of=Dates.Year) == QDates.QDate(2017,1,6)

    #tolast
    @test Dates.tolast(qdt, QDates.先勝) == QDates.QDate(2017,5,27)
    @test Dates.tolast(qdt, QDates.友引) == QDates.QDate(2017,5,28)
    @test Dates.tolast(qdt, QDates.先負) == QDates.QDate(2017,5,29)
    @test Dates.tolast(qdt, QDates.仏滅) == QDates.QDate(2017,5,24)
    @test Dates.tolast(qdt, QDates.大安) == QDates.QDate(2017,5,25)
    @test Dates.tolast(qdt, QDates.赤口) == QDates.QDate(2017,5,26)

    @test Dates.tolast(qdt, QDates.先勝, of=Dates.Year) == QDates.QDate(2017,12,26)
    @test Dates.tolast(qdt, QDates.友引, of=Dates.Year) == QDates.QDate(2017,12,27)
    @test Dates.tolast(qdt, QDates.先負, of=Dates.Year) == QDates.QDate(2017,12,28)
    @test Dates.tolast(qdt, QDates.仏滅, of=Dates.Year) == QDates.QDate(2017,12,29)
    @test Dates.tolast(qdt, QDates.大安, of=Dates.Year) == QDates.QDate(2017,12,30)
    @test Dates.tolast(qdt, QDates.赤口, of=Dates.Year) == QDates.QDate(2017,12,25)
end

end


text:
using Test
using PolyFields

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

add_jl(x) = endswith(x, ".jl") ? x : x*".jl"

if length(ARGS) > 0
    tests = map(add_jl, ARGS)
else
    tests = [
        "test_cell.jl",
        "test_species.jl",
        "test_interactions.jl",
        "test_system.jl",
        "test_updaters.jl"
    ]
end

println("Testing...")

for test in tests
    include(test)
end


text:
#! A simple Psi4 input script to compute a SCF reference using Psi4's libJK

psi4 = pyimport("psi4")
np   = pyimport("numpy") # used only to cast to Psi4 arrays

build_superfunctional = nothing
if VersionNumber(psi4.__version__) >= v"1.3a1"
    build_superfunctional = psi4.driver.dft.build_superfunctional
else
    build_superfunctional = psi4.driver.dft_funcs.build_superfunctional
end

function psi4view(psi4matrix)
   # Assumes Float64 type, C ordering
   if !hasproperty(psi4matrix,:__array_interface__)
      throw(ArgumentError("Input matrix cannot be accessed. Cannot be an rvalue"))
   end
   array_interface = psi4matrix.__array_interface__
   array_interface["data"][2] == false   || @warn "Not writable"
   array_interface["strides"] == nothing || @warn "Different ordering than C"
   array_interface["typestr"] == "<f8"   || @warn "Not little-endian Float64 eltype"
   ptr   = array_interface["data"][1]
   shape = reverse(array_interface["shape"])
   ndims = length(shape)
   permutedims(unsafe_wrap(Array{Float64,ndims}, Ptr{Float64}(ptr), shape), reverse(1:ndims))
end

# Diagonalize routine
function build_orbitals(diag, A, ndocc)
    Fp = psi4.core.triplet(A, diag, A, true, false, true)

    #A_view = psi4view(A)
    #nbf = size(A_view,1)
    nbf = A.shape[1]
    Cp = psi4.core.Matrix(nbf, nbf)
    eigvecs = psi4.core.Vector(nbf)
    Fp.diagonalize(Cp, eigvecs, psi4.core.DiagonalizeOrder.Ascending)

    C = psi4.core.doublet(A, Cp, false, false)
    C_jl = np.asarray(C)

    #Cocc = psi4.core.Matrix(nbf, ndocc)
    #Cocc_view = psi4view(Cocc) # returns C ordering
    #Cocc_view .= C_jl[:,1:ndocc] # or reverse because C ordering
    #Cocc.np[:] = C.np[:, 1:ndocc] # It is not going to work
    Cocc_jl  = zeros(nbf,ndocc)
    Cocc_jl .= C_jl[:,1:ndocc] # or reverse because C ordering
    Cocc = psi4.core.Matrix.from_array(Cocc_jl)

    D = psi4.core.doublet(Cocc, Cocc, false, true)

    C, Cocc, D, eigvecs
end

function ks_solver(alias, mol, options, V_builder::Function,
                   jk_type="DF", output="output.dat", restricted=true)

    # Build our molecule
    mol = mol.clone()
    mol.reset_point_group("c1")
    mol.fix_orientation(true)
    mol.fix_com(true)
    mol.update_geometry()

    # Set options
    psi4.set_output_file(output)

    psi4.core.prepare_options_for_module("SCF")
    psi4.set_options(options)
    psi4.core.set_global_option("SCF_TYPE", jk_type)

    maxiter = 20
    E_conv = psi4.core.get_option("SCF", "E_CONVERGENCE")
    D_conv = psi4.core.get_option("SCF", "D_CONVERGENCE")

    # Integral generation from Psi4's MintsHelper
    wfn = psi4.core.Wavefunction.build(mol, psi4.core.get_global_option("BASIS"))
    mints = psi4.core.MintsHelper(wfn.basisset())
    S = mints.ao_overlap()

    # Build the V Potential
    sup = build_superfunctional(alias, restricted)[1]
    sup.set_deriv(2)
    sup.allocate()

    vname = "RV"
    if !restricted
        vname = "UV"
    end
    Vpot = psi4.core.VBase.build(wfn.basisset(), sup, vname)
    Vpot.initialize()

    # Get nbf and ndocc for closed shell molecules
    nbf = wfn.nso()
    ndocc = wfn.nalpha()
    if wfn.nalpha() != wfn.nbeta()
        error("Only valid for RHF wavefunctions!")
    end

    println("\nNumber of occupied orbitals: ", ndocc)
    println("Number of basis functions:   ", nbf)

    # Build H_core
    V = mints.ao_potential()
    T = mints.ao_kinetic()
    H = T.clone()
    H.add(V)

    # Orthogonalizer A = S^(-1/2)
    A = mints.ao_overlap()
    A.power(-0.5, 1.e-14)

    # Build core orbitals
    C, Cocc, D, eigs = build_orbitals(H, A, ndocc)

    # Setup data for DIIS
    #t = @elapsed begin
       E = 0.0
       Enuc = mol.nuclear_repulsion_energy()
       Eold = 0.0

       # Initialize the JK object
       jk = psi4.core.JK.build(wfn.basisset())
       jk.set_memory(Int(1.25e8))  # 1GB
       jk.initialize()
       jk.print_header()

       diis_obj = psi4.p4util.solvers.DIIS(max_vec=3, removal_policy="largest")
    #end

    #println("\nTotal time taken for setup: $t seconds")

    println("\nStarting SCF iterations:")

    println("\n    Iter            Energy             XC E         Delta E        D RMS\n")
    SCF_E = 0.0
    #t = @elapsed begin

       for SCF_ITER in 1:maxiter

           # Compute JK
           jk.C_left_add(Cocc)
           jk.compute()
           jk.C_clear()

           # Build Fock matrix
           F = H.clone()
           F.axpy(2.0, jk.J()[1])
           F.axpy(-Vpot.functional().x_alpha(), jk.K()[1])

           # Build V
           ks_e = 0.0

           Vpot.set_D([D])
           Vpot.properties()[1].set_pointers(D)
           V = V_builder(D, Vpot)
           if isnothing(V)
               ks_e = 0.0
           else
               ks_e, V = V
               V = psi4.core.Matrix.from_array(V)
           end

           F.axpy(1.0, V)

           # DIIS error build and update
           diis_e = psi4.core.triplet(F, D, S, false, false, false)
           diis_e.subtract(psi4.core.triplet(S, D, F, false, false, false))
           diis_e = psi4.core.triplet(A, diis_e, A, false, false, false)

           diis_obj.add(F, diis_e)

           dRMS = diis_e.rms()

           # SCF energy and update
           SCF_E  = 2H.vector_dot(D)
           SCF_E += 2jk.J()[1].vector_dot(D)
           SCF_E -= Vpot.functional().x_alpha() * jk.K()[1].vector_dot(D)
           SCF_E += ks_e
           SCF_E += Enuc

           printfmt("SCF Iter{1:3d}: {2:18.14f}   {3:11.7f}   {4:1.5E}   {5:1.5E}\n",
                          SCF_ITER, SCF_E, ks_e, (SCF_E - Eold), dRMS)
           if abs(SCF_E - Eold) < E_conv && dRMS < D_conv
               break
           end

           Eold = SCF_E

           # DIIS extrapolate
           F = diis_obj.extrapolate()

           # Diagonalize Fock matrix
           C, Cocc, D, eigs = build_orbitals(F, A, ndocc)

           if SCF_ITER == maxiter
               error("Maximum number of SCF cycles exceeded.")
           end
       end
    #end

    #println("\nTotal time for SCF iterations: $t seconds.")

    printfmt("\nFinal SCF energy: {:.8f} hartree \n", SCF_E)

    data = Dict()
    data["Da"] = D
    data["Ca"] = C
    data["eigenvalues"] = eigs
    SCF_E, data
end


text:
#=
MIP verification of neural networks using MIPVerify.jl library.
Note: L1, L2 and Linf are supported.

Example on how to run:
`julia verify.jl "model.mat" "model eps=0.1/" mnist cnn_lenet_small inf 0.1 1000 120 5 lp`
- `model.mat`: the name of the .mat file containing model weights (note, mip_verify.py contains a
script on how to transform the weights from the original format to the one supported by MIPVerify.jl)
- `model eps=0.1/`: the name of the folder that will contain the output of the verifier
- `mnist`: the dataset
- `cnn_lenet_small`: the CNN architecture (the definition is given below)
- `inf`: the Lp-norm used to verify the existence of adversarial examples
- `0.1`: the radius of the Lp-norm
- `1000`: the number of points to verify
- `120`: the timeout of the verifier (in seconds)
- `5`: the timeout of the presolver (in seconds)
- `lp`: the type of the presolver (lp - linear programming, ia - interval arithmetics)
=#

using MIPVerify
using Gurobi
using Memento
using MAT

mat_mip_file = ARGS[1]
out_mip_file_dir = ARGS[2]
dataset_name = ARGS[3]
nn_type = ARGS[4]
norm_str = ARGS[5]
eps = parse(Float64, ARGS[6])
n_eval = parse(Int, ARGS[7])
time_limit = parse(Int, ARGS[8])
time_limit_bounds = parse(Int, ARGS[9])
presolver = ARGS[10]

if norm_str == "inf"
	norm = Inf
elseif norm_str == "2"
	norm = 2
else
	error("wrong norm_str")
end

if presolver == "ia"
	tightening_alg = interval_arithmetic
elseif presolver == "lp"
	tightening_alg = lp
else
	error("wrong presolver")
end

param_dict = mat_mip_file |> matread;
dataset = read_datasets(dataset_name)

n_h = size(dataset.test.images, 2)
n_w = size(dataset.test.images, 3)
n_c = size(dataset.test.images, 4)
n_out = length(unique(dataset.test.labels))
n_in = n_h * n_w * n_c  # needed only for FC networks
if nn_type == "fc1"
	n_hidden = 1024
	fc1 = get_matrix_params(param_dict, "fc1", (n_in, n_hidden))
	softmax = get_matrix_params(param_dict, "softmax", (n_hidden, n_out))

	nnparams = Sequential([
		Flatten(4),
		fc1, ReLU(interval_arithmetic),
		softmax
		],
		"$(out_mip_file_dir)"
	)
elseif nn_type == "cnn_lenet_small"
	conv1 = get_conv_params(param_dict, "conv1", (4, 4, n_c, 16), expected_stride = 2)

	conv2 = get_conv_params(param_dict, "conv2", (4, 4, 16, 32), expected_stride = 2)

	fc1 = get_matrix_params(param_dict, "fc1", (div(div(n_h, 2), 2)^2 * 32, 100))
	softmax = get_matrix_params(param_dict, "logits", (100, n_out))

	nnparams = Sequential(
		[   conv1, ReLU(interval_arithmetic),
		conv2, ReLU(),
# 		Flatten([1, 3, 2, 4]),  # this works for PyTorch standard
		Flatten([1, 4, 3, 2]),  # this works for TF standard
		fc1, ReLU(),
		softmax],
		"$(out_mip_file_dir)"
	)
else
	error("wrong model_name")
end

# Important to check whether we imported the model correctly
f = frac_correct(nnparams, dataset.test, 100)
println("Fraction correct of the first 100: $(f)")

# On which points to run MIP
target_indexes = 1:n_eval

MIPVerify.setloglevel!("info")  # "info" to get the number of unstable ReLUs

if norm == Inf
	MIPVerify.batch_find_untargeted_attack(
		nnparams,
		dataset.test,
		target_indexes,
		GurobiSolver(Gurobi.Env(), BestObjStop=eps, TimeLimit=time_limit),
		norm_order=Inf,
		tightening_algorithm=tightening_alg,
		rebuild=true,
		cache_model=false,
		solve_if_predicted_in_targeted=false,
		tightening_solver=GurobiSolver(Gurobi.Env(), TimeLimit=time_limit_bounds, OutputFlag=0),
	    pp = MIPVerify.LInfNormBoundedPerturbationFamily(eps),
		solve_rerun_option = MIPVerify.never
	)
elseif norm == 2
	MIPVerify.batch_find_untargeted_attack(
		nnparams,
		dataset.test,
		target_indexes,
		GurobiSolver(Gurobi.Env(), BestObjStop=eps, TimeLimit=time_limit),
		norm_order=2,
		tightening_algorithm=tightening_alg,
		rebuild=true,
		cache_model=false,
		solve_if_predicted_in_targeted=false,
		tightening_solver=GurobiSolver(Gurobi.Env(), TimeLimit=time_limit_bounds, OutputFlag=0),
	    pp = MIPVerify.L2NormBoundedPerturbationFamily(eps),
		solve_rerun_option = MIPVerify.never
	)
else
	error("wrong norm of perturbations")
end


text:
#class defintion
#type Person
#    bar
#    name::String
 #   aage::Int
#    other
#end
struct Foo
    bar
    baz
end
zhangsan = Foo("张三",32)
println(zhangsan)
typeof(zhangsan)
#println(zhangsan.name)

max_stars_repo_name: zqsheng/snippet

max_issues_repo_name: zqsheng/snippet

max_forks_repo_name: zqsheng/snippet


text:
using FTCTests
using Test
using Transducers
using LinearAlgebra
using DataFrames

@testset "evaluate" begin
    t0, tf = 0.0, 1.0
    ts = t0:0.01:tf
    d = 3
    sol = ts |> Map(t -> (; plant= (;state = (; p=ones(d))))) |> collect
    df = DataFrame(time=ts, sol=sol)
    traj_des(t) = zeros(d)
    jld2 = Dict("df" => df, "traj_des" => traj_des, "t0" => t0, "tf" => tf)
    t1 = 0.5
    error_norm = norm(ones(d))
    for threshold in [error_norm-0.1, error_norm+0.1]
        is_success = evaluate(jld2, t1, threshold)
        if error_norm < threshold
            @test is_success
        else
            @test !is_success
        end
    end
end


text:
using Base.Test
using QuantumOptics

particlenumber = 3
modenumber = 3
b_spin = SpinBasis(1//2)

m = complex(randn(modenumber, modenumber))
b = BosonicNParticleBasis(particlenumber, modenumber)
op = DenseOperator(GenericBasis([modenumber]), m)
op2 = SparseOperator(GenericBasis([modenumber]), sparse(m))
op_ = nparticleoperator_1(b, op)
op2_ = nparticleoperator_1(b, op2)
tic()
op_ = nparticleoperator_1(b, op)
toc()
tic()
op2_ = nparticleoperator_1(b, op2)
toc()

@test tracedistance(op_, full(op2_)) < 1e-12

m = complex(randn(modenumber^2, modenumber^2))
b = BosonicNParticleBasis(particlenumber, modenumber)
op = DenseOperator(GenericBasis([modenumber^2]), m)
op2 = SparseOperator(GenericBasis([modenumber^2]), sparse(m))
op_ = nparticleoperator_2(b, op)
op2_ = nparticleoperator_2(b, op2)
tic()
op_ = nparticleoperator_2(b, op)
toc()
tic()
op2_ = nparticleoperator_2(b, op2)
toc()

@test tracedistance(op_, full(op2_)) < 1e-12

# print(op_)
# print(op2_)

# println(length(b.occupations))
# println(length(b))


text:
@testset "Test basis" begin

using jFEMTools
using Tensors
import Tensors: Vec, gradient, norm
jF = jFEMTools

#Test Dubiner 3d basis
@show "Testing Dubiner 3D basis"
interpolation = jF.Dubiner{jF.Tetrahedron,3}()
exact_val = [0.8660254037844386,2.7386127875258306,0.8660254037844386,1.118033988749895,
            7.685213074469699,5.477225575051661,6.48074069840786,1.71846588560844,
            3.7416573867739413,-0.33071891388307384]
for i in 1:10
    @test jF.dubiner_basis(0.5,0.5,0.5,i,2) ≈ exact_val[i]
end

#Test Dubiner 2d base
interpolation = jF.Dubiner{jF.Triangle,4}()
quad_rule = jF.QuadratureRule{jF.Triangle}(jF.Strang(),5)

# Test nodal points for Lagrange base
nodal_points, topology = jF.get_nodal_points(jF.Triangle,1)
@test  nodal_points == [Vec{2, Float64}((0.0, 0.0)),
                        Vec{2, Float64}((1.0, 0.0)),
                        Vec{2, Float64}((0.0, 1.0))]
@test topology == Dict(0=>3,1=>0,2=>0)
nodal_points, topology = jF.get_nodal_points(jF.Triangle,2)
@test nodal_points == [Vec{2, Float64}((0.0, 0.0)),
                        Vec{2, Float64}((1.0, 0.0)),
                        Vec{2, Float64}((0.0, 1.0)),
                        Vec{2, Float64}((0.5, 0.5)),
                        Vec{2, Float64}((0.0, 0.5)),
                        Vec{2, Float64}((0.5, 0.0))]
@test topology == Dict(0=>3,1=>3,2=>0)
nodal_points, topology = jF.get_nodal_points(jF.Triangle, 3)
@test length(nodal_points) == 10
@test topology == Dict(0=>3,1=>6,2=>1)
@test nodal_points[10] ≈ Vec{2, Float64}((1/3, 1/3))
# Test Lagrange base
function ref_value(i::Int, ξ::Vec{2})
    ξ_x = ξ[1]
    ξ_y = ξ[2]
    i == 2 && return ξ_x
    i == 3 && return ξ_y
    i == 1 && return 1. - ξ_x - ξ_y
    throw(ArgumentError("no shape function $i for interpolation $ip"))
end

interpolation = jF.Lagrange(jF.Triangle,1)
ξ_ref = Vec{2}((0.5,0.5))
for j = 1:3
    @test abs(ref_value(j,ξ_ref) - jF.value(interpolation,j,ξ_ref)) < eps(Float64)
    @test norm(gradient(ξ -> ref_value(j, ξ), ξ_ref) - jF.gradient_value(interpolation,j,ξ_ref), Inf) < eps(Float64)
end

function ref_value1(i::Int, ξ::Vec{1})
    ξ_x = ξ[1]
    i == 1 && return 1-ξ_x
    i == 2 && return ξ_x
    throw(ArgumentError("no shape function $i for interpolation $ip"))
end
ξ_ref = Vec{1}((0.5,))
interpolation = jF.Lagrange(jF.Segment,1)
for j = 1:2
    @test abs(ref_value1(j,ξ_ref) - jF.value(interpolation,j,ξ_ref)) < eps(Float64)
    @test norm(gradient(ξ -> ref_value1(j, ξ), ξ_ref) - jF.gradient_value(interpolation,j,ξ_ref), Inf) < eps(Float64)
end

#Test Legendre basis
function ref_value(i::Int, x::Real)
    @assert i <= 3
    if i == 0
        return 1.0
    elseif i == 1
        return sqrt(3)*(2*x - 1)
    elseif i == 2
        return sqrt(5)*(6*x^2 - 6*x + 1)
    else
        return sqrt(7)*(2*x - 1)*(10*x^2 - 10*x + 1)
    end
end
function ref_dvalue(i::Int, x::Real)
    @assert i <= 3
    if i == 0
        return 0
    elseif i == 1
        return sqrt(3)*2
    elseif i == 2
        return sqrt(5)*(12*x - 6)
    else
        return sqrt(7)*12*(5*x^2 - 5*x + 1)
    end
end

# interpolation = jF.Legendre(jF.Segment,3)
# quad_rule = QuadratureRule{1,RefTetrahedron}(GaussLegendre(),3)
# for i = 0:3
#     @test integrate(x->(ref_value(i,x[1]) - value(interpolation, i+1, x))^2,quad_rule) < eps()
#     @test integrate(x->(ref_dvalue(i,x[1]) - gradient_value(interpolation, i+1, x)[1])^2,quad_rule) < eps()
# end

end

max_stars_repo_name: Paulms/jFEMToools

max_issues_repo_name: Paulms/jFEMToools

max_forks_repo_name: Paulms/jFEMToools


text:
#=
  Hidden Markov Models

  This is from the BLOG distribution hmm.dblog
  """
  A hidden Markov model with four states and four output symbols.
  The parameters of this model are completely made-up.
  DBLOG model
  """

A:1 C:2 G:3 T:4

Distributions of variable S1 (num:0)
2.00000 =>    7524  (0.752400)
4.00000 =>     853  (0.085300)
3.00000 =>     813  (0.081300)
1.00000 =>     810  (0.081000)

Distributions of variable S2 (num:0)
1.00000 =>    9778  (0.977800)
2.00000 =>     176  (0.017600)
3.00000 =>      40  (0.004000)
4.00000 =>       6  (0.000600)

Distributions of variable S3 (num:0)
1.00000 =>    9226  (0.922600)
4.00000 =>     640  (0.064000)
2.00000 =>      84  (0.008400)
3.00000 =>      50  (0.005000)

Distributions of variable S4 (num:0)
1.00000 =>    8574  (0.857400)
4.00000 =>    1006  (0.100600)
2.00000 =>     300  (0.030000)
3.00000 =>     120  (0.012000)

Distributions of variable S5 (num:0)
3.00000 =>    7383  (0.738300)
4.00000 =>    1365  (0.136500)
1.00000 =>     724  (0.072400)
2.00000 =>     528  (0.052800)

Distributions of variable S6 (num:0)
4.00000 =>    3381  (0.338100)
2.00000 =>    2506  (0.250600)
1.00000 =>    2320  (0.232000)
3.00000 =>    1793  (0.179300)

Distributions of variable S7 (num:0)
1.00000 =>    2828  (0.282800)
4.00000 =>    2816  (0.281600)
2.00000 =>    2467  (0.246700)
3.00000 =>    1889  (0.188900)

Distributions of variable S8 (num:0)
2.00000 =>    3246  (0.324600)
3.00000 =>    2758  (0.275800)
4.00000 =>    2169  (0.216900)
1.00000 =>    1827  (0.182700)

  Cf ~/webppl/hmm1.wppl

=#

using Turing, StatsPlots, DataFrames
include("jl_utils.jl")

@model function hmm1()
    n = 8

    A = 1
    C = 2
    G = 3
    T = 4
    states = [A,C,G,T]

    AA ~ Uniform(0,1)
    AC ~ Uniform(0,1)
    AG ~ Uniform(0,1)
    AT ~ Uniform(0,1)

    CA ~ Uniform(0,1)
    CC ~ Uniform(0,1)
    CG ~ Uniform(0,1)
    CT ~ Uniform(0,1)

    GA ~ Uniform(0,1)
    GC ~ Uniform(0,1)
    GG ~ Uniform(0,1)
    GT ~ Uniform(0,1)

    TA ~ Uniform(0,1)
    TC ~ Uniform(0,1)
    TG ~ Uniform(0,1)
    TT ~ Uniform(0,1)

    S = tzeros(n)
    for t in 1:n
        if t == 1
            S[t] ~ Categorical([0.3,0.2,0.1,0.4])
        else
            prev_t = S[t-1]
            S[t] ~
                (prev_t == A) ? Categorical(simplex([AA,AC,AG,AT])) :
                (prev_t == C) ? Categorical(simplex([CA,CC,CG,CT])) :
                (prev_t == G) ? Categorical(simplex([GA,GC,GG,GT])) :
                (prev_t == T) ? Categorical(simplex([TA,TC,TG,TT])) : Dirac(5)
        end
    end

    ResultA = 1
    ResultC = 2
    ResultG = 3
    ResultT = 4
    O = tzeros(n)
    for t in 1:n
        St = S[t]
        O[t] ~ (St == A) ? Categorical([0.85,0.05,0.05, 0.05]) : # [ResultA,ResultC,ResultG,ResultT]
               (St == C) ? Categorical([0.05,0.85,0.05, 0.05]) :
               (St == G) ? Categorical([0.05,0.05,0.85, 0.05]) :
               (St == T) ? Categorical([0.05,0.05,0.05, 0.85]) : Dirac(5)
    end

    # Evidence for the Hidden Markov Model.
    true ~ Dirac(O[1] == ResultC)
    true ~ Dirac(O[2] == ResultA)
    true ~ Dirac(O[3] == ResultA)
    true ~ Dirac(O[4] == ResultA)
    true ~ Dirac(O[5] == ResultG)

    # """
    # Queries for the Hiddem Markov Model, given the evidence.
    # Note that we can query S(5) even though our observations only
    # went up to time 4.
    # """
    # Note: I changed from BLOG/WebPPL's 0-based to 1-based.
    S1 ~ Dirac(S[1])
    S2 ~ Dirac(S[2])
    S3 ~ Dirac(S[3])
    S4 ~ Dirac(S[4])
    S5 ~ Dirac(S[5])
    S6 ~ Dirac(S[6])
    S7 ~ Dirac(S[7])
    S8 ~ Dirac(S[8])

end

model = hmm1()

num_chains = 4

# chns = sample(model, Prior(), 10_000)
# chns = sample(model, MH(), 10_000)
chns = sample(model, PG(15), 10_000)
# chns = sample(model, SMC(), 10_000)
# chns = sample(model, IS(), 10_000)

display(chns)
# display(plot(chns))

println("A:1 C:2 G:3 T:4")
show_var_dist_pct(chns,:S1)
show_var_dist_pct(chns,:S2)
show_var_dist_pct(chns,:S3)
show_var_dist_pct(chns,:S4)
show_var_dist_pct(chns,:S5)
show_var_dist_pct(chns,:S6)
show_var_dist_pct(chns,:S7)
show_var_dist_pct(chns,:S8)

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
#--define directory of .l and .html files
cd(raw"C:\julia-playgroung\calculadora")

#--teste de Aplicação web
using Blink;println("using Blink")
Blink.AtomShell.install();
#varinfo(Blink)#whos(Blink)

#--publish window with a massage
w = Window();
f = open("test.html") do file
            read(file,String)
    end

#publish window w with html file f
body!(w,f)

max_stars_repo_name: lugabiel/calculadora

max_issues_repo_name: lugabiel/calculadora

max_forks_repo_name: lugabiel/calculadora


text:
module Lp
	using Shared_funcs
	using Dcd
	using Conv

	:Created
	:Infeasible
	:Optimal
	:Unbounded
	:Maxit
	:Error

	:Cf0
	:Cf2b

	type SingleVar_translation{T}
		index::Int
		op1_mul::T
		op2_add::T
		SingleVar_translation() = (x = new(); index = -1; op1_mul = one(T); op2_add = zero(T); return x;)
	end

	type Cf0_problem{T} #Minimize
		conv::Conv.Converter
		n::Int
		m::Int
		c::Vector{T}
		A::Matrix{T}
		b::Vector{T}
		params::Params
		z_transl::SingleVar_translation{T}

		Cf0_problem() = new()
	end

	type Cf2b_problem{T} #Minimize
		cf0::Cf0_problem{T}
		l::Vector{T}
		h::Vector{T}

		Cf2b_problem() = new()
	end

	type Solution{T}
		solved::Bool
		status::Symbol
		z::T
		x::Vector{T}
		iter::Int

		dcd::Dcd.Session

		function Solution(params::Params)
			x = new()
			x.solved = false
			x.status = :Created
			x.dcd = Dcd.Session(get(params, "Dcd", false))
			x.iter = 0
			return x
		end
	end

	function get_form(prob::Cf0_problem) return :CF0 end
	function get_form(prob::Cf2b_problem) return :CF2b end
	function get_t(prob::Cf0_problem) return prob.conv.t end
	function get_t(prob::Cf2b_problem) return get_t(prob.cf0) end
	function get_n(prob::Cf0_problem) return prob.n end
	function get_m(prob::Cf0_problem) return prob.m end
	function get_n(prob::Cf2b_problem) return get_n(prob.cf0) end
	function get_m(prob::Cf2b_problem) return get_m(prob.cf0) end

	function comp_density(prob::Cf0_problem) return Shared_funcs.comp_density(prob.A[1:prob.m, 1:prob.n]) end
	function comp_density(prob::Cf2b_problem) return comp_density(prob.cf0) end

	function construct_solution(type_t::DataType, params::Params)
		return Solution{type_t}(params)
	end

	function transl_single{T}(transl::SingleVar_translation{T}, v::T)
		return (v*transl.op1_mul)+transl.op2_add
	end

	function translate_solution{T}(prob::Cf0_problem{T}, raw_sol::Solution{T})
		sol = construct_solution(prob.conv.t, prob.params)
		sol.solved = raw_sol.solved
		sol.status = raw_sol.status
		sol.dcd = raw_sol.dcd
		sol.iter = raw_sol.iter
		if (sol.solved)
			sol.z = transl_single(prob.z_transl, raw_sol.z)
			sol.x = raw_sol.x
		else
			sol.z = zero(T)
			sol.x = Array(T, 0)
		end

		return sol
	end

	function translate_solution{T}(prob::Cf2b_problem{T}, raw_sol::Solution{T})
		return translate_solution(prob.cf0, raw_sol)
	end

	function construct_Cf0_problem(params::Params)
		conv = Conv.converter(params["type"])
		prob = Cf0_problem{conv.t}()
		prob.conv = conv
		prob.z_transl = SingleVar_translation{prob.conv.t}()
		prob.params = params
		return prob
	end

	function construct_Cf2b_problem(params::Params)
		cf0 = construct_Cf0_problem(params)
		prob = Cf2b_problem{cf0.conv.t}()
		prob.cf0 = cf0
		return prob
	end

	# min: z = cx, subj: Ax <= b, x >= 0
	function fill_problem{T}(params::Params, prob::Cf0_problem{T}, c::Vector, A::Matrix, b::Vector)
		prob.n = length(c)
		prob.m = length(b)
		prob.c = vcat( Conv.vector(prob.conv, c), zeros(T, prob.m) )
		prob.A = hcat( Conv.matrix(prob.conv, A), Conv.matrix(prob.conv, eye(prob.m)) )
		prob.b = Conv.vector(prob.conv, b)
	end

	# min: z = cx, subj: Ax <= b, x >= 0
	function create_min_Cf0_problem(params::Params, c::Vector, A::Matrix, b::Vector)
		prob = construct_Cf0_problem(params)
		fill_problem(params, prob, c, A, b)
		return prob
	end

	# max: z = cx, subj: Ax <= b, x >= 0
	function create_max_Cf0_problem(params::Params, c::Vector, A::Matrix, b::Vector)
		mone = Conv.conv(params["type"], -1)
		prob = create_min_Cf0_problem(params, mone*(c), A, b)
		prob.z_transl.op1_mul = mone
		return prob
	end

	# min: z = cx, subj: zl <= Ax <= zh, xl <= x <= xh
	function fill_problem{T}(params::Params, prob::Cf2b_problem{T}, c::Vector, A::Matrix,
			zl::Vector, zh::Vector, xl::Vector, xh::Vector)
		fill_problem(params, prob.cf0, c, A, Array(T, 0))
		prob.l = vcat( Conv.vector(prob.cf0.conv, xl), Conv.vector(prob.cf0.conv, zl) )
		prob.h = vcat( Conv.vector(prob.cf0.conv, xh), Conv.vector(prob.cf0.conv, zh) )
	end

	# min: z = cx, subj: zl <= Ax <= zh, xl <= x <= xh
	function create_min_Cf2b_problem(params::Params, c::Vector, A::Matrix,
			zl::Vector, zh::Vector, xl::Vector, xh::Vector)
		prob = construct_Cf2b_problem(params)
		fill_problem(params, prob, c, A, zlh, xlh)
		return prob
	end

	# max: z = cx, subj: zl <= Ax <= zh, xl <= x <= xh
	function create_max_Cf2b_problem(params::Params, c::Vector, A::Matrix,
			zl::Vector, zh::Vector, xl::Vector, xh::Vector)
		mone = Conv.conv(params["type"], -1)
		prob = create_min_Cf2b_problem(params, mone*(c), A, zlh, xlh)
		prob.cf0.z_transl.op1_mul = mone
		return prob
	end

	function dcd_var_info(n, i) return i <= n ? ("x", i) : ("w", i-n); end
	function dcd_var(n, i) s,i = dcd_var_info(n, i); return @sprintf "%s%d" s i; end
	function dcd_iR(n, m, iB)
		iR = Int[]
		for i = 1:m+n
			if findfirst(iB, i) == 0 push!(iR, i) end
		end
		return iR
	end

	function dcd_pivots_impl(sol, typ)

		function dcd_print_basis(n, iB)
			for i = 1:length(iB)
				vs,vi = dcd_var_info(n, iB[i])
				print(vs, vi)
			end
			@printf "\n"
		end

		function dcd_print_pivot_type(n, iB, iR, typ)
			if (typ == 2) dcd_print_basis(n, iB)
			elseif (typ == 3) dcd_print_basis(n, iR)
			elseif (typ == 4) println(sort(iB))
			elseif (typ == 5) println(sort(iR))
			end
		end

		function dcd_print_pivot(i, iter, n, iB, iR, typ)
			if (haskey(iter, "pivot"))
				c,r = iter["pivot"]
				if (typ == 1) println( i, ". ",  dcd_var(n, iR[c]), ",", dcd_var(n, iB[r]) ) end
				iR[c], iB[r] = iB[r], iR[c]
				if (typ != 1) dcd_print_pivot_type(n, iB, iR, typ) end
			end
		end

		if isempty(sol.dcd.iters) return; end
		n = sol.prob.n
		iB = deepcopy(sol.dcd.iters[1]["iB"])
		iR = dcd_iR(sol.prob.n, sol.prob.m, iB)
		if (typ != 1) dcd_print_pivot_type(n, iB, iR, typ) end
		for i = 1:length(sol.dcd.iters)
			iter = sol.dcd.iters[i]
			dcd_print_pivot(i, iter, n, iB, iR, typ)
		end
	end

	function dcd_pivots(sol) dcd_pivots_impl(sol, 1) end
	function dcd_basis(sol) dcd_pivots_impl(sol, 2) end
	function dcd_nbasis(sol) dcd_pivots_impl(sol, 3) end
	function dcd_ibasis(sol) dcd_pivots_impl(sol, 4) end
	function dcd_inbasis(sol) dcd_pivots_impl(sol, 5) end
	function dcd_iters(sol) println(length(sol.dcd.iters)) end
	function dcd_key(sol, key)
		for i = 1:length(sol.dcd.iters)
			iter = sol.dcd.iters[i]
			if (haskey(iter, key))
				val = iter[key]
				println(i, ". ", val)
			end
		end
	end
end

max_stars_repo_name: jadnohra/mlcp-lib

max_issues_repo_name: jadnohra/mlcp-lib

max_forks_repo_name: jadnohra/mlcp-lib


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

import Base.LinAlg, Base.LinAlg.BlasReal, Base.LinAlg.BlasComplex

srand(100)
## BLAS tests - testing the interface code to BLAS routines
@testset for elty in [Float32, Float64, Complex64, Complex128]
    @testset "syr2k!" begin
        U = randn(5,2)
        V = randn(5,2)
        if elty == Complex64 || elty == Complex128
            U = complex.(U, U)
            V = complex.(V, V)
        end
        U = convert(Array{elty, 2}, U)
        V = convert(Array{elty, 2}, V)
        @test tril(LinAlg.BLAS.syr2k('L','N',U,V)) ≈ tril(U*V.' + V*U.')
        @test triu(LinAlg.BLAS.syr2k('U','N',U,V)) ≈ triu(U*V.' + V*U.')
        @test tril(LinAlg.BLAS.syr2k('L','T',U,V)) ≈ tril(U.'*V + V.'*U)
        @test triu(LinAlg.BLAS.syr2k('U','T',U,V)) ≈ triu(U.'*V + V.'*U)
    end

    if elty in (Complex64, Complex128)
        @testset "her2k!" begin
            U = randn(5,2)
            V = randn(5,2)
            U = complex.(U, U)
            V = complex.(V, V)
            U = convert(Array{elty, 2}, U)
            V = convert(Array{elty, 2}, V)
            @test tril(LinAlg.BLAS.her2k('L','N',U,V)) ≈ tril(U*V' + V*U')
            @test triu(LinAlg.BLAS.her2k('U','N',U,V)) ≈ triu(U*V' + V*U')
            @test tril(LinAlg.BLAS.her2k('L','C',U,V)) ≈ tril(U'*V + V'*U)
            @test triu(LinAlg.BLAS.her2k('U','C',U,V)) ≈ triu(U'*V + V'*U)
        end
    end

    o4 = ones(elty, 4)
    z4 = zeros(elty, 4)

    I4 = eye(elty, 4)
    I43 = eye(elty, 4, 3)
    L4 = tril(ones(elty, (4,4)))
    U4 = triu(ones(elty, (4,4)))
    Z4 = zeros(elty, (4,4))

    elm1 = convert(elty, -1)
    el2 = convert(elty, 2)
    v14 = convert(Vector{elty}, [1:4;])
    v41 = convert(Vector{elty}, [4:-1:1;])

    let n = 10
        @testset "dot products" begin
            if elty <: Real
                x1 = convert(Vector{elty}, randn(n))
                x2 = convert(Vector{elty}, randn(n))
                @test BLAS.dot(x1,x2) ≈ sum(x1.*x2)
                @test_throws DimensionMismatch BLAS.dot(x1,rand(elty, n + 1))
            else
                z1 = convert(Vector{elty}, complex.(randn(n),randn(n)))
                z2 = convert(Vector{elty}, complex.(randn(n),randn(n)))
                @test BLAS.dotc(z1,z2) ≈ sum(conj(z1).*z2)
                @test BLAS.dotu(z1,z2) ≈ sum(z1.*z2)
                @test_throws DimensionMismatch BLAS.dotc(z1,rand(elty, n + 1))
                @test_throws DimensionMismatch BLAS.dotu(z1,rand(elty, n + 1))
            end
        end
        @testset "iamax" begin
            if elty <: Real
                x = convert(Vector{elty}, randn(n))
                @test BLAS.iamax(x) == indmax(abs.(x))
            else
                z = convert(Vector{elty}, complex.(randn(n),randn(n)))
                @test BLAS.iamax(z) == indmax(map(x -> abs(real(x)) + abs(imag(x)), z))
            end
        end
        @testset "axpy" begin
            if elty <: Real
                x1 = convert(Vector{elty}, randn(n))
                x2 = convert(Vector{elty}, randn(n))
                α  = rand(elty)
                @test BLAS.axpy!(α,copy(x1),copy(x2)) ≈ x2 + α*x1
                @test_throws DimensionMismatch BLAS.axpy!(α, copy(x1), rand(elty, n + 1))
                @test_throws DimensionMismatch BLAS.axpy!(α, copy(x1), 1:div(n,2), copy(x2), 1:n)
                @test_throws ArgumentError BLAS.axpy!(α, copy(x1), 0:div(n,2), copy(x2), 1:(div(n, 2) + 1))
                @test_throws ArgumentError BLAS.axpy!(α, copy(x1), 1:div(n,2), copy(x2), 0:(div(n, 2) - 1))
                @test BLAS.axpy!(α,copy(x1),1:n,copy(x2),1:n) ≈ x2 + α*x1
            else
                z1 = convert(Vector{elty}, complex.(randn(n), randn(n)))
                z2 = convert(Vector{elty}, complex.(randn(n), randn(n)))
                α  = rand(elty)
                @test BLAS.axpy!(α, copy(z1), copy(z2)) ≈ z2 + α * z1
                @test_throws DimensionMismatch BLAS.axpy!(α, copy(z1), rand(elty, n + 1))
                @test_throws DimensionMismatch BLAS.axpy!(α, copy(z1), 1:div(n, 2), copy(z2), 1:(div(n, 2) + 1))
                @test_throws ArgumentError BLAS.axpy!(α, copy(z1), 0:div(n,2), copy(z2), 1:(div(n, 2) + 1))
                @test_throws ArgumentError BLAS.axpy!(α, copy(z1), 1:div(n,2), copy(z2), 0:(div(n, 2) - 1))
                @test BLAS.axpy!(α,copy(z1),1:n,copy(z2),1:n) ≈ z2 + α*z1
            end
        end
        @testset "nrm2, iamax, and asum for StridedVectors" begin
            a = rand(elty,n)
            b = view(a,2:2:n,1)
            @test BLAS.nrm2(b) ≈ norm(b)
            if elty <: Real
                @test BLAS.asum(b) ≈ sum(abs.(b))
                @test BLAS.iamax(b) ≈ indmax(abs.(b))
            else
                @test BLAS.asum(b) ≈ sum(abs.(real(b))) + sum(abs.(imag(b)))
                @test BLAS.iamax(b) == indmax(map(x -> abs(real(x)) + abs(imag(x)), b))
            end
        end
        # scal
        α = rand(elty)
        a = rand(elty,n)
        @test BLAS.scal(n,α,a,1) ≈ α * a

        @testset "trsv" begin
            A = triu(rand(elty,n,n))
            @testset "Vector and SubVector" for x in (rand(elty, n), view(rand(elty,2n),1:2:2n))
                @test A\x ≈ BLAS.trsv('U','N','N',A,x)
                @test_throws DimensionMismatch BLAS.trsv('U','N','N',A,ones(elty,n+1))
            end
        end
        @testset "ger, her, syr" for x in (rand(elty, n), view(rand(elty,2n), 1:2:2n)),
            y in (rand(elty,n), view(rand(elty,3n), 1:3:3n))

            A = rand(elty,n,n)
            α = rand(elty)

            @test BLAS.ger!(α,x,y,copy(A)) ≈ A + α*x*y'
            @test_throws DimensionMismatch BLAS.ger!(α,ones(elty,n+1),y,copy(A))

            A = rand(elty,n,n)
            A = A + A.'
            @test issymmetric(A)
            @test triu(BLAS.syr!('U',α,x,copy(A))) ≈ triu(A + α*x*x.')
            @test_throws DimensionMismatch BLAS.syr!('U',α,ones(elty,n+1),copy(A))

            if elty <: Complex
                A = rand(elty,n,n)
                A = A + A'
                α = real(α)
                @test triu(BLAS.her!('U',α,x,copy(A))) ≈ triu(A + α*x*x')
                @test_throws DimensionMismatch BLAS.her!('U',α,ones(elty,n+1),copy(A))
            end
        end
        @testset "copy" begin
            x1 = convert(Vector{elty}, randn(n))
            x2 = convert(Vector{elty}, randn(n))
            BLAS.copy!(x2, 1:n, x1, 1:n)
            @test x2 == x1
            @test_throws DimensionMismatch BLAS.copy!(x2, 1:n, x1, 1:(n - 1))
            @test_throws ArgumentError BLAS.copy!(x1, 0:div(n, 2), x2, 1:(div(n, 2) + 1))
            @test_throws ArgumentError BLAS.copy!(x1, 1:(div(n, 2) + 1), x2, 0:div(n, 2))
        end
        # trmv
        A = triu(rand(elty,n,n))
        x = rand(elty,n)
        @test BLAS.trmv('U','N','N',A,x) ≈ A*x
        @testset "symmetric/Hermitian multiplication" begin
            x = rand(elty,n)
            A = rand(elty,n,n)
            Aherm = A + A'
            Asymm = A + A.'
            @testset "symv and hemv" begin
                @test BLAS.symv('U',Asymm,x) ≈ Asymm*x
                @test_throws DimensionMismatch BLAS.symv!('U',one(elty),Asymm,x,one(elty),ones(elty,n+1))
                @test_throws DimensionMismatch BLAS.symv('U',ones(elty,n,n+1),x)
                if elty <: BlasComplex
                    @test BLAS.hemv('U',Aherm,x) ≈ Aherm*x
                    @test_throws DimensionMismatch BLAS.hemv('U',ones(elty,n,n+1),x)
                    @test_throws DimensionMismatch BLAS.hemv!('U',one(elty),Aherm,x,one(elty),ones(elty,n+1))
                end
            end

            @testset "symm error throwing" begin
                @test_throws DimensionMismatch BLAS.symm('L','U',ones(elty,n,n-1),rand(elty,n,n))
                @test_throws DimensionMismatch BLAS.symm('R','U',ones(elty,n-1,n),rand(elty,n,n))
                @test_throws DimensionMismatch BLAS.symm!('L','U',one(elty),Asymm,ones(elty,n,n),one(elty),rand(elty,n-1,n))
                @test_throws DimensionMismatch BLAS.symm!('L','U',one(elty),Asymm,ones(elty,n,n),one(elty),rand(elty,n,n-1))
                if elty <: BlasComplex
                    @test_throws DimensionMismatch BLAS.hemm('L','U',ones(elty,n,n-1),rand(elty,n,n))
                    @test_throws DimensionMismatch BLAS.hemm('R','U',ones(elty,n-1,n),rand(elty,n,n))
                    @test_throws DimensionMismatch BLAS.hemm!('L','U',one(elty),Aherm,ones(elty,n,n),one(elty),rand(elty,n-1,n))
                    @test_throws DimensionMismatch BLAS.hemm!('L','U',one(elty),Aherm,ones(elty,n,n),one(elty),rand(elty,n,n-1))
                end
            end
        end
        @testset "trmm error throwing" begin
            @test_throws DimensionMismatch BLAS.trmm('L','U','N','N',one(elty),triu(rand(elty,n,n)),ones(elty,n+1,n))
            @test_throws DimensionMismatch BLAS.trmm('R','U','N','N',one(elty),triu(rand(elty,n,n)),ones(elty,n,n+1))
        end

        #trsm
        A = triu(rand(elty,n,n))
        B = rand(elty,(n,n))
        @test BLAS.trsm('L','U','N','N',one(elty),A,B) ≈ A\B

        #will work for SymTridiagonal,Tridiagonal,Bidiagonal!
        @testset "banded matrix mv" begin
            @testset "gbmv" begin
                TD  = Tridiagonal(rand(elty,n-1),rand(elty,n),rand(elty,n-1))
                x   = rand(elty,n)
                #put TD into the BLAS format!
                fTD = zeros(elty,3,n)
                fTD[1,2:n] = TD.du
                fTD[2,:] = TD.d
                fTD[3,1:n-1] = TD.dl
                @test BLAS.gbmv('N',n,1,1,fTD,x) ≈ TD*x
            end
            #will work for SymTridiagonal only!
            @testset "sbmv" begin
                if elty <: BlasReal
                    ST  = SymTridiagonal(rand(elty,n),rand(elty,n-1))
                    x   = rand(elty,n)
                    #put TD into the BLAS format!
                    fST = zeros(elty,2,n)
                    fST[1,2:n] = ST.ev
                    fST[2,:] = ST.dv
                    @test BLAS.sbmv('U',1,fST,x) ≈ ST*x
                else
                    dv = real(rand(elty,n))
                    ev = rand(elty,n-1)
                    bH = zeros(elty,2,n)
                    bH[1,2:n] = ev
                    bH[2,:] = dv
                    fullH = diagm(dv) + diagm(conj(ev),-1) + diagm(ev,1)
                    @test BLAS.hbmv('U',1,bH,x) ≈ fullH*x
                end
            end
        end
    end

    @testset "gemv" begin
        @test all(BLAS.gemv('N', I4, o4) .== o4)
        @test all(BLAS.gemv('T', I4, o4) .== o4)
        @test all(BLAS.gemv('N', el2, I4, o4) .== el2 * o4)
        @test all(BLAS.gemv('T', el2, I4, o4) .== el2 * o4)
        @test_throws DimensionMismatch BLAS.gemv('N',I43,o4)
        o4cp = copy(o4)
        @test_throws DimensionMismatch BLAS.gemv!('T',one(elty),I43,o4,one(elty),o4cp)
        @test_throws DimensionMismatch BLAS.gemv!('C',one(elty),I43,o4,one(elty),o4cp)
        @test all(BLAS.gemv!('N', one(elty), I4, o4, elm1, o4cp) .== z4)
        @test all(o4cp .== z4)
        o4cp[:] = o4
        @test all(BLAS.gemv!('T', one(elty), I4, o4, elm1, o4cp) .== z4)
        @test all(o4cp .== z4)
        @test all(BLAS.gemv('N', U4, o4) .== v41)
        @test all(BLAS.gemv('N', U4, o4) .== v41)
    end
    @testset "gemm" begin
        @test all(BLAS.gemm('N', 'N', I4, I4) .== I4)
        @test all(BLAS.gemm('N', 'T', I4, I4) .== I4)
        @test all(BLAS.gemm('T', 'N', I4, I4) .== I4)
        @test all(BLAS.gemm('T', 'T', I4, I4) .== I4)
        @test all(BLAS.gemm('N', 'N', el2, I4, I4) .== el2 * I4)
        @test all(BLAS.gemm('N', 'T', el2, I4, I4) .== el2 * I4)
        @test all(BLAS.gemm('T', 'N', el2, I4, I4) .== el2 * I4)
        @test all(LinAlg.BLAS.gemm('T', 'T', el2, I4, I4) .== el2 * I4)
        I4cp = copy(I4)
        @test all(BLAS.gemm!('N', 'N', one(elty), I4, I4, elm1, I4cp) .== Z4)
        @test all(I4cp .== Z4)
        I4cp[:] = I4
        @test all(BLAS.gemm!('N', 'T', one(elty), I4, I4, elm1, I4cp) .== Z4)
        @test all(I4cp .== Z4)
        I4cp[:] = I4
        @test all(BLAS.gemm!('T', 'N', one(elty), I4, I4, elm1, I4cp) .== Z4)
        @test all(I4cp .== Z4)
        I4cp[:] = I4
        @test all(BLAS.gemm!('T', 'T', one(elty), I4, I4, elm1, I4cp) .== Z4)
        @test all(I4cp .== Z4)
        @test all(BLAS.gemm('N', 'N', I4, U4) .== U4)
        @test all(BLAS.gemm('N', 'T', I4, U4) .== L4)
        @test_throws DimensionMismatch BLAS.gemm!('N','N', one(elty), I4, I4, elm1, eye(elty,5))
        @test_throws DimensionMismatch BLAS.gemm!('N','N', one(elty), I43, I4, elm1, I4)
        @test_throws DimensionMismatch BLAS.gemm!('T','N', one(elty), I43, I4, elm1, I43)
        @test_throws DimensionMismatch BLAS.gemm!('N','T', one(elty), I43, I43, elm1, I43)
        @test_throws DimensionMismatch BLAS.gemm!('T','T', one(elty), I43, I43, elm1, I43')
    end
    @testset "gemm compared to (sy)(he)rk" begin
        if eltype(elm1) <: Complex
            @test all(triu(BLAS.herk('U', 'N', U4)) .== triu(BLAS.gemm('N', 'T', U4, U4)))
            @test all(tril(BLAS.herk('L', 'N', U4)) .== tril(BLAS.gemm('N', 'T', U4, U4)))
            @test all(triu(BLAS.herk('U', 'N', L4)) .== triu(BLAS.gemm('N', 'T', L4, L4)))
            @test all(tril(BLAS.herk('L', 'N', L4)) .== tril(BLAS.gemm('N', 'T', L4, L4)))
            @test all(triu(BLAS.herk('U', 'C', U4)) .== triu(BLAS.gemm('T', 'N', U4, U4)))
            @test all(tril(BLAS.herk('L', 'C', U4)) .== tril(BLAS.gemm('T', 'N', U4, U4)))
            @test all(triu(BLAS.herk('U', 'C', L4)) .== triu(BLAS.gemm('T', 'N', L4, L4)))
            @test all(tril(BLAS.herk('L', 'C', L4)) .== tril(BLAS.gemm('T', 'N', L4, L4)))
            ans = similar(L4)
            @test all(tril(BLAS.herk('L','C', L4)) .== tril(BLAS.herk!('L', 'C', real(one(elty)), L4, real(zero(elty)), ans)))
            @test all(Base.LinAlg.copytri!(ans, 'L') .== LinAlg.BLAS.gemm('T', 'N', L4, L4))
            @test_throws DimensionMismatch BLAS.herk!('L','N',real(one(elty)),eye(elty,5),real(one(elty)),eye(elty,6))
        else
            @test all(triu(BLAS.syrk('U', 'N', U4)) .== triu(BLAS.gemm('N', 'T', U4, U4)))
            @test all(tril(BLAS.syrk('L', 'N', U4)) .== tril(BLAS.gemm('N', 'T', U4, U4)))
            @test all(triu(BLAS.syrk('U', 'N', L4)) .== triu(BLAS.gemm('N', 'T', L4, L4)))
            @test all(tril(BLAS.syrk('L', 'N', L4)) .== tril(BLAS.gemm('N', 'T', L4, L4)))
            @test all(triu(BLAS.syrk('U', 'T', U4)) .== triu(BLAS.gemm('T', 'N', U4, U4)))
            @test all(tril(BLAS.syrk('L', 'T', U4)) .== tril(BLAS.gemm('T', 'N', U4, U4)))
            @test all(triu(BLAS.syrk('U', 'T', L4)) .== triu(BLAS.gemm('T', 'N', L4, L4)))
            @test all(tril(BLAS.syrk('L', 'T', L4)) .== tril(BLAS.gemm('T', 'N', L4, L4)))
            ans = similar(L4)
            @test all(tril(BLAS.syrk('L','T', L4)) .== tril(BLAS.syrk!('L', 'T', one(elty), L4, zero(elty), ans)))
            @test all(Base.LinAlg.copytri!(ans, 'L') .== BLAS.gemm('T', 'N', L4, L4))
            @test_throws DimensionMismatch BLAS.syrk!('L','N',one(elty),eye(elty,5),one(elty),eye(elty,6))
        end
    end
end

@testset "syr for eltype $elty" for elty in (Float32, Float64, Complex{Float32}, Complex{Float64})
    A = rand(elty, 5, 5)
    @test triu(A[1,:] * A[1,:].') ≈ BLAS.syr!('U', one(elty), A[1,:], zeros(elty, 5, 5))
    @test tril(A[1,:] * A[1,:].') ≈ BLAS.syr!('L', one(elty), A[1,:], zeros(elty, 5, 5))
    @test triu(A[1,:] * A[1,:].') ≈ BLAS.syr!('U', one(elty), view(A, 1, :), zeros(elty, 5, 5))
    @test tril(A[1,:] * A[1,:].') ≈ BLAS.syr!('L', one(elty), view(A, 1, :), zeros(elty, 5, 5))
end

@testset "her for eltype $elty" for elty in (Complex{Float32}, Complex{Float64})
    A = rand(elty, 5, 5)
    @test triu(A[1,:] * A[1,:]') ≈ BLAS.her!('U', one(real(elty)), A[1,:], zeros(elty, 5, 5))
    @test tril(A[1,:] * A[1,:]') ≈ BLAS.her!('L', one(real(elty)), A[1,:], zeros(elty, 5, 5))
    @test triu(A[1,:] * A[1,:]') ≈ BLAS.her!('U', one(real(elty)), view(A, 1, :), zeros(elty, 5, 5))
    @test tril(A[1,:] * A[1,:]') ≈ BLAS.her!('L', one(real(elty)), view(A, 1, :), zeros(elty, 5, 5))
end

max_stars_repo_name: caodi/julia

max_issues_repo_name: caodi/julia

max_forks_repo_name: caodi/julia


text:
# maybe put this script with the other one together?
export sudoku_1, sudoku_2
sudoku_1 = [0 0 4 7 2 0 9 0 0; 0 3 9 0 0 8 0 0 5; 0 0 1 5 0 6 0 0 4; 0 4 0 0 1 0 5 2 0; 0 2 8 0 5 0 1 7 0; 0 1 6 0 3 0 0 9 0; 4 0 0 9 0 1 3 0 0; 1 0 0 3 0 0 8 4 0; 0 0 7 0 8 5 6 0 0];

sudoku_2 = [5 3 0 0 7 0 0 0 0; 6 0 0 1 9 5 0 0 0; 0 9 8 0 0 0 0 6 0; 8 0 0 0 6 0 0 0 3; 4 0 0 8 0 3 0 0 1; 7 0 0 0 2 0 0 0 6;0 6 0 0 0 0 2 8 0; 0 0 0 4 1 9 0 0 5; 0 0 0 0 8 0 0 7 9];


text:
using DataFrames, Dates, Parameters

# Abstract type for dispatch
abstract type TreatmentPanel end

# Types for number of treatment units and periods
# Treatment duration - this is either continuous or discontinuous
abstract type TreatmentDurationType end
struct Continuous <: TreatmentDurationType end
struct Discontinuous <: TreatmentDurationType end

# Treatment timing - relevant only for MultiUnitTreatments, can be simultaneous or staggered
abstract type TreatmentTimingType end
struct Simultaneous{T <: TreatmentDurationType} <: TreatmentTimingType end
struct Staggered{T <: TreatmentDurationType} <: TreatmentTimingType end

# Unit type - either single or multiple treated units
abstract type TreatmentType end
struct SingleUnitTreatment{T <: TreatmentDurationType} <: TreatmentType end
struct MultiUnitTreatment{T <: TreatmentTimingType}  <: TreatmentType end

# BalancedPanel will have an N×T matrix of treatment assigment and outcomes
"""
    BalancedPanel{TreatmentType}

A TreatmentPanel in which all N treatment units are observed for the same T periods.

The object is constructed from a `DataFrame` which contains all outcome and covariate data. The
constructor requires passing information on which column in the `DataFrame` holds the subject and
time period identifiers, as well as information on the timing of treatment. Treatments are specified
as a `Pair` of either a `String` or `Symbol` identifying the unit treated, and a `Date` or
`Int` value, identifying the treatment period. Where treatments have an end point, the
treatment period os specified as a `Tuple` of either `Date` or `Int`, indicating the first
and last period of treatment. Where individual units have multiple treatment periods,
these are specified as a `Pair{Union{String, Symbol}, Vector{Union{Int, Date}}}` of treatment
unit identifier and vector of timings. Finally, single continuous, single periodic, and
multiple period treatments can be generalized to multiple treated units.
The following table provides an overview of the types of treatment pattern supported:

|                     |  Only starting point        |   Start and end point                     |   Multiple start & end points                       |
|---------------------|-----------------------------|-------------------------------------------|-----------------------------------------------------|
| **one unit**        |  Pair{String, Date}         |   Pair{String, Tuple{Date, Date}}         |  Pair{String}, Vector{Tuple{Date, Date}}}           |
| **multiple units**  |  Vector{Pair{String, Date}} |   Vector{Pair{String, Tuple{Date, Date}}} |  Vector{Pair{String}, Vector{Tuple{Date, Date}}}}   |
"""
@with_kw struct BalancedPanel{UTType} <: TreatmentPanel where UTType <: TreatmentType
    W::Union{Matrix{Bool}, Matrix{Union{Missing, Bool}}}
    Y::Matrix{Float64}
    df::DataFrame
    id_var::Union{String, Symbol}
    t_var::Union{String, Symbol}
    outcome_var::Union{String, Symbol}
    ts::Vector{T1} where T1 <: Union{Date, Int64}
    is::Vector{T2} where T2 <: Union{Symbol, String, Int64}
end

# Check that ID, time, and outcome variable are provided
function check_id_t_outcome(df, outcome_var, id_var, t_var)
    # Check relevant info has been provided
    !isnothing(outcome_var) || error(ArgumentError(
        "Please specify outcome_var, the name of the column in your dataset holding the "*
        "outcome variable of interest in your dataset."
        ))
    !isnothing(id_var) || error(ArgumentError(
            "Please specify id_var, the name of the column in your data set holding the "*
            "identifier of your units of observation (panel dimension)."
            ))
    !isnothing(t_var) || error(ArgumentError(
            "Please specify t_var, the name of the column in your dataset holding the "*
            "time dimension."
        ))

    # Ensure columns exist in data
    f = in(names(df))
    f(string(id_var)) || throw("Error: ID variable $id_var is not present in the data.")
    f(string(t_var)) || throw("Error: Time variable $t_var is not present in the data.")
    f(string(outcome_var)) || throw("Error: ID variable $outcome_var is not present in the data.")
end

# Functions to get all treatment periods
function treatment_periods(ta::Pair{T1, S1}) where T1 where S1
    [last(ta)]
end

function treatment_periods(ta::Pair{T1, S1}) where T1 where S1 <: Union{Pair{Int, Int}, Pair{Date, Date}}
    collect(last(ta))
end

function treatment_periods(ta::Vector{Pair{T1, S1}}) where T1 where S1
    last.(ta)
end

function treatment_periods(ta::Vector{Pair{T1, S1}}) where T1 where S1 <: Union{Pair{Int, Int}, Pair{Date, Date}}
    unique(reduce(vcat, collect.(last.(ta))))
end

# Functions to construct treatment assignment matrix
function construct_W(ta::Pair{T1, S1}, N, T, is, ts) where T1 where S1
    W = [false for i = 1:N, j = 1:T]
    W[findfirst(==(ta[1]), is), findfirst(==(ta[2]), ts):end] .= true

    return W
end

function construct_W(ta::Pair{T1, S1}, N, T, is, ts) where T1 where S1 <: Union{Pair{Int, Int}, Pair{Date, Date}}
    W = [false for i = 1:N, j = 1:T]
    W[findfirst(==(ta[1]), is), findfirst(==(ta[2][1]), ts):findfirst(==(ta[2][2]), ts)] .= true

    return W
end

function construct_W(tas::Vector{Pair{T1, S1}}, N, T, is, ts) where T1 where S1
    W = [false for i = 1:N, j = 1:T]
    for ta ∈ tas
        W[findfirst(==(ta[1]), is), findfirst(==(ta[2]), ts):end] .= true
    end

    return W
end

function construct_W(tas::Vector{Pair{T1, S1}}, N, T, is, ts) where T1 where S1 <: Union{Pair{Int, Int}, Pair{Date, Date}}
    W = [false for i = 1:N, j = 1:T]
    for ta ∈ tas
        W[findfirst(==(ta[1]), is), findfirst(==(ta[2][1]), ts):findfirst(==(ta[2][2]), ts)] .= true
    end

    return W
end

# Constructor for single continuous treatment - returns BalancedPanel{SingleUnitTreatment{Continuous}}
function BalancedPanel(df::DataFrame, treatment_assignment::Pair{T1, T2};
    id_var = nothing, t_var = nothing, outcome_var = nothing,
    sort_inplace = false) where T1 where T2 <: Union{Date, Int}

    # Get all units and time periods
    is = sort(unique(df[!, id_var])); i_set = Set(is)
    ts = sort(unique(df[!, t_var])); t_set = Set(ts)

    # Dimensions
    N = length(is)
    T = length(ts)

    # Get all treatment units and treatment periods
    treated_i = first(treatment_assignment)
    treated_t = last(treatment_assignment)

    # Sanity checks
    check_id_t_outcome(df, outcome_var, id_var, t_var)
    in(treated_i, i_set) || throw("Error: Treatment unit $treated_i is not in the list of unit identifiers $id_var")
    in(treated_t, t_set) || throw("Error: Treatment period $treated_t is not in the list of time identifiers $t_var")

    # Sort data if necessary, in place if required
    df = ifelse(issorted(df, [id_var, t_var]), df,
                                               ifelse(sort_inplace, sort!(df, [id_var, t_var]),
                                                                    sort(df, [id_var, t_var])))

    # Treatment matrix
    W = construct_W(treatment_assignment, N, T, is, ts)

    # Outcome matrix
    Y = zeros(eltype(df[!, outcome_var]), size(W))
    for (row, i) ∈ enumerate(is), (col, t) ∈ enumerate(ts)
        Y[row, col] = only(df[(df[!, id_var] .== i) .& (df[!, t_var] .== t), outcome_var])
    end

    BalancedPanel{SingleUnitTreatment{Continuous}}(W, Y, df, id_var, t_var, outcome_var, ts, is)
end

# Getter functions
"""
    treated_ids(x <: BalancedPanel)

    Returns the indices of treated units in the panel, so that Y[treated_ids(x), :] returns a
    (Nₜᵣ×T) matrix of outcomes for treated units in all periods.
"""
function treated_ids(x::BalancedPanel{SingleUnitTreatment{T}}) where T
    for i ∈ 1:size(x.Y, 1)
        for t ∈ 1:size(x.Y, 2)
            if x.W[i, t]
                return i
            end
        end
    end
end

function treated_ids(x::BalancedPanel{MultiUnitTreatment{T}}) where T
    findall(>(0), vec(sum(Y, dims = 2)))
end

"""
    treated_labels(x <: BalancedPanel)

    Returns the labels of treated units as given by the `id_var` column in the underlying data set.
"""
function treated_labels(x::BalancedPanel{SingleUnitTreatment{T}}) where T
    x.is[treated_ids(x)]
end

"""
    first_treated_period_ids(x <: BalancedPanel)

    Returns the indices of the first treated period for each treated units, that is, a Vector{Int}
    of length Nₜᵣ, where each element is the index of the first 1 in the row of treatment matrix W
    corresonding to the treatment unit.
"""
function first_treated_period_ids(x::BalancedPanel{SingleUnitTreatment{T}}) where T
    findfirst(vec(x.W[treated_ids(x), :]))
end

"""
    first_treated_period_labels(x <: BalancedPanel)

    Returns the labels of the first treated period for each treated units, that is, a Vector{T}
    of length Nₜᵣ, where T is the eltype of the `t_var` column in the underlying data.
"""
function first_treated_period_labels(x::BalancedPanel{SingleUnitTreatment{T}}) where T
    x.ts[first_treated_period_ids(x)]
end

"""
    length_T₀(x <: BalancedPanel)

    Returns the number of pre-treatment periods for each treated unit.
"""
function length_T₀(x::BalancedPanel{SingleUnitTreatment{Continuous}})
    first_treated_period_ids(x) - 1
end

"""
    length_T₁(x <: BalancedPanel)

    Returns the number of treatment periods for each treated unit.
"""
function length_T₁(x::BalancedPanel{SingleUnitTreatment{Continuous}})
    size(x.Y, 2) .- first_treated_period_ids(x) + 1
end

"""
    get_y₁₀(x <: BalancedPanel)

    Returns the pre-treatment outcomes for the treated unit(s). For SingleUnitTreatment designs,
    this is a vector of length T₀, while for MultiUnitTreatment designs, it is a (Nₜᵣ×T₀) matrix
"""
function get_y₁₀(x::BalancedPanel{SingleUnitTreatment{Continuous}})
    x.Y[treated_ids(x), 1:first_treated_period_ids(x)-1]
end

"""
    get_y₁₁(x <: BalancedPanel)

    Returns the post-treatment outcomes for the treated unit(s). For SingleUnitTreatment designs,
    this is a vector of length T₁, while for MultiUnitTreatment designs, it is a (Nₜᵣ×T₁) matrix
"""
function get_y₁₁(x::BalancedPanel{SingleUnitTreatment{Continuous}})
    x.Y[x.W]
end

""" sc
    get_y₀₀(x <: BalancedPanel)

    Returns the pre-treatment outcomes for the untreated units, an (Nₖₒ×T₀) matrix
"""
function get_y₀₀(x::BalancedPanel{SingleUnitTreatment{Continuous}})
    x.Y[Not(treated_ids(x)), 1:first_treated_period_ids(x)-1]
end

"""
    get_y₀₀(x <: BalancedPanel)

    Returns the post-treatment outcomes for the untreated units, an (Nₖₒ×T₁) matrix
"""
function get_y₀₁(x::BalancedPanel{SingleUnitTreatment{Continuous}})
    x.Y[Not(treated_ids(x)), first_treated_period_ids(x):end]
end

"""
    decompose_y(x <: BalancedPanel)

    Decomposes the outcome matrix Y into four elements:

    * Pre-treatment outcomes for treated units (y₁₀)
    * Post-treatment outcomes for treated units (y₁₁)
    * Pre-treatment outcomes for control units (y₀₀)
    * Post-treatment outcomes for treated units (y₀₁)

    and returns a tuple (y₁₀, y₁₁, y₀₀, y₀₁)
"""
function decompose_y(x)
    get_y₁₀(x), get_y₁₁(x), get_y₀₀(x), get_y₀₁(x)
end

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

# Constructor for single start/end treatment - returns BalancedPanel{SingleUnitTreatment{Discontinuous}}
function BalancedPanel(df::DataFrame, treatment_assignment::Pair{T1, T2};
    id_var = nothing, t_var = nothing, outcome_var = nothing,
    sort_inplace = false) where T1 where T2 <: Union{Pair{Date, Date}, Pair{Int, Int}}

    # Get all units and time periods
    is = sort(unique(df[!, id_var])); i_set = Set(is)
    ts = sort(unique(df[!, t_var])); t_set = Set(ts)

    # Dimensions
    N = length(is)
    T = length(ts)

    # Get all treatment units and treatment periods
    treated_i = first(treatment_assignment)
    treatment_start_end = last(treatment_assignment)

    # Sanity checks
    check_id_t_outcome(df, outcome_var, id_var, t_var)
    in(treated_i, i_set) || throw("Error: Treatment unit $treated_i is not in the list of unit identifiers $id_var")
    in(treatment_start_end[1], t_set) || throw("Error: Treatment start point $(treatment_start_end[1]) is not in the list of time identifiers $t_var")
    in(treatment_start_end[2], t_set) || throw("Error: Treatment end point $(treatment_start_end[2]) is not in the list of time identifiers $t_var")

    # Sort data if necessary, in place if required
    df = ifelse(issorted(df, [id_var, t_var]), df,
                                               ifelse(sort_inplace, sort!(df, [id_var, t_var]),
                                                                    sort(df, [id_var, t_var])))

    # Treatment matrix
    W = construct_W(treatment_assignment, N, T, is, ts)

    # Outcome matrix
    Y = zeros(size(W))
    for (row, i) ∈ enumerate(is), (col, t) ∈ enumerate(ts)
        try
            Y[row, col] = only(df[(df[!, id_var] .== i) .& (df[!, t_var] .== t), outcome_var])
        catch ArgumentError
            throw("$(nrow(df[(df[!, id_var] .== i) .& (df[!, t_var] .== t), :])) outcomes present in the data for unit $i in period $t")
        end
    end

    BalancedPanel{SingleUnitTreatment{Discontinuous}}(W, Y, df, id_var, t_var, outcome_var, ts, is)
end

# Fallback method - if the length of treatment assignment is one use single treatment method above
function BalancedPanel(df::DataFrame, treatment_assignment;
    id_var = nothing, t_var = nothing, outcome_var = nothing,
    sort_inplace = false) where NType where TType

    if length(treatment_assignment) == 1
        return BalancedPanel(df, only(treatment_assignment); id_var = id_var, t_var = t_var,
                                outcome_var = outcome_var, sort_inplace = sort_inplace)
    end

    # Get all units and time periods
    is = sort(unique(df[!, id_var])); i_set = Set(is)
    ts = sort(unique(df[!, t_var])); t_set = Set(ts)

    # Get all treatment units and treatment periods
    treated_is = first.(treatment_assignment)
    treated_is = typeof(treated_is) <: AbstractArray ? treated_is : [treated_is]
    treated_ts = treatment_periods(treatment_assignment)

    # Dimensions
    N = length(is)
    T = length(ts)

    ### SANITY CHECKS ###
    check_id_t_outcome(df, outcome_var, id_var, t_var)
    for ti ∈ treated_is
        in(ti, i_set) || throw("Error: Treatment unit $ti is not in the list of unit identifiers $id_var")
    end

    for tt ∈ treated_ts
        in(tt, t_set) || throw("Error: Treatment period $tt is not in the list of time identifiers $t_var")
    end

    # Sort data if necessary, in place if required
    df = ifelse(issorted(df, [id_var, t_var]), df,
                                               ifelse(sort_inplace, sort!(df, [id_var, t_var]),
                                                                    sort(df, [id_var, t_var])))

    # Treatment matrix
    W = construct_W(treatment_assignment, N, T, is, ts)

    # Outcome matrix
    Y = zeros(eltype(df[!, outcome_var]), size(W))
    for (row, i) ∈ enumerate(is), (col, t) ∈ enumerate(ts)
        Y[row, col] = only(df[(df[!, id_var] .== i) .& (df[!, t_var] .== t), outcome_var])
    end

    # Determine TreatmentType and TreatmentDurationType
    uttype = if all(==(treatment_assignment[1][2]), last.(treatment_assignment))
        Simultaneous
    else
        Staggered
    end

    tdtype = if typeof(treatment_assignment) <: Pair
        if typeof(treatment_assignment[2]) <: Pair
            Discontinuous
        else
            Continuous
        end
    else
        if typeof(treatment_assignment[1][2]) <: Pair
            Discontinuous
        else
            Continuous
        end
    end

    BalancedPanel{MultiUnitTreatment{uttype{tdtype}}}(W, Y, df, id_var, t_var, outcome_var, ts, is)
end

## UnblancedPanel - N observations but not all of them for T periods

#!# Not yet implemented


text:
#=
The Vectors defined here are necessary to
ClapeyronParam
=#
struct Compressed4DMatrix{T,V<:AbstractVector{T}}
    values::V
    outer_indices::Vector{Tuple{Int,Int}} #index of components
    inner_indices::Vector{Tuple{Int,Int}} #index of sites
    outer_size::Tuple{Int,Int} #size of component matrix
    inner_size::Tuple{Int,Int} #size of sites matrices
end

function Base.show(io::IO,mime::MIME"text/plain",m::Compressed4DMatrix{T}) where T
    n = length(m.values)
    println(io,typeof(m)," with ",n," entr",(n == 1 ? "y:" : "ies:"))
    for (idx,(i,j),(a,b)) in indices(m)
        if idx != 1
        println(io)
        end
        print(io," ",(i,a)," >=< ",(j,b),": ",m.values[idx])
    end
end

function Base.show(io::IO,m::Compressed4DMatrix{T}) where T
    print(io,typeof(m))
    print(io,m.values)
end
function Compressed4DMatrix{T}() where T
    return Compressed4DMatrix(T[],Tuple{Int,Int}[],Tuple{Int,Int}[],(0,0),(0,0))
end

const MatrixofMatrices{T} = AbstractMatrix{<:AbstractMatrix{T}} where T

function Compressed4DMatrix(x::MatrixofMatrices{T}) where T
    outer_size = size(x)
    is1, is2 = 0, 0
    os1, os2 = outer_size
    for i in eachindex(x)
        _is1,_is2 = size(x[i])
        is1,is2 = max(_is1,is1),max(_is2,is2)
    end
    inner_size = (is1,is2)

    values = T[]
    inner_indices = Tuple{Int,Int}[]
    outer_indices = Tuple{Int,Int}[]

    if iszero(os1) & iszero(os2)
        return Compressed4DMatrix(values,outer_indices,inner_indices,outer_size,inner_size)
    end

    #self association

    for i in 1:os1
        for j in i:os2 #there can be self association
            xi = x[i,j]
            if iszero(length(xi))
                continue
            end
            a1,a2 = size(xi)
            for a in 1:a1
                start = ifelse(i == j,a,1)
                for b in start:a2 #this includes (i,i)(a,a) (sCKSAFT)
                    if !_iszero(xi[a,b])
                        push!(values,xi[a,b])
                        push!(outer_indices,(i,j))
                        push!(inner_indices,(a,b))
                    end
                end
            end
        end
    end
    return Compressed4DMatrix{T,Vector{T}}(values,outer_indices,inner_indices,outer_size,inner_size)
end

function Base.getindex(m::Compressed4DMatrix,i::Int,j::Int)
    i,j = minmax(i,j)
    @inbounds begin
    idx = searchsorted(m.outer_indices,(i,j))
    return AssocView(view(m.values,idx),view(m.inner_indices,idx),m.inner_size)
    end
end

#Base.eltype(m::Compressed4DMatrix{T}) where T = T

function Base.setindex!(m::Compressed4DMatrix,val,i::Int)
    @inbounds begin
        m.values[i] = val
    end
end

struct AssocView{T,V,I} <: AbstractMatrix{T}
    values::V
    indices::I
    size::Tuple{Int,Int}
end

function AssocView(V::Vi,I::Ii,s) where Ii where Vi <:AbstractVector{<:T} where T
    return AssocView{T,Vi,Ii}(V,I,s)
end

Base.size(m::AssocView) = m.size
Base.eltype(m::AssocView{T}) where T = T

function Base.getindex(m::AssocView{T},i::Int,j::Int) where T
    @inbounds begin
        idxs = searchsorted(m.indices,(i,j))
        if iszero(length(idxs))
            idxs = searchsorted(m.indices,(j,i))
            iszero(length(idxs)) &&  return zero(T)
        end
        idx = first(idxs)
        return m.values[idx]
    end
end

function zero_assoc(m::Compressed4DMatrix,::Type{T} = Float64) where T <:Number
    newvalues = zeros(T,length(m.values))
    return Compressed4DMatrix(newvalues,m.outer_indices,m.inner_indices,m.outer_size,m.inner_size)
end

function zero_assoc(m::Matrix{<:Matrix},a::Type{T} = Float64) where T <:Number
    return zero_assoc(Compressed4DMatrix(m),a)
end

#=
function zero_assoc(m::Compressed4DMatrix,x::T = 0.0) where T <:Number
    newvalues = fill(x,length(m.values))
    return Compressed4DMatrix(newvalues,m.outer_indices,m.inner_indices,m.outer_size,m.inner_size)
end
=#
function indices(x::Compressed4DMatrix)
    return zip(1:length(x.values),x.outer_indices,x.inner_indices)
end

function indices(x::PackedVofV)
    return x.p
end

struct SparsePackedMofV{E,P<:PackedVofV}<:SparseArrays.AbstractSparseMatrixCSC{E,Int}
    storage::P
    idx::SparseMatrixCSC{Int,Int}
end

function SparsePackedMofV(storage,idx)
    E = eltype(storage)
    P = typeof(storage)
    return SparsePackedMofV{E,P}(storage,idx)
end

_findnz(x::SparseMatrixCSC) = SparseArrays.findnz(x)
function _findnz(x::AbstractMatrix{<:AbstractVector})
    idxs = findall(z->!iszero(length(z)),x)
    i = first.(idxs)
    j = last.(idxs)
    vals = x[idxs]
    return i,j,vals
end

function SparsePackedMofV(m_of_v::AbstractMatrix{<:AbstractVector})
    i,j,unpack_sparse_vals = _findnz(m_of_v)
    #unpack_sparse_vals = m_of_v[sparse_idx]
    pack_sparse_vals = PackedVectorsOfVectors.pack(unpack_sparse_vals)
    len_linear = length(unpack_sparse_vals)
    idx_linear = collect(1:len_linear)
    m,n = size(m_of_v)
    idx = sparse(i,j,idx_linear,m,n)
    E = eltype(pack_sparse_vals)
    P = typeof(pack_sparse_vals)
    return SparsePackedMofV{E,P}(pack_sparse_vals,idx)
end

function Base.getindex(x::SparsePackedMofV,i::Int,j::Int)
    _idx = x.idx[i,j]
    iszero(_idx) && (return view(x.storage.v,1:0))

    return x.storage[_idx]
end
@inline Base.size(x::SparsePackedMofV) = size(x.idx)
@inline SparseArrays.nnz(x::SparsePackedMofV) = length(x.storage.p)
@inline function SparseArrays.findnz(x::SparsePackedMofV)
    i,j,_ = SparseArrays.findnz(x.idx)
    return i,j,x.storage
end
@inline SparseArrays.nonzeros(x::SparsePackedMofV) = x.storage
@inline SparseArrays.rowvals(x::SparsePackedMofV) = SparseArrays.rowvals(x.idx)
@inline SparseArrays.nzrange(A::SparsePackedMofV, col::Integer) = SparseArrays.nzrange(A.idx,col)
@inline SparseArrays.getcolptr(A::SparsePackedMofV) = SparseArrays.getcolptr(A.idx)
export SparsePackedMofV

function Base.show(io::IO,::MIME"text/plain",A::SparsePackedMofV)
    m,n = size(A)
    println(io,"$(m)×$(n) Sparse Packed Matrix of Vectors of eltype $(eltype(A.storage)) with $(length(A.storage)) non empty values:")
    vals = A.storage
    rows = rowvals(A)
    for j in 1:n
        for ii ∈ nzrange(A, j)
            i = rows[ii]
            val= vals[ii]
            println(io,"  ($i,$j) => $val")
        end
    end
end


text:
include("/root/.julia/v0.3/Gadfly/src/Gadfly.jl")
importall Gadfly

max_stars_repo_name: Kaggle/docker-julia

max_issues_repo_name: AI-For-Rural/docker-julia

max_forks_repo_name: AI-For-Rural/docker-julia


text:
"""
    compute_frames!(elements::Vector{UA}; parent=nothing)
        where UA<:Union{AbstractObject,AbstractAction}

Set elem.frames.frames to the computed frames for each elem in elements.
"""
function compute_frames!(
    elements::Vector{UA};
    parent = nothing,
    parent_counter = 0,
) where {UA<:Union{AbstractObject,AbstractAction}}
    available_subframes = typemin(Int):typemax(Int)
    if parent !== nothing
        last_frames = get_frames(parent)
        available_subframes = 1:length(get_frames(parent))
    else
        last_frames = nothing
    end
    is_first = true
    counter = 1
    for elem in elements
        if last_frames === nothing && get_frames(elem) === nothing
            throw(ArgumentError("Frames need to be defined explicitly in the initial
                Object/Background or Action."))
        end
        if get_frames(elem) === nothing
            set_frames!(parent, elem, last_frames; is_first = is_first)
        end
        last_frames = get_frames(elem)
        if !(get_frames(elem) ⊆ available_subframes)
            @warn("Action defined outside the frame range of the parent object.
            Action #$counter for Object #$parent_counter is defined for frames
            $(get_frames(elem)) but Object #$parent_counter exists only for $(available_subframes).
            (Info: Background is counted as Object #1)")
        end
        is_first = false
        counter += 1
    end
end

"""
    get_current_setting()

Return the current setting of the current object
"""
function get_current_setting()
    object = CURRENT_OBJECT[1]
    return object.current_setting
end

"""
    get_interpolation(frames::UnitRange, frame)

Return a value between 0 and 1 which represents the relative `frame` inside `frames`.
"""
function get_interpolation(frames::UnitRange, frame)
    frame == last(frames) && return 1.0
    t = (frame - first(frames)) / (length(frames) - 1)
    # makes sense to only allow 0 ≤ t ≤ 1
    t = min(1.0, t)
end

"""
    get_interpolation(action::AbstractAction, frame)

Return the value of the `action.anim` Animation based on the relative frame given by
`get_interpolation(get_frames(action), frame)`
"""
function get_interpolation(action::AbstractAction, frame)
    t = get_interpolation(get_frames(action), frame)
    if !(action.anim.frames[end].t ≈ 1)
        @warn "Animations should be defined from 0.0 to 1.0"
    end
    return interpolation_to_transition_val(at(action.anim, t), action.transition)
end

"""
    interpolation_to_transition_val(interpolation_val, Transition)

Returns the transition value for the given `interpolation_val`.
If the interpolation value is already of the correct form it just gets returned.
Otherwise the Transition function like `get_position` is called and the interpolated value
is calculated.
"""
interpolation_to_transition_val(t, ::Nothing) = t
interpolation_to_transition_val(t::Point, trans::Translation) = t
interpolation_to_transition_val(t::Float64, trans::Rotation) = t
interpolation_to_transition_val(t::Scale, trans::Scaling) = t

function interpolation_to_transition_val(t, trans::Translation)
    from = get_position(trans.from)
    to = get_position(trans.to)
    return from + t * (to - from)
end

function interpolation_to_transition_val(t, trans::Scaling)
    from = get_scale(trans.from)
    to = get_scale(trans.to)
    return from + t * (to - from)
end

function isapprox_discrete(val; atol = 1e-4)
    return isapprox(val, round(val); atol = atol)
end

function polywh(polygon::Vector{Vector{Point}})
    T = typeof(polygon[1][1].x)
    min_x = typemax(T)
    min_y = typemax(T)
    max_x = typemin(T)
    max_y = typemin(T)
    for poly in polygon
        for p in poly
            min_x = min(min_x, p.x)
            min_y = min(min_y, p.y)
            max_x = max(max_x, p.x)
            max_y = max(max_y, p.y)
        end
    end
    return max_x - min_x, max_y - min_y
end

function get_polypoint_at(points, t; pdist = polydistances(points))
    if t ≈ 0
        return points[1]
    end
    ind, surplus = nearestindex(pdist, t * pdist[end])

    nextind = mod1(ind + 1, length(points))
    overshootpoint = between(
        points[ind],
        points[nextind],
        surplus / distance(points[ind], points[nextind]),
    )
    return overshootpoint
end