text:
# Forked from Julia Jupyter notebook 2019-12-PolaronMultipleBranches.ipynb

using PolaronMobility
using Test

using Plots
using Printf

const hbar = const ħ = 1.05457162825e-34;          # kg m2 / s
const Boltzmann = const kB =  1.3806504e-23;                  # kg m2 / K s2
const eV = const q = const ElectronVolt = 1.602176487e-19;                         # kg m2 / s2
const MassElectron = 9.10938188e-31;                          # kg

# ((freq THz)) ((IR Activity / e^2 amu^-1))
# These data from MAPbI3-Cubic_PeakTable.csv
# https://github.com/WMD-group/Phonons/tree/master/2015_MAPbI3/SimulatedSpectra
# Data published in Brivio2015 (PRB)
# https://doi.org/10.1103/PhysRevB.92.144308
MAPI= [
#96.20813558773261 0.4996300522819191
#93.13630357703363 1.7139631746083817
#92.87834578121567 0.60108592692181
#92.4847918585963 0.0058228799414729
#92.26701437594754 0.100590086574602
#89.43972834606603 0.006278895133832249
#46.89209141511332 0.2460894564364346
#46.420949316788 0.14174282581124137
#44.0380222871706 0.1987196948553428
#42.89702947649343 0.011159939465770681
#42.67180170168193 0.02557751102757614
#41.46971205834201 0.012555230726601503
#37.08982543385215 0.00107488277468418
#36.53555265689563 0.02126940080871224
#30.20608114002676 0.009019481779712388
#27.374810898415028 0.03994453721421388
#26.363055017011728 0.05011922682554448
#9.522966890022039 0.00075631870522737
4.016471586720514 0.08168931020200264
3.887605410774121 0.006311654262282101
3.5313112232401513 0.05353548710183397
2.755392921480459 0.021303020776321225
2.4380741812443247 0.23162784335484837
2.2490917637719408 0.2622203718355982
2.079632190634424 0.23382298607799906
2.0336707697261187 0.0623239656843172
1.5673011873879714 0.0367465760261409
1.0188379384951798 0.0126328938653956
1.0022960504442775 0.006817361620021601
0.9970130778462072 0.0103757951973341
0.9201781906386209 0.01095811116040592
0.800604081794174 0.0016830270365341532
0.5738689505255512 0.00646428491253749
#0.022939578929507105 8.355742795827834e-05   # Acoustic modes!
#0.04882611767873102 8.309858592685e-06
#0.07575149723846182 2.778248540373041e-05
]

vol=(6.29E-10)^3
ϵ_o=4.5
meff=0.12
mb=meff*MassElectron

ϵ_i=IRtoDielectric(MAPI,vol)
ϵ_s=sum(ϵ_i)+ϵ_o # total (static) dielectric = sum of ionic, and optical

IRtoalpha(MAPI, volume=vol, ϵ_o=ϵ_o, ϵ_s=ϵ_s, meff=meff)

ϵ_polar_modes=DielectricFromIRmode.(eachrow(MAPI), volume=vol)
f_dielectric=hcat( MAPI[:,1], ϵ_polar_modes)
alphas=frohlichPartial.(eachrow(f_dielectric), ϵ_o = ϵ_o, ϵ_s = ϵ_o+sum(ϵ_polar_modes), meff=meff)

# First idea: solve for each alpha_i as a separate variational problem, then
# aggregate solutions
function separate_variational_solutions(alphas, MAPI)
    mobilityproblem=hcat(alphas, feynmanvw.(alphas), MAPI[:,1])
    inverse_μ=Hellwarth1999mobilityRHS.(eachrow(mobilityproblem), meff, 300)
    μ=sum(inverse_μ)^-1
    @printf("\n\tμ(Hellwarth1999)= %f m^2/Vs \t= %.2f cm^2/Vs",μ,μ*100^2)

    # as published in FrostPolaronMobility2017 PRB
    polaronmobility(300, 4.5, 24.1, 2.25E12, 0.12)

    # re-use these data, directly from the polaron problem
    polaronmobility(300, ϵ_o, ϵ_s, 2.25E12, 0.12)

    for T in 10:10:500
        inverse_μ=Hellwarth1999mobilityRHS.(eachrow(mobilityproblem), meff, T)
        μ=sum(inverse_μ)^-1
        @printf("\nT %f \tμ(Hellwarth1999)= %f m^2/Vs \t= %.2f cm^2/Vs",T,μ,μ*100^2)
    end

    normal=polaronmobility(10:10:500, ϵ_o, ϵ_s, 2.25E12, 0.12)

    hcat(normal.T, normal.Hμ)

    # OK, we have a problem I think in that we are optimising different things here
end

# Second idea: solve for variational problem with explicit set of interacting
# modes in the 'B' component of the free energy

# copy + paste from the code
using QuadGK
# Define Osaka's free-energies (Hellwarth1999 version) as Julia functions
# Equation numbers follow above Hellwarth et al. 1999 PRB
# 62b
A(v,w,β)=3/β*( log(v/w) - 1/2*log(2*π*β) - log(sinh(v*β/2)/sinh(w*β/2)))
# 62d
Y(x,v,β)=1/(1-exp(-v*β))*(1+exp(-v*β)-exp(-v*x)-exp(v*(x-β)))
# 62c integrand
#   Nb: Magic number 1e-10 adds stablity to optimisation; v,w never step -ve
f(x,v,w,β)=(exp(β-x)+exp(x))/sqrt(1e-10+ w^2*x*(1-x/β)+Y(x,v,β)*(v^2-w^2)/v)
# 62c
B(v,w,β,α) = α*v/(sqrt(π)*(exp(β)-1)) * quadgk(x->f(x,v,w,β),0,β/2)[1]
# 62e
C(v,w,β)=3/4*(v^2-w^2)/v * (coth(v*β/2)-2/(v*β))
# 62a
F(v,w,β,α)=-(A(v,w,β)+B(v,w,β,α)+C(v,w,β))

β(T,ω)=ħ*ω/(kB*T)

ω_hellwarth=2π*1E12*HellwarthBScheme(MAPI) # just low frequency branches, otherwise produces too high a value

F(v,w,ω_hellwarth,T,α_s,ω_s)=-(A(v,w,β(T,ω_hellwarth))+sum(B.(v,w,β.(T,ω_s),α_s))+C(v,w,β(T,ω_hellwarth)))

# OK, let's try this with a single value
@test F(7.6,6.5, ω_hellwarth, 300, 2.4, ω_hellwarth) ≈ F(7.6,6.5, β(300,ω_hellwarth), 2.4)
# And compare to the original route

# Split the strength, should still give same result
@test F(7.6,6.5, ω_hellwarth, 300, [1.2, 1.2], [ω_hellwarth, ω_hellwarth]) ≈ F(7.6,6.5, β(300,ω_hellwarth), 2.4)

# repack / rename variables from above
α_s=alphas
ω_s=2π*1E12*MAPI[:,1]

sum(α_s)

import PolaronMobility.feynmanvw
using Optim

function feynmanvw(T, ω_hellwarth, α_s, ω_s; v=7.1, w=6.5, verbose::Bool=true) # v,w defaults
    # Initial v,w to use
    initial=[v,w]
    # Main use of these bounds is stopping v or w going negative, at which you get a NaN error as you are evaluating log(-ve Real)
    lower=[0.1,0.1]
    upper=[100.0,100.0]

    myf(x) = F(x[1],x[2],ω_hellwarth, T, α_s, ω_s)
    # Wraps the function so just the two variational params are exposed, so Optim can call it

    # Now updated to use Optim > 0.15.0 call signature (Julia >0.6 only)
    res=optimize(OnceDifferentiable(myf, initial; autodiff = :forward),
                 lower, upper, initial, Fminbox( BFGS() ))
#                ,Optim.Options(g_tol=1e-15, allow_f_increases=true))
    # specify Optim.jl optimizer. This is doing all the work.

    if Optim.converged(res) == false
        print("\tWARNING: Failed to converge to v,w soln? : ",Optim.converged(res) )
    end

    if verbose # pretty print Optim solution
        println()
        show(res)
    end

    v,w=Optim.minimizer(res)
    return v,w
end

function HellwarthMobility(v,w,βred, α, ω)
# Hellwarth1999 - directly do contour integration in Feynman1962, for
        # finite temperature DC mobility
        # Hellwarth1999 Eqn (2) and (1) - These are going back to the general
        # (pre low-T limit) formulas in Feynman1962.  to evaluate these, you
        # need to do the explicit contour integration to get the polaron
        # self-energy
        R=(v^2-w^2)/(w^2*v) # inline, page 300 just after Eqn (2)

        #b=R*βred/sinh(b*βred*v/2) # This self-references b! What on Earth?
        # OK! I now understand that there is a typo in Hellwarth1999 and
        # Biaggio1997. They've introduced a spurious b on the R.H.S. compared to
        # the original, Feynman1962:
        b=R*βred/sinh(βred*v/2) # Feynman1962 version; page 1010, Eqn (47b)

        a=sqrt( (βred/2)^2 + R*βred*coth(βred*v/2))
        k(u,a,b,v) = (u^2+a^2-b*cos(v*u))^(-3/2)*cos(u) # integrand in (2)
        K=quadgk(u->k(u,a,b,v),0,Inf)[1] # numerical quadrature integration of (2)

        #Right-hand-side of Eqn 1 in Hellwarth 1999 // Eqn (4) in Baggio1997
        RHS= α/(3*sqrt(π)) * βred^(5/2) / sinh(βred/2) * (v^3/w^3) * K
        μ=RHS^-1 * (q)/(ω*mb)
        @printf("\n\tμ(Hellwarth1999)= %f m^2/Vs \t= %.2f cm^2/Vs",μ,μ*100^2)
        μ
end

fout=open("MultipleBranches-MAPI.dat","w")
@printf(fout,"T ExplicitBranches-v w mob HellwarthBScheme-v w mob")

Trange=10:1:500
for T in Trange
    # With explicit modes
    v,w=feynmanvw(T, ω_hellwarth, α_s, ω_s)
    μ=HellwarthMobility(v,w, β(T, ω_hellwarth), 2.4, ω_hellwarth)
    @printf(fout,"%d    %f %f %f",T,v,w,10_000*μ)

    # With Hellwarth effective mode frequency
    v,w=feynmanvw(T, ω_hellwarth, 2.4, ω_hellwarth)
    μ=HellwarthMobility(v,w, β(T, ω_hellwarth), 2.4, ω_hellwarth)
    @printf(fout," %f %f %f\n",v,w,10_000*μ)

    println("T: $(T) v: $(v) w:$(w)") # Just for STDOUT so you can monitor what's going on...
end

close(fout)


text:
"""
    Fill <: Imputor

Fills in the missing data with a specific value.

# Fields
* `value::Any`: A scalar missing value or a function that returns the a scalar if
    passed the data with missing data removed (e.g, `mean`)
"""
struct Fill{T} <: Imputor
    value::T
end

"""
    Fill() -> Fill

By default `Fill()` will use the mean of the existing values as the fill value.
"""
Fill() = Fill(mean)

"""
    impute!(imp::Fill, ctx::Context, data::AbstractVector)

Computes the fill value if `imp.value` is a `Function` (i.e., `imp.value(drop(copy(data)))`)
and replaces all missing values in the `data` with that value.
"""
function impute!(imp::Fill, ctx::Context, data::AbstractVector)
    fill_val = if isa(imp.value, Function)
        imp.value(Iterators.drop(copy(data)))
    else
        imp.value
    end

    for i in 1:length(data)
        if ismissing(ctx, data[i])
            data[i] = fill_val
        end
    end

    return data
end

ext: jl


text:
## An example showing how the QtWebKit functionality can be used within julia
## This sets up two web views, one holding an editor and one for output
## The output command can be clicked on which calls back into julia to update the
## code in the editor.

module SG
using PySide

type SimpleGui
    block
    ace
    output
    function SimpleGui(parent)
        parent = PySide.project(parent)
        ace = PySide.QtWebKit[:QWebView](parent)
        f = "file:///" * Pkg.dir("PySide", "tpl", "ace.html")
        qinvoke(ace, :load, QtCore[:QUrl](f))

        output =  PySide.QtWebKit[:QWebView](parent)
        f = "file:///" * Pkg.dir("PySide", "tpl", "output.html")
        qinvoke(output, :load, QtCore[:QUrl](f))

        sp = Splitter("Vertical", parent)
        addWidget(sp, ace)
        addWidget(sp, output)
        new(sp, ace, output)
        end
end

function call_ace(x::SimpleGui, cmd, args)
    page = x.ace[:page]()
    frame = page[:currentFrame]()
    el = frame[:findFirstElement]("#editor")
    el[:evaluateJavaScript]("ace.edit('editor').$cmd($args)")
end

getValue(x::SimpleGui) = call_ace(x, "getValue", "")
setValue(x::SimpleGui, value::String) = call_ace(x, "setValue", "'$value', 0")
getSelection(x::SimpleGui) = call_ace(x, "getSelection", "")

function prependOutput(x::SimpleGui, txt)
    outel = x.output[:page]()[:currentFrame]()[:findFirstElement]("#output")
    outel[:prependInside](txt)
end
function execJsOutput(x::SimpleGui, cmd)
    if cmd == nothing return() end
    outel = x.output[:page]()[:currentFrame]()[:findFirstElement]("#output")
    outel[:evaluateJavaScript](cmd)
end

end

 ## using QtWebKit to make a simple interface
 using Mustache,  JSON
# using GoogleCharts
using PySide
reload(Pkg.dir("PySide", "examples", "Evaluate.jl"))
using Evaluate
using SG

##
##################################################
prepare_out(x::Any) = ("
$(x)
", nothing)
prepare_out(x::Nothing) = ("", nothing)
prepare_out(x::Function) = prepare_out(nothing)

function eval_parsed_block(expr)
    ## warp Evaluate.exec_cmd in try/catch
    ## return ???
    try
        Evaluate.exec_cmd(expr)
    catch e
        ("Error: ($e)...", [])
    end
end

function eval_block()
    txt = SG.getValue(sg)
    println(txt)
    for (cmd, expr) in Evaluate.parseit(txt)
        result, output = eval_parsed_block(expr)
        println(result)
        (html, js) = prepare_out(result) # what to do with output (printed to STDOUT)
        SG.prependOutput(sg, html)
        SG.execJsOutput(sg, js)
        ## extra bit calls back into julia on a click event
        SG.prependOutput(sg, "
$cmd
")
    end
    SG.setValue(sg, "")
end

w = MainWindow()
setWindowTitle(w, "Simple GUI")
sg = SG.SimpleGui(w)
setCentralWidget(w, sg.block)

## ## Define some actions
## evaluate ... Qt.QKeySequence("Ctrl+X, Ctrl+C");
## reload ...
evaluate_action = Qt[:QAction](PySide.project(w))
qinvoke(evaluate_action, :setText, "evaluate")
qinvoke(evaluate_action, :setShortcut, Qt[:QKeySequence]("Open")) ## This fails!!
qinvoke(evaluate_action, :setStatusTip, "Evaluate buffer contents")
qconnect(evaluate_action, :triggered, eval_block)

tb = Qt[:QToolBar](PySide.project(w))
qinvoke(tb, :addAction, evaluate_action)
qinvoke(w, :addToolBar, tb)

## ## Integration of JavaScript -> julia
qnew_class("Julia", "QtCore.QObject", meths=[Dict([(:meth_dfn, "evaluate = QtCore.Signal(str)")])])
julia = qnew_class_instance("Julia")
qconnect(julia, :evaluate) do cmd
  SG.setValue(sg, cmd)
end

## How to call JavaScript from within Qt
## This sets up a function `eval_julia` for evaluating the contents of a div tag
## see the formatting of the output above
frame = sg.output[:page]()[:currentFrame]()
frame[:addToJavaScriptWindowObject]("julia", julia)
frame[:evaluateJavaScript]("window.eval_julia = function(self) {julia.evaluate(\$(self).text());}")

raise(w)


text:
toStrTuple(bpe::GenericBPE{T}, x::AbstractString) where T = Tuple(T.(merges(bpe, x)))

struct Statistic{B<:GenericBPE}
  bpe::B
  vocab::Dict{Tuple, Int}
  vkeys::Vector{Tuple}
  pair_freq::Dict{Pair{String, String}, Int}
  pair_index::Dict{Pair{String, String}, Vector{Int}}
end

Statistic(bper::BPELearner) = Statistic(bper.bpe, bper.vocabs)
function Statistic(bpe::GenericBPE, v::Dict{String, Int})
  vocab = Dict{Tuple, Int}((toStrTuple(bpe, k), v) for (k, v) ∈ pairs(v))
  vkeys = collect(keys(vocab))
  pair_freq = Dict{Pair{String, String}, Int}()
  pair_index = Dict{Pair{String, String}, Vector{Int}}()

  for (i, tp) ∈ enumerate(vkeys)
    for p ∈ bi_pairs(tp)
      pair_freq[p] = get(pair_freq, p, 0) + vocab[tp]
      indices = get!(pair_index, p, Int[])
      push!(indices, i)
    end
  end

  Statistic(bpe, vocab, vkeys, pair_freq, pair_index)
end

"find adjacent element. return a list of Pair"
function bi_pairs(stp)::Vector{Pair{String, String}}
    init, rstp = Iterators.peel(stp)
    map((x)->x[1]=>x[2], zip(stp, rstp))
end

"if only length one, return []"
bi_pairs(stp::Tuple{String}) = Vector{Pair{String, String}}()

"most frequently pair"
most_freq(stats::Statistic) = argmax(stats.pair_freq)

"get pair freq"
get_freq(stats::Statistic, p::Pair) = get(stats.pair_freq, p, 0)

"get word freq"
get_freq(stats::Statistic, str::String) = get_freq(stats, toStrTuple(stats.bpe, str))
get_freq(stats::Statistic, tp) = get(stats.vocab, tp, 0)

"merge the gived pair and update Statistic"
function merge_pair!(stats::Statistic, pair::Pair{String, String})
    bpi = collect(stats.pair_index[pair])
    sort!(bpi); unique!(bpi)
    for index ∈ bpi
        wtp = stats.vkeys[index]
        freq = stats.vocab[wtp]
        delete!(stats.vocab, wtp)

        bps = bi_pairs(wtp)
        nwtp, indices = merged_pairs(bps, pair)
        adjacents = adjacent_pairs(bps, indices)
        mps = setdiff!(bi_pairs(nwtp), bps)

        stats.vkeys[index] = nwtp
        stats.vocab[nwtp] = freq

        update!(stats, index, nwtp, mps)
        update!(stats, index, adjacents)
    end

    delete!(stats.pair_freq, pair)
    delete!(stats.pair_index, pair)
    stats
end

"update Statistic of the new merged pair"
function update!(stats::Statistic, index::Int, nwtp::Tuple, mps::Vector)
    for mp ∈ mps
        stats.pair_freq[mp] = get_freq(stats, mp) + stats.vocab[nwtp]
        indices = get!(stats.pair_index, mp, Int[])
        push!(indices, index)
    end
end

"update Statistic of the removed adjacent pairs"
function update!(stats::Statistic, index::Int, adjacents::Vector)
    for ap ∈ adjacents
        stats.pair_freq[ap] -= stats.vocab[stats.vkeys[index]]
        dindex = findfirst(isequal(index), stats.pair_index[ap])
        dindex !== nothing && deleteat!(stats.pair_index[ap], dindex)
    end
    stats
end

"update a tuple"
function update!(stats::Statistic, index::Int, nwtp::Tuple)
    nbps = bi_pairs(nwtp)
    for np ∈ nbps
        stats.pair_freq[np] = get_freq(stats, np) + stats.vocab[nwtp]
        indices = get!(stats.pair_index, np, Int[])
        push!(indices, index)
    end
    stats
end

"update a word"
function update!(stats::Statistic, str::String; freq::Int = 1)
    stp = toStrTuple(stats.bpe, str)
    if haskey(stats.vocab, stp)
        stats.vocab[stp] += freq
        for sbp ∈ bi_pairs(stp)
            stats.pair_freq[sbp] += freq
        end
    else
        stats.vocab[stp] = get_freq(stats.vocab, stp) + freq
        push!(stats.vkeys, stp)
        update!(stats, length(stats.vkeys), stp)
    end
    stats
end

"update a frequency vocab"
function update!(stats::Statistic, vocab::Dict{String, Int})
    for (word, freq) ∈ vocab
        update!(stats, word; freq=freq)
    end
    stats
end

"get the adjacent pair of merged pair"
adjacent_pairs(wtp::Tuple, pair::Pair{String, String}) = adjacent_pairs(bi_pairs(wtp), pair)
function adjacent_pairs(_bps::Array{Pair{String, String}}, pair::Pair{String, String})
    bps = collect(_bps)
    indice = Vector{Int}()

    merged = intern(join(pair))
    for (i, x) ∈ enumerate(bps)
        if 1<i && bps[i-1] == pair
            bps[i] = merged=>x.second
        elseif x == pair
            push!(indice, i)
        end
    end
    adjacent_pairs(_bps, indice)
end

function adjacent_pairs(bps::Array{Pair{String, String}}, indices::Vector{Int})
    adjs = filter!((i)->1<=i<=length(bps), unique!(foldl((l,r)->(push!(l, r-1);push!(l, r);push!(l, r+1);l), indices, init=Int[])))
    bps[adjs]
end

"merge the pair in the word tuple"
merged_pairs(wtp::Tuple, pair::Pair{String, String}) = merged_pairs(bi_pairs(wtp), pair)
function merged_pairs(bps::Array{Pair{String, String}}, pair::Pair{String, String})
    bps = collect(bps)
    nwss = Vector{String}()
    indice = Vector{Int}()

    merged = intern(join(pair))
    for (i, x) ∈ enumerate(bps)
        if 1<i && bps[i-1] == pair
            bps[i] = merged=>x.second
        elseif x == pair
            push!(nwss, merged)
            push!(indice, i)
            continue
        else
            push!(nwss, x.first)
        end

        if i == length(bps)
            push!(nwss, x.second)
        end
    end
    Tuple(nwss), indice
end


text:
function load_sequences( file_path )
    seqs = ASCIIString[]
    scores = Float64[]
    open( file_path ) do f
        for line in eachline( f )
            spl = split( line )
            push!( seqs, spl[1] )
            push!( scores, float64( spl[2] ) )
        end
    end
    return seqs,scores
end

max_stars_repo_name: tbepler/ets-analysis-scripts

max_issues_repo_name: tbepler/ets-analysis-scripts

max_forks_repo_name: tbepler/ets-analysis-scripts


text:
#!/usr/bin/env julia
function import_hop(fname,sw_ndegen)
    tmp=open(fname*"/irvec.txt","r") do fp
        readlines(fp)
    end
    nr=length(tmp)
    rvec=zeros(Float64,3,nr)
    for i in 1:nr
        rvec[:,i]=parse.(Float64,split(tmp[i]))
    end
    tmp=open(fname*"/ham_r.txt","r") do fp
        readlines(fp)
    end
    no=Int(sqrt(length(tmp)/nr))
    tmp2=[parse.(Float64,split(strip(l,['\n','(',')']),',')) for l in tmp]
    tmp=[t[1]+1im*t[2] for t in tmp2]
    ham_r=zeros(Complex,nr,no,no)
    for i in 1:no
        for j in 1:no
            for k in 1:nr
                ham_r[k,j,i]=tmp[i+no*(j-1)+no*no*(k-1)]
            end
        end
    end

    if sw_ndegen
        tmps=open(fname*"/ndegen.txt","r") do fp
            readlines(fp)
        end
        ndegen=parse.(Float64,tmp)
    else
        ndegen=ones(Float64,(1,nr))
    end
    return(rvec,ndegen,ham_r,no,nr)
end

function import_out(fname)
    tmp=open(fname,"r") do fp
        readlines(fp)
    end
    tmp1=[parse.(Float64,split(l)) for l in tmp]
    r0=tmp1[1][1:3]
    con=0
    for rr in tmp1
        if rr[1:3]==r0
            con+=1
        end
    end
    no=Int(sqrt(con))
    nr=div(length(tmp1),con)
    rvec=zeros(Float64,3,nr)
    for i in 1:nr
        for j in 1:3
            rvec[j,i]=tmp1[i][j]
        end
    end
    tmp=[tp[4]+1im*tp[5] for tp in tmp1]
    ham_r=zeros(Complex,nr,no,no)
    for i in 1:no
        for j in 1:no
            for k in 1:nr
                ham_r[k,j,i]=tmp[k+nr*(j-1)+no*nr*(i-1)]
            end
        end
    end
    ndegen=ones(Float64,1,nr)
    return(rvec,ndegen,ham_r,no,nr)
end

function import_hr(fname)
    tmp=open(fname*"_hr.dat","r") do fp
        readlines(fp)
    end
    no,nr=int(tmp[2]),int(tmp[3])
    return(no,nr)
end

function import_Hopping()
    tmp=open("Hopping.dat","r") do fp
        readlines(fp)
    end
    tmp1=[split(l) for l in tmp[2:5]]
    axis=[[parse(Float64,l) for l in tp] for tp in tmp1]
    tmp1=split(tmp[5])
    no=parse(Int,tmp1[1])
    nr=parse(Int,tmp1[2])
    nor=parse(Int,tmp1[3])
    tmp1=[tmp[no*no*i][:3] for i in 1:nr]
    rvec=[[parse(Float64,l) for l in tp] for tp in tmp1]
    tmp1=[split(l) for l in tmp[(7+no):end]]
    ham_r=reshape([parse(Float64,tp[9])+1im*parse(Float64,tp[10]) for tp in tmp1],nr,no,no)
    return(rvec,ndegen,ham_r,no,nr,axis)
end


text:
using Images, GtkReactive, ImageView
using TestImages
using Gtk.ShortNames  # for signal_connect
using Random  # for the init game
import Base.show # Used for adding struct specific prints

# used example from    https://discourse.julialang.org/t/imageview-gtk-unastable-results/25267

mutable struct GameState
    img::Array{T,2} where T <:Colorant
    shuf_img::Array{T,2} where T <:Colorant
    guidict::Dict{String,Any}
    isActive::Bool
    isFirstClick::Bool
    nof_blocks::Int
    block_size::Int

    box_handler::Union{Nothing, ImageView.AnnotationHandle{ImageView.AnchoredAnnotation{AnnotationBox}} }
    blocks::Union{Nothing, Array{Tuple{UnitRange{Int64},UnitRange{Int64}},2}}
    swaps::Int64

    function GameState(img::Array{T,2}, nof_blocks::Int) where T<:Colorant
        shuf_img, blocks, block_size = initGame(img, nof_blocks)
        guidict = imshow(shuf_img)
        new(img, shuf_img, guidict, true, true, nof_blocks, block_size, nothing,blocks,0)
    end

    function GameState(img::Array{T,2}) where T<:Colorant
        #default if 4 block per row
        shuf_img, blocks, block_size = initGame(img, 4)
        guidict = imshow(shuf_img)
        new(img, shuf_img, guidict, true, true, 4, block_size, nothing,blocks,0)
    end
end
"""
Dedicated show function for the GameState struct
"""
function show(io::IO, ::MIME"text/plain", game::GameState)
    active_str = game.isActive ? "Active" : "Not Active"
    window_name = get_gtk_property(game.guidict["gui"]["window"], :name, String)
    ann_cnt = length(game.guidict["annotations"].value)
    print(io, "Game State: $(active_str).")
          #In window '$(window_name)'' has $(ann_cnt) annotations and $(game.swaps) swaps (nof_blocks = $(game.nof_blocks))")
end

"""
Override the print version as well
"""
function show(io::IO, game::GameState)
    active_str = game.isActive ? "Active" : "Not Active"
    window_name = get_gtk_property(game.guidict["gui"]["window"], :name, String)
    ann_cnt = length(game.guidict["annotations"].value)
    print(io, "Game State: $(active_str).",
          "In window '$(window_name)'' has $(ann_cnt) annotations and $(game.swaps) swaps (nof_blocks = $(game.nof_blocks))")
end

"""
Shuffles the input image to an unorderded versoin
"""
function initGame(img::Array{T,2}, nof_blocks::Int) where T<:Colorant
    block_size = Int(size(img)[1] / nof_blocks)
    # i* block_size + 1 : (i+1) * block_size
    row_range = [i* block_size + 1 : (i+1) * block_size for i in 0:nof_blocks-1]

    blocks = [(x,y) for x in row_range for y in row_range]
    blocks = reshape(blocks, nof_blocks, nof_blocks)  # orginze in a matrix for bebugging

    new = zeros(RGB, size(img))
    blocks_new = shuffle(blocks)
    for b in 1:length(blocks)
        new[blocks_new[b][1], blocks_new[b][2]] = img[blocks[b][1], blocks[b][2]]
    end
    return new, blocks, block_size
end

function close_window_signal(c::Base.GenericCondition{Base.AlwaysLockedST}, game::GameState)
    println("Window closed - finishing game")
    game.isActive = false
    notify(c)
end

"""
return box coordinates
if mode == "box" - coordinage in left, top, right. bottom
if mode == "ind" - indeceis matching the box index
"""
function get_box_coords(coord::Tuple{AbstractFloat,AbstractFloat}; box_size::Int=32, mode::String="box")
    hor_offset = rem(coord[1], box_size)
    ver_offset = rem(coord[2], box_size)

    row_ind = Int(div(coord[1], box_size))
    col_ind = Int(div(coord[2], box_size))

    left = coord[1] - hor_offset + 1
    top = coord[2] - ver_offset + 1

    right = coord[1] - hor_offset + box_size
    bottom = coord[2] - ver_offset + box_size

    if mode == "box"
        println("box: left=$left top=$top right=$right bottom=$bottom")
        return Int.(round.([left,top,right,bottom]))
    elseif mode == "ind"
        println("box index row=$row_ind col=$col_ind")
        return row_ind,col_ind
    else
        nothing
    end
end

function mark_box(left,top,right,bottom; color::RGB=RGB(0,0,1))
    global guidict
    box_ann = annotate!(guidict, AnnotationBox(left,top,right,bottom, linewidth=2.5, color=color))
end

function swap_boxes(game_state::GameState, boxA, boxB)
    img = game_state.shuf_img
    boxA = Int.(round.(boxA))
    boxB = Int.(round.(boxB))
    game_state.box_handler.ann.data.left
    tmp = copy(img[boxA[2]:boxA[4], boxA[1]:boxA[3]])
    img[boxA[2]:boxA[4], boxA[1]:boxA[3]] = img[boxB[2]:boxB[4], boxB[1]:boxB[3]]
    img[boxB[2]:boxB[4], boxB[1]:boxB[3]] = tmp
    img
end

# ===============
# Main game Flow
# ===============

nof_blocks = 2
img = testimage("lena")

# open window and show the image
game_state = GameState(img, nof_blocks)
guidict = game_state.guidict

# Create a condition object
c = Condition()

# shorthandes for some of the objects
win = game_state.guidict["gui"]["window"]
can = game_state.guidict["gui"]["canvas"]

gui = game_state.guidict["gui"]

# Notify the condition object when the window closes
signal_connect(close_window_signal, win, :destroy)

# controlling the button press signal - https://juliagizmos.github.io/GtkReactive.jl/stable/drawing.html
# NOTE: using a simple map will work in interacive way , but when used in a script we will need to prevent it from getting
#       collected by the garbage collector - the way is to assign it to a varable and use prevent gc
sig_click = map(can.mouse.buttonpress) do btn  # Reactive.Signal{GtkReactive.MouseButton{GtkReactive.UserUnit}}
    global game_state
    x_curr = Float64(btn.position.x)  #row
    y_curr = Float64(btn.position.y)  # col
    println("Found a button clicked @ [x=$(round(x_curr, digits=2)), y=$(round(y_curr, digits=2))]")
    println(game_state)

    # invalid points occur at startup
    if x_curr > 0 && y_curr > 0
        if game_state.isFirstClick
            println("first click")  # TODO: add hihlight of the block
            game_state.isFirstClick = false
            cc = get_box_coords((x_curr, y_curr), box_size=game_state.block_size)
            #box_ann = mark_box(cc[1], cc[2], cc[3], cc[4])
            game_state.box_handler = annotate!(game_state.guidict, AnnotationBox(cc[1], cc[2], cc[3], cc[4], linewidth=2.5, color=colorant"blue"))

            println("cc meta data:",typeof(cc), cc)
        else
            println("second click - need to swap") # TODO: add swap functionality
            game_state.isFirstClick = true
            game_state.swaps += 1
            left,top,right,bottom = get_box_coords((x_curr, y_curr), box_size=game_state.block_size)
            # FIXME try showing a very short period of selected box - and then sleep and remove
            #       the refresh does not work

            swap_boxes(game_state,
                       [game_state.box_handler.ann.data.left, game_state.box_handler.ann.data.top,
                        game_state.box_handler.ann.data.right, game_state.box_handler.ann.data.bottom],
                       [left,top,right,bottom])

            delete!(game_state.guidict, game_state.box_handler)
        end
    end
    if game_state.img == game_state.shuf_img
        println("Game WON!!!  finishing game")
        notify(c)
    end
end

GtkReactive.gc_preserve(win, sig_click)  # prevetns the mapping from being garbage collected

# Wait for the notification before proceeding ...
wait(c)


text:
@testset "$TEST $G" begin

e1 = E(1,2)
re1 = E(2,1)
@test e1.src == src(e1) == 1
@test e1.dst == dst(e1) == 2
@test reverse(e1) == re1
@test sort(e1) == e1
@test sort(re1) == e1

e = first(edges(G(10,20)))
@test E == typeof(e)
if E <: AIndexedEdge
    @test idx(e) > 0
end

e = first(edges(DG(10,20)))
@test E == typeof(e)
if E <: AIndexedEdge
    @test idx(e) > 0
end

end # testset


text:
import accrue.cryptoerase.runtime.Condition;

class Foo {
    public Foo() {}
}

class C {
  void m() {
      int{H} x = 42;
      if (x > 0) {
	  Foo{L} f = new Foo();
      }
  }

  public static void main(String[] args) {
    new C().m();
  }
}

max_stars_repo_name: HarvardPL/cryptoerase

max_issues_repo_name: HarvardPL/cryptoerase

max_forks_repo_name: HarvardPL/cryptoerase


text:
struct TaxonNode{T <: AbstractTaxon}
    taxon::T
    rank::TaxonomicRank
    parent::Union{TaxonNode, Nothing}
    children::Vector{TaxonNode}
end

"""
    TaxonNode(nctx::NCBITaxon)

    Constructor for TaxonNode with no parent (meaning it is the root of a tree).
"""
TaxonNode(nctx::AbstractTaxon, rank::TaxonomicRank)  = TaxonNode(nctx, rank, nothing, TaxonNode[])

"""
    TaxonNode(nctx::NCBITaxon,  parent::TaxonNode)

    Constructor for TaxonNode with a parent
"""
TaxonNode(nctx::AbstractTaxon, rank::TaxonomicRank, parent::TaxonNode) = TaxonNode(nctx, rank, parent, TaxonNode[])

taxon(node::TaxonNode{NCBITaxonWrapper}) = node.taxon.val

struct FlexibleRankClassificationTree <: AbstractTaxonomy
    root::TaxonNode
end
struct FixedRankClassificationTree <: AbstractTaxonomy
    root::TaxonNode
end

function AbstractTrees.children(node::TaxonNode)
    if isdefined(node, :children)
        return node.children
    end
    return []
end

AbstractTrees.printnode(io::IO, node::TaxonNode) = print(io, Taxonomy.taxon(node))

function _addchild!(parent::TaxonNode, child::TaxonNode )
    push!(parent.children, child)
end

function _addnode!(newtaxon::AbstractTaxon, rank::TaxonomicRank, parent::TaxonNode)
    newnode = TaxonNode(newtaxon, rank, parent)
    _addchild!(parent, newnode)

    return newnode
end


text:
function parseConvergence(fileName::String)

  file = open(fileName)
  lines = readlines(file)
  close(file)

  frequencyMask = occursin.(
    "HelmholtzSolve:  Frequency (Hz):",
    lines
  )

  itStartMask = occursin.(
    "IterSolver: Calling complex iterative solver",
    lines
  )

  itEndMask = occursin.(
    "ComputeChange: NS (ITER=1)",
    lines
  ) .| occursin.(
    "NUMERICAL ERROR:: IterSolve: Too many iterations was needed.",
    lines
  )

  lineNumber = 1:length(lines)
  nSolutions = length(lines[frequencyMask])

  data = Dict()

  for n in 1:nSolutions

    freq = parse(Float64, split(lines[frequencyMask][n])[end])

    convLines = split.(
      lines[
        (lineNumber[itStartMask][n] + 1):(lineNumber[itEndMask][n] - 1)
      ]
    )

    it = [parse.(Float64, convLines[n])[1] for n in 1:length(convLines)]
    v = [parse.(Float64, convLines[n])[2] for n in 1:length(convLines)]

    convData = [it v]

    data[freq] = convData

  end

  return data

end

max_stars_repo_name: CrocoDuckoDucks/Acoustics-Models

max_issues_repo_name: CrocoDuckoDucks/Acoustics-Models

max_forks_repo_name: CrocoDuckoDucks/Acoustics-Models


text:
function CorrelationId()
    CorrelationId(UInt32(0), CorrelationValue())
end

function CorrelationValue(val::UInt=UInt(0))
    CorrelationValue(val, (UInt(0), UInt(0), UInt(0), UInt(0)), UInt(0))
end

encode_value(val::UInt) = CorrelationValue(val)

Base.:(==)(c1::CorrelationId, c2::CorrelationId) = c1.header == c2.header && c1.value == c2.value
Base.hash(c::CorrelationId) = hash(c.header) + hash(c.value)

struct CorrelationIdHeader
    struct_size::UInt32
    correlation_type::CorrelationType
    class_id::UInt32
    reserved::UInt32

    function CorrelationIdHeader(struct_size::UInt32, correlation_type::CorrelationType, class_id::UInt32, reserved::UInt32)
        @assert struct_size <= 0xff
        @assert class_id <= 0xffff
        @assert reserved <= 0xf

        return new(struct_size, correlation_type, class_id, reserved)
    end
end

function CorrelationIdHeader(c::CorrelationId)
    CorrelationIdHeader(c.header)
end

function CorrelationIdHeader(header::UInt32)
    CorrelationIdHeader(
            # [ - 8 - ][ - 4 - ][ - 16 - ][ - 4 - ] >> 24
            # [ - zeros - ][ - 8 - ]
            header >> 24,

            # [ - 8 - ][ - 4 - ][ - 16 - ][ - 4 - ] << 8
            # [ - 4 - ][ - 16 - ][ - 4 - ][ - 8 zeros - ] >> 28
            # [ - zeros - ][ - 4 - ]
            CorrelationType(header << 8 >> 28),

            # [ - 8 - ][ - 4 - ][ - 16 - ][ - 4 - ] << 12
            # [ - 16 - ][ - 4 - ][ - 12 zeros - ] >> 16
            # [ - zeros - ][ - 16 - ]
            header << 12 >> 16,

            # [ - 8 - ][ - 4 - ][ - 16 - ][ - 4 - ] << 28
            # [ - 4 - ][ - 28 - ] >> 28
            # [ - zeros - ][ - 4 - ]
            header << 28 >> 28
        )
end

function encode_header(header::CorrelationIdHeader) :: UInt32
    (((((header.struct_size << 4) + UInt32(header.correlation_type)) << 16) + header.class_id) << 4) + header.reserved
end

function Base.show(io::IO, c::CorrelationId; debug::Bool=false)
    if debug
        header = CorrelationIdHeader(c)
        println(io, "CorrelationId")
        println(io, "  size $(header.struct_size)")
        println(io, "  value type $(header.correlation_type)")
        println(io, "  class id $(header.class_id)")
        println(io, "  reserved bytes $(header.reserved)")
        print(io, "  value $(c.value)")
    else
        print(io, "CorrelationId($(c.value.ptr))")
    end
end

const DEFAULT_CLASS_ID = UInt32(0)

function CorrelationId(value::UInt64; class_id::UInt32=DEFAULT_CLASS_ID)
    header = CorrelationIdHeader(UInt32(sizeof(CorrelationId)), BLPAPI_CORRELATION_TYPE_INT, class_id, UInt32(0))
    return CorrelationId(header, value)
end

CorrelationId(value::Integer; class_id::UInt32=DEFAULT_CLASS_ID) = CorrelationId(UInt64(value), class_id=class_id)
CorrelationId(header::CorrelationIdHeader, value::UInt64) = CorrelationId(encode_header(header), encode_value(value))

max_stars_repo_licenses: OpenSSL MIT

max_issues_repo_licenses: OpenSSL MIT

max_forks_repo_licenses: OpenSSL MIT


text:
include("../src/HelperFunctions.jl")
include("../src/InitialConditions.jl")
include("../src/KalmanFilter.jl")
include("../src/EM.jl")
include("../src/EstimateDFM.jl")
include("../src/Update.jl")

using Suppressor

using DataFrames
using Statistics
using Dates
using LinearAlgebra
using SparseArrays
using CSV

# data set up
sample_data = CSV.File("test_data.csv") |> DataFrame!
y_tmp = sample_data[!, Not(:date)]
dates = sample_data[!, :date]
q = zeros(4)
p = 1
monthly_quarterly_array = gen_monthly_quarterly(dates, y_tmp)
R_mat = [2 -1 0 0 0; 3 0 -1 0 0; 2 0 0 -1 0; 1 0 0 0 -1]
blocks = DataFrame(a=ones(size(y_tmp)[2]), b=ones(size(y_tmp)[2]))
init_conds = initialize_conditions(y_tmp; dates=dates, p=1, blocks=blocks, R_mat=R_mat)
A = init_conds[:A]; C = init_conds[:C]; Q = init_conds[:Q]; R = init_conds[:R]; Z0 = init_conds[:Z0]; V0 = init_conds[:V0]
y_est = y_tmp[[sum(.!ismissing.(Array(x))) > 0 for x in eachrow(y_tmp)], :] |> Array |> transpose
nM = sum(.!monthly_quarterly_array)

# creating an artificial lagged dataset
data_lag = create_lag(sample_data, 1)
old_y = copy(data_lag)
new_y = copy(sample_data)
target_variable = Symbol("x_world.sa")
target_period = Dates.Date(2020,6,1)

@testset "Helper functions" begin

	# testing standardize function
	test_data = DataFrame(Dict(:a => ["a", "a", "a"], :b => [1,2,3]))
	@test standardize(test_data) == DataFrame(Dict(:a => ["a", "a", "a"], :b => [-1.0,0.0,1.0]))

	#testing spline_fill function
	@test isequal(spline_fill([1,2,missing,4]), [1.0,2.0,3.0,4.0])
	@test isequal(spline_fill([1,2,missing,5,missing]), [1.0,2.0,3.375,5.0,missing])

	# testing digital_filter function
	@test (sum((digital_filter([missing, 1,2,5,3,missing]) .- [2.16, 1.0, 2.0, 5.0, 3.0, 2.33]))) < 0.01

	# testing fill_na function
	@test ((DataFrame(Dict(:x=>[missing, 1, 2, missing, 3.5], :y=>[24, 42, 76, 89.0, missing])) |> fill_na)[:output] .- DataFrame(Dict(:x=>[2.07725, 1.0, 2.0, 2.8125, 3.5], :y=>[24.0, 42.0, 76.0, 89.0, 61.9689])) |> x-> sum(x.x) + sum(x.y)) < 0.001
	@test (DataFrame(Dict(:x=>[missing, 1, 2, missing, 3.5], :y=>[24, 42, 76, 89.0, missing])) |> fill_na)[:na_indices] == DataFrame(Dict(:x=>[1,0,0,1,0], :y=>[0,0,0,0,1]))

	# testing gen_monthly_quarterly_matrix function
	test_dates = [Dates.Date(2020,i,1) for i in 1:12]
	test_data = DataFrame(
		Dict(
			:a => [0 for i in 1:12],
			:b => [missing,missing,0,missing,missing,0,missing,missing,0,missing,missing,0]
		)
	)
	@test gen_monthly_quarterly(test_dates, test_data) == [0,1]

	# testing date_col_name function
	test_data = DataFrame(a=[Dates.Date(2020,1,1), Dates.Date(2020,2,1)], b=[1,2])
	@test date_col_name(test_data) == :a
	error_out = @capture_out date_col_name(test_data[!, Not(:a)])
	@test error_out == "No column of type Dates.Date\n"
end

@testset "Initial conditions" begin
	Y = DataFrame(Dict(:a=>[1,2,4,3,2], :b=>[2,1,3,1,5]))
	dates = [Dates.Date(2020,1,1), Dates.Date(2020,2,1), Dates.Date(2020,3,1), Dates.Date(2020,4,1), Dates.Date(2020,5,1)]
	blocks = DataFrame(Dict(:a=>[1.0, 1.0], :b=>[1.0, 1.0]))
	R_mat = [2 -1 0 0 0; 3 0 -1 0 0; 2 0 0 -1 0; 1 0 0 0 -1]
	tmp = initialize_conditions(Y; dates=dates, p=p, blocks=blocks, R_mat=R_mat)

	@test sum(tmp[:A]) - 2.438 < 0.001
	@test sum(tmp[:C]) - 4.2649 < 0.001
	@test isequal(sum(tmp[:Q]), NaN)
	@test tmp[:R] == [0.0001 0.0; 0.0 0.0001]
	@test tmp[:Z0] == zeros(12)
	@test isequal(sum(tmp[:V0]), NaN)
end

@testset "Kalman filter" begin
	output = kalman_filter(y_est, A, C, Q, R, Z0, V0)

	@test sum(output[:Zsmooth]) ≈ 17.108467354695524
	@test sum(output[:Vsmooth]) ≈ 808.0226237078814
	@test sum(output[:VVsmooth]) ≈ 232.87470774743034
	@test sum(output[:loglik]) ≈ 541.3938089668698

	output_dfm = estimate_dfm(sample_data; blocks=blocks, p=p, max_iter=10, threshold=1e-5)
	constparams = kalman_filter_constparams(y_tmp; output_dfm=output_dfm, lag=0)

	@test sum(constparams[:Plag][1]) ≈ 79.91396625742085
	@test sum(constparams[:X_smooth]) ≈ 14.352390345012175
	@test sum(constparams[:Vsmooth]) ≈ 80.02804204917477
	@test sum(constparams[:F]) ≈ 109.36315084230503

	constparams = kalman_filter_constparams(y_tmp; output_dfm=output_dfm, lag=2)
	@test sum(skipmissing(constparams[:Plag][2])) ≈ 22.067361604719792
	@test sum(skipmissing(constparams[:Plag][3])) ≈ 23.35769893782239
end

@testset "EM" begin
	em_output = EM_step(y_est; A=A, C=C, Q=Q, R=R, Z0=Z0, V0=V0, p=p, blocks=blocks, R_mat=R_mat, q=q, nM=nM, monthly_quarterly_array=monthly_quarterly_array)

	@test sum(em_output[:A_new]) ≈ 11.552663255407325
	@test sum(em_output[:C_new]) ≈ 23.45337243218924
	@test sum(em_output[:Q_new]) ≈ 1.594071735564734
	@test sum(em_output[:Z0]) ≈ 0.031060939148336537
	@test sum(em_output[:V0]) ≈ 48.99573305597052
	@test sum(em_output[:loglik]) ≈ 541.3938089668698

	loglik = 10; prev_loglik = 1e-6;
	converged, decrease = values(EM_convergence(loglik, prev_loglik, 1e-5))
	@test converged == 0
	@test decrease == 0

	loglik = 1e-8; prev_loglik = 5;
	converged, decrease = values(EM_convergence(loglik, prev_loglik, 1e-1))
	@test converged == 0
	@test decrease == 1

	loglik = 2; prev_loglik = 5;
	converged, decrease = values(EM_convergence(loglik, prev_loglik, 1))
	@test converged == 1
	@test decrease == 1
end

@testset "EstimateDFM" begin
	output = estimate_dfm(sample_data; blocks=blocks, p=p, max_iter=10, threshold=1e-5)

	@test sum(output[:Xsmooth_std]) ≈ 10.808521737527098
	@test sum(output[:C]) ≈ 23.099766843345556
	@test sum(output[:R]) ≈ 0.0011
	@test sum(output[:A]) ≈ 12.745531615217024
	@test sum(output[:Q]) ≈ 0.10217159203541483
	@test sum(output[:means]) ≈ 0.060732284509809845
	@test sum(output[:stds]) ≈ 0.5829554656784255
	@test sum(output[:Z0]) ≈ 0.7589260732889963
	@test sum(output[:V0]) ≈ 48.99573305597052
	@test sum(output[:loglik]) ≈ 6544.007956191773
	@test sum(output[:LL]) ≈ 54935.04005844516

	predictions = predict_dfm(sample_data; output_dfm=output, months_ahead=3, lag=0)
	@test sum([sum(i) for i in eachcol(predictions)[2:end]]) ≈ 14.124786767120007
end

@testset "Update" begin
	output_dfm = estimate_dfm(sample_data; blocks=blocks, p=p, max_iter=10, threshold=1e-5)
	news = news_dfm(;old_y=data_lag, new_y=sample_data, output_dfm=output_dfm, target_variable=Symbol("x_world.sa"), target_period=Dates.Date(2020,6,1))

	@test sum(news[:actual]) ≈ -0.6080187052009832
	@test sum(news[:forecast]) ≈ -0.8269307754720284
	@test sum(news[:weight]) ≈ 1.9326984705441284
	@test sum(news[:innov]) ≈ 15.778245847006348
	@test sum(news[:row_miss]) ≈ 2427
	@test sum(news[:col_miss]) ≈ 66
	@test sum(news[:singlenews]) ≈ 0.04320382154267746
	@test sum(news[:y_old]) ≈ -0.3871829288269627
	@test sum(news[:y_new]) ≈ -0.33385843789784786

	updated_nowcast = gen_news(;old_y=create_lag(new_y, 1), new_y=new_y, output_dfm=output_dfm, target_variable=Symbol("x_world.sa"), target_period=Dates.Date(2020,6,1))[:news_table]

	@test sum(skipmissing(updated_nowcast[!, :forecast])) ≈ -82.69307754720283
	@test sum(skipmissing(updated_nowcast[!, :actual])) ≈ -60.80187052009832
	@test sum(skipmissing(updated_nowcast[!, :weight])) ≈ 1.9326984705441284
	@test sum(skipmissing(updated_nowcast[!, :impact_releases])) ≈ 4.3203821542677465
	@test sum(skipmissing(updated_nowcast[!, :impact_total])) ≈ 4.3203821542677465
	@test sum(skipmissing(updated_nowcast[!, :data_release])) ≈ 2
end

@testset "Full functionality" begin
	calc  = estimate_dfm(sample_data; blocks=blocks, p=1, max_iter=100, threshold=1e-5)
	export_dfm(output_dfm=calc, out_path="output")
	imported  = import_dfm(path="output")

	predcalc = predict_dfm(sample_data; output_dfm=calc, months_ahead=3, lag=0)
	@test sum([sum(i) for i in eachcol(predcalc)[2:end]]) ≈ 15.531695602544705

	predimported = predict_dfm(sample_data; output_dfm=imported, months_ahead=3, lag=0)
	@test sum([sum(i) for i in eachcol(predimported)[2:end]]) ≈ 15.531695602544705

	news = gen_news(old_y=sample_data, new_y=sample_data, output_dfm=imported, target_variable=Symbol("x_world.sa"), target_period=Dates.Date(2020,6,1))
	@test sum([sum(i) for i in eachcol(news[:news_table])[2:end]]) ≈ -136.49942413496785
end


text:
# Sufficient reduction projections

"""
srproj calculates the MNIR Sufficient Reduction projection from text counts on
to the attribute dimensions of interest (covars in mnlm). In particular, for
counts C, with row sums m, and mnlm coefficients φ_j corresponding to attribute
j, z_j = C'φ_j/m is the SR projection in the direction of j.
The MNIR paper explains how V=[v_1 ... v_K],
your original covariates/attributes, are independent of text counts C given SR
projections Z=[z_1 ... z_K].
dir == nothing returns projections in all directions.
"""
srproj(m::DMRPaths, counts, dir::Union{Nothing,Int}=nothing; focusj=axes(counts,2), select=defsegselect) =
  srproj(coef(m, select), counts, dir; intercept=hasintercept(m), focusj=focusj)

srproj(m::DMRCoefs, counts, dir::Union{Nothing,Int}=nothing; focusj=axes(counts,2)) =
  srproj(coef(m), counts, dir; intercept=hasintercept(m), focusj=focusj)

"""
srproj calculates the MNIR Sufficient Reduction projection from text counts on
to the attribute dimensions of interest (covars in mnlm). In particular, for
counts C, with row sums m, and mnlm coefficients φ_j corresponding to attribute
j, z_j = C'φ_j/m is the SR projection in the direction of j.
The MNIR paper explains how V=[v_1 ... v_K],
your original covariates/attributes, are independent of text counts C given SR
projections Z=[z_1 ... z_K].
dir == nothing returns projections in all directions.
"""
function srproj(coefs::AbstractMatrix{T}, counts, dir::Union{Nothing,Int}=nothing; intercept=true, focusj=axes(counts,2)) where T
   ixoffset = intercept ? 1 : 0 # omitting the intercept
   if dir==nothing
     Φ=coefs[ixoffset+1:end,focusj]' # all directions
   else
     Φ=coefs[ixoffset+dir:ixoffset+dir,focusj]' # keep only desired directions
   end
   m = sum(counts, dims=2) # total counts per observation
   z = counts[:,focusj]*Φ ./ (m+(m.==0)) # scale to get frequencies
   [z m] # m is part of the sufficient reduction
end

"""
Like srproj but efficiently interates over a sparse counts matrix, and
only projects in a single direction (dir).
"""
function srproj(coefs::AbstractMatrix{T}, counts::SparseMatrixCSC, dir::Int; intercept=true, focusj=axes(counts,2)) where T
   ixoffset = intercept ? 1 : 0 # omitting the intercept
   n,d = size(counts)
   zm = zeros(n,2)
   φ = vec(coefs[ixoffset+dir,:]) # keep only desired directions
   rows = rowvals(counts)
   vals = nonzeros(counts)
   for j = 1:d
      for i in nzrange(counts, j)
         row = rows[i]
         val = vals[i]
         if in(j,focusj)
           zm[row,1] += val*φ[j]  # projection part
         end
         zm[row,2] += val       # m = total count
      end
   end
   for i=1:n
     mi = zm[i,2]
     if mi > 0
       # scale to get frequencies
       zm[i,1] /= mi
     end
   end
   zm # m is part of the sufficient reduction
end

"""
  Builds the design matrix X for predicting covar in direction projdir
  dmr version
  inz=[1] and testrank=false always for dmr, so variables are ignored and only here for convinence
    of unified calling function
"""
function srprojX(coefs::AbstractMatrix{T},counts,covars,projdir; includem=true, srprojargs...) where T
  # dims
  n,p = size(covars)
  # ixnotdir = 1:p .!= projdir
  ixnotdir = setdiff(1:p,[projdir])

  # design matrix w/o counts data
  X_nocounts = [ones(n) getindex(covars,:,ixnotdir)]

  # add srproj of counts data to X
  Z = srproj(coefs,counts,projdir; srprojargs...)
  X = [X_nocounts Z]

  if !includem
    # drop last column with total counts m
    X = X[:,1:end-1]
  end

  inz = [1]

  X, X_nocounts, inz
end
srprojX(m::DMRCoefs,counts,covars,projdir; inz=[1], testrank=false,
  kwargs...) = srprojX(coef(m),counts,covars,projdir; kwargs...)
srprojX(m::DMRPaths,counts,covars,projdir; select=defsegselect, inz=[1], testrank=false,
  kwargs...) = srprojX(coef(m, select),counts,covars,projdir; kwargs...)

"""
srproj for hurdle dmr takes two coefficent matrices
coefspos, coefszero, and a two specific directions
and returns an n-by-4 matrix Z = [zpos zzero m l].
dirpos = 0 omits positive counts projections and
dirzero = 0 omits zero counts projections.
Setting any of these to nothing will return projections in all directions.
"""
srproj(m::HDMRCoefs, counts, dirpos::D=0, dirzero::D=0; includel=includelinX(m), kwargs...) where {D<:Int}=
  srproj(coef(m)..., counts, dirpos, dirzero; intercept=hasintercept(m), includel=includel, kwargs...)
srproj(m::HDMRPaths, counts, dirpos::D=0, dirzero::D=0; includel=includelinX(m), select=defsegselect, kwargs...) where {D<:Int}=
  srproj(coef(m, select)..., counts, dirpos, dirzero; intercept=hasintercept(m), includel=includel, kwargs...)

"""
srproj for hurdle dmr takes two coefficent matrices
coefspos, coefszero, and a two specific directions
and returns an n-by-4 matrix Z = [zpos zzero m l].
dirpos = 0 omits positive counts projections and
dirzero = 0 omits zero counts projections.
Setting any of these to nothing will return projections in all directions.
"""
function srproj(coefspos::C, coefszero::C, counts, dirpos::D, dirzero::D;
  includem=true,        # whether to include total counts in Z
  includel=false,      # whether to include total 1(counts) in Z
  kwargs...) where {T, C<:AbstractMatrix{T}, D<:Int}

  if dirpos>0 && dirzero>0
    zpos = srproj(coefspos, counts, dirpos; kwargs...)
    zzero = srproj(coefszero, posindic(counts), dirzero; kwargs...)
    # second element should be same m in both, but because zero model
    # only sums indicators it generates smaller totals, so use the one
    # from the pos model
    # TODO: this needs to be fleshed out better in the theory to guide this choice
    # revised paper says to use l too
    Z = [zpos[:,1] zzero[:,1] zpos[:,2] zzero[:,2]]
    Z[:, [true, true, includem, includel]]
  elseif dirpos>0
    Z = srproj(coefspos, counts, dirpos; kwargs...)
    Z[:, [true, includem]]
  elseif dirzero>0
    Z = srproj(coefszero, posindic(counts), dirzero; kwargs...)
    Z[:, [true, includel]]
  else
    error("No direction to project to (dirpos=$dirpos,dirzero=$dirzero)")
  end
end

function ixcovars(dir::Int, inpos, inzero)
  # @assert dir ∈ inzero "projection direction $dir must be included in coefzero estimation!"
  # @assert dir ∈ inpos "projection direction $dir must be included in coefpos estimation!"

  # findfirst returns 0 if not found
  dirpos = something(findfirst(isequal(dir), inpos), 0)
  dirzero = something(findfirst(isequal(dir), inzero), 0)

  ineither = union(inzero,inpos)
  ixnotdir = setdiff(ineither,[dir])

  dirpos,dirzero,ineither,ixnotdir
end

"""
  Builds the design matrix X for predicting covar in direction projdir
  hdmr version
  Assumes that covars include all variables for both positives and zeros models
  and indicates which variables are where with the index arrays inpos and inzero.
  inz=[1,2] if both zpos and zzero are included
  inz=[2] if zpos is dropped due to collinearity
"""
function srprojX(coefspos::M, coefszero::M, counts, covars, projdir::Int;
  inpos=1:size(covars,2), inzero=1:size(covars,2),
  includem=true,        # whether to include total counts in Z
  includel=false,      # whether to include total 1(counts) in Z
  inz=[1,2], testrank=true, srprojargs...) where {T, M<:AbstractMatrix{T}}

  # dims
  n,p = size(covars)

  # get pos and zero subset indices
  dirpos,dirzero,ineither,ixnotdir = ixcovars(projdir, inpos, inzero)

  # design matrix w/o counts data
  X_nocounts = [ones(n) getindex(covars,:,ixnotdir)]

  includezpos = 1 ∈ inz
  # add srproj of counts data to X
  if includezpos
    Z = srproj(coefspos, coefszero, counts, dirpos, dirzero;
      includem=includem, includel=includel, srprojargs...)
    X = [X_nocounts Z]
  end

  if !includezpos || (testrank && rank(X) < size(X,2))
    if includezpos
      @info("rank(X) = $(rank(X)) < $(size(X,2)) = size(X,2). dropping zpos.")
    else
      @info("includezpos == false. dropping zpos.")
    end
    Z = srproj(coefszero, posindic(counts), dirzero; srprojargs...)
    X = [X_nocounts Z]
    inz = [2]
  end

  X, X_nocounts, inz
end

"Set the default value for whether to include obs-specific lexicon (l) in the srproj"
includelinX(m::HDMRCoefs{<:Hurdle}) = false
includelinX(m::HDMRPaths{<:Union{Missing, Hurdle}}) = false
includelinX(m::HDMRPaths{<:Hurdle}) = false

includelinX(m::HDMRCoefs{<:InclusionRepetition}) = true
includelinX(m::HDMRPaths{<:InclusionRepetition}) = true
includelinX(m::HDMRPaths{<:Union{Missing, InclusionRepetition}}) = true
includelinX(m::HDMRPaths{Missing}) = true

srprojX(m::HDMRCoefs,counts,covars,projdir; includel=includelinX(m), kwargs...) = srprojX(coef(m)...,counts,covars,projdir; inpos=m.inpos, inzero=m.inzero, includel=includel, kwargs...)
srprojX(m::HDMRPaths,counts,covars,projdir; includel=includelinX(m), select=defsegselect, kwargs...) = srprojX(coef(m, select)...,counts,covars,projdir; inpos=m.inpos, inzero=m.inzero, includel=includel, kwargs...)


text:
#Construct correlator from umbrella window time series.
function corr_i(W,mean)
    dt = 2 #2 fs timestep
    tCorr = 3000 #time over which correlations expected to die off.
    z =  W[:,2] .- mean # z time series (centered)
    nSamp = length(z)
    corr = zeros(tCorr+1,2) #initilize array that will contain correlation averges for a given t.
    norm = zeros(tCorr+1,1) #array of normalizations for different t.
    corr[1,1] = 0.0
    var0 = var(z)
    corr[1,2] = var0
    for i in 1:(nSamp) #loop over all samples (as separate time origins)
        tMax = nSamp
        if i+tCorr<nSamp
            tMax = i + tCorr
        end
        idx = 2 #shift index to account for  <δz(0)δ(0)> term.
        for j in (i+1):tMax
            t = (j - i)
            corr[idx,2] = corr[idx,2] + z[j]*z[i]
            corr[idx,1] = t*dt/1000 #conv to picoseconds.
            norm[idx] = norm[idx] + 1
            idx = idx + 1
        end
        #if i % tCorr == 0
        #    println(norm[2])
        #end
    end
    for k = 2:length(corr[:,1])
        corr[k,2] = corr[k,2]/norm[k]
    end

    #println(corr[length(corr[:,2]),1]," ",corr[length(corr[:,2]),2])

    ##################################
    #Below I make the resolution for longer times more coarse, to get rid of some of the noise
    #that might impact numerical integration.
    corr_cut = zeros(trunc(Int,1+750+2250/150),2)
    idx = 1
    for kk = 1:750
        corr_cut[kk,2] = corr[kk,2]
        corr_cut[kk,1] = corr[kk,1]
        idx = idx + 1
    end
    for jj in 751:150:(tCorr+1)
        corr_cut[idx,2] = corr[jj,2]
        corr_cut[idx,1] = corr[jj,1]
        idx = idx + 1
    end
    ###################################

    return corr, corr_cut, var0, tCorr
end

#Integrate correlation function from 0 → "∞". Trapezoid rule.
function Nint_corr(δδ)
    l = length(δδ[:,2])
    s = 0
    for i in 1:(l-1)
        s = s + (δδ[i+1,1] - δδ[i,1])*0.5*(δδ[i+1,2] + δδ[i,2])
    end
    return s
end


text:
module VND
using StaticArrays

import Base: +, *, /, ==, <, >, norm, dot
export VecN, norm, dist

type VecN
    elem::SVector;
end

sizeof(a::VecN) = length(a.elem)
sOK(a::VecN, b::VecN) = (sizeof(a) == sizeof(b)) ? true :
                         throw(BoundsError("Vector of different lengths"));

(+)(a::VecN, b::VecN) = [a.elem[i] + b.elem[i] for i in 1:sizeof(a) if sOK(a,b)]
(*)(x::Real, a::VecN) = [a.elem[i]*x for i in 1:sizeof(a)]
(*)(a::VecN, x::Real) = x*a
(/)(a::VecN, x::Real) = [a.elem[i]/x for i in 1:sizeof(a)]

(==)(a::VecN, b::VecN) = any([(a.elem[i] == b.elem[i]) for i in 1:sizeof(a) if sOK(a,b)])

dot(a::VecN, b::VecN) = sum([a.elem[i]*b.elem[i] for i in 1:sizeof(a) if sOK(a,b)])
norm(a::VecN) = sqrt(dot(a,a));
dist(a::VecN, b::VecN) = sum(map(x -> x*x,[a.elem[i]-b.elem[i] for i in 1:sizeof(a) if sOK(a,b)]))

end

#
# ----------
#  Examples
# ----------
#

using VND, StaticArrays

a = @SVector [1.1,2.2,3.3];
b = @SVector [2.1,3.2,4.3];

aa = VecN(a);
bb = VecN(b);

aa2 = 2.0*aa
bb3 = bb/3.0

cc = aa + bb
aa == bb

aanorm = norm(aa)
abdist = dist(aa,bb)

#
# Π can be calculated in higher dimensions
# For example the unit 4-ball has a volume of (Π * Π)/2.0
#
# see: https://en.wikipedia.org/wiki/Volume_of_an_n-ball
#


text:
# Autogenerated wrapper script for TensorRT_jll for x86_64-w64-mingw32-cuda+11.3
export libnvinfer, libnvinfer_plugin, libnvonnxparser, libnvparsers, trtexec

using CUDNN_jll
JLLWrappers.@generate_wrapper_header("TensorRT")
JLLWrappers.@declare_library_product(libnvinfer, "nvinfer.dll")
JLLWrappers.@declare_library_product(libnvinfer_plugin, "nvinfer_plugin.dll")
JLLWrappers.@declare_library_product(libnvonnxparser, "nvonnxparser.dll")
JLLWrappers.@declare_library_product(libnvparsers, "nvparsers.dll")
JLLWrappers.@declare_executable_product(trtexec)
function __init__()
    JLLWrappers.@generate_init_header(CUDNN_jll)
    JLLWrappers.@init_library_product(
        libnvinfer,
        "bin\\nvinfer.dll",
        nothing,
    )

    JLLWrappers.@init_library_product(
        libnvinfer_plugin,
        "bin\\nvinfer_plugin.dll",
        nothing,
    )

    JLLWrappers.@init_library_product(
        libnvonnxparser,
        "bin\\nvonnxparser.dll",
        nothing,
    )

    JLLWrappers.@init_library_product(
        libnvparsers,
        "bin\\nvparsers.dll",
        nothing,
    )

    JLLWrappers.@init_executable_product(
        trtexec,
        "bin\\trtexec.exe",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
module SmoothedParticles
import WriteVTK
import SparseArrays
import StaticArrays
import Match

include("structs.jl")
export AbstractParticle,
       ParticleSystem,
       ParticleField,
       DataField,
       RealVector, Vec2, Vec3,
       RealMatrix, Mat2, Mat3,
       VECX, VECY, VECZ, VEC0, MAT0, MAT1

include("kernels.jl")
export wendland2, Dwendland2, rDwendland2,
       wendland3, Dwendland3, rDwendland3, DDwendland3,
       spline24, Dspline24, rDspline24,
       spline23, Dspline23, rDspline23

include("core.jl")
export apply!,
       apply_unary!,
       apply_binary!,
       create_cell_list!,
       dist,
       assemble_matrix,
       assemble_vector,
       dot,
       norm

include("grids.jl")
export Grid,
       Squaregrid,
       Hexagrid,
       CubicGrid

include("geometry.jl")
export generate_particles!,
       Shape,
       Rectangle,
       Circle,
       Ellipse,
       Ball,
       Box,
       BooleanUnion,
       BooleanIntersection,
       BooleanDifference,
       Specification,
       BoundaryLayer,
       is_inside

include("IO.jl")
export save_frame!,
       new_pvd_file,
       save_pvd_file

end


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

using Markdown
using InteractiveUtils

# ╔═╡ a6648952-5a6d-11eb-2c81-99d1514983c0
using Econometrics,
RDatasets,
GLM,
GLMNet,
Statistics,
Plots,
DataFrames,
StatsPlots,
StatsBase

# ╔═╡ 8f1873f2-63b0-11eb-289f-4d880cee7e5a
md"Импортируем нужные библиотеки"

# ╔═╡ b842ca4a-63b1-11eb-2190-0b6c834f85e0
md"### Первый взгляд на датасет"

# ╔═╡ 2f38687c-63b3-11eb-384e-c1847038f4c5
md"Для построения регрессионых моделей воспользуемся датасетом **Auto** из `RDatasets`"

# ╔═╡ 093d42ea-5a72-11eb-19b1-6b0abb122099
auto = dataset("ISLR", "Auto")

# ╔═╡ ad79eb3c-63b5-11eb-3d99-1f2b7aab9a18
md"Выведем некоторую статистическую информацию по колонкам"

# ╔═╡ 22840156-5a72-11eb-0e57-ff3178e5ee6f
describe(auto)

# ╔═╡ c24ee52a-63d4-11eb-1815-0d184c506af2
md"**MPG** - это так называемые литры на 100 км. пути"

# ╔═╡ 61db537c-63d9-11eb-013c-9f74989da4d7
md"Сперва давайте посмотрим на графиках на зависимость между **MPG** и четыремя первыми переменными (после **Cylinders**)
"

# ╔═╡ fa02d9c0-63d9-11eb-1867-2f9780b312fe
scatter(auto.Displacement, auto.MPG, xaxis="Displacement", yaxis="MPG", legend = false)

# ╔═╡ 757e36ae-63d9-11eb-277a-9b2ac8f3b483
scatter(auto.Horsepower, auto.MPG, xaxis="Horsepower", yaxis="MPG", legend = false)

# ╔═╡ d46d4030-63d9-11eb-32d0-1780207fae50
scatter(auto.Weight, auto.MPG, xaxis="Weight", yaxis="MPG", legend = false)

# ╔═╡ e2e70718-63d9-11eb-02e5-0dff07a9c6c7
scatter(auto.Acceleration, auto.MPG, xaxis="Acceleration", yaxis="MPG", legend = false)

# ╔═╡ b26f1aac-63da-11eb-2e66-27cd5597e043
md"
Важно также посмотреть на *outliers* - так называемые выбросы. Подробности и тонкости работы с данными представлены в другом туториале, тут же мы посмотрим только на один важный момент :)

Можно построить **box plot** и после удалить выбросы, которые видны на графике
"

# ╔═╡ d987355c-63da-11eb-1170-910084cc4274
boxplot!(auto.Horsepower, title = "Box Plot - Horsepower", ylabel = "Horsepower", legend = false)

# ╔═╡ f0845d38-63db-11eb-0be2-eb5217bf6170
md"### Модели в `GLM` и `Econometrics`"

# ╔═╡ 17e60466-63c4-11eb-1f97-072f5539cb38
md"Для начала воспользуемся пакетом `GLM`. Начнем с простой парной регрессии: $$\widehat{MPG} = \widehat{\beta_0} + \widehat{\beta_1} \cdot Horsepower$$"

# ╔═╡ 8b294048-63bd-11eb-3df5-959b55069357
linear_model = lm(@formula(MPG ~ Horsepower), auto)

# ╔═╡ e86909f8-63c4-11eb-3f85-9dd2450e8c5f
md"Выведем коэффициент детерминации $$R^2$$"

# ╔═╡ 1fc60af8-63c1-11eb-0814-054e7b9be6f3
GLM.r2(linear_model)

# ╔═╡ f42c722a-63c4-11eb-262c-514bf90aa295
md"Horsepower объясняет примерно 60,6% дисперсии разброса MPG. Вполне неплохо, построим график с данными, регрессией и доверительным интервалом для предсказания с уровнем доверия 95%"

# ╔═╡ 4161b2ca-63c1-11eb-1f47-3f66d74863c3
begin
	pred = DataFrame(Horsepower = 45:0.1:230);
	pr = GLM.predict(linear_model, pred, interval = :prediction, level = 0.95);
	plot(xlabel="Horsepower", ylabel="MPG", legend=:best)
	plot!(auto.Horsepower, auto.MPG, label="data", seriestype=:scatter)
	plot!(pred.Horsepower, pr.prediction, label="model", linewidth=3,
        ribbon = (pr.prediction .- pr.lower, pr.upper .- pr.prediction))
end

# ╔═╡ 35163060-63c7-11eb-2013-d1878e8e2c0b
md"По графику видим, что связь возможно гиперболическая. Рассмотрим модель $$\widehat{MPG} = \widehat{\beta_0} + \widehat{\beta_1} \cdot \dfrac{1}{Horsepower}$$"

# ╔═╡ 65664ad0-63c6-11eb-3e2f-7bf9841d5265
hyperbolic_model = lm(@formula(MPG ~ 1 / Horsepower), auto)

# ╔═╡ 7582b77a-63c6-11eb-3ac2-615c4e294aab
GLM.r2(hyperbolic_model)

# ╔═╡ 5ba889e6-63c7-11eb-14a1-216dee432abe
md"Эта модель, похоже, и правда лучше объясняет **MPG**. Построим аналогичный график"

# ╔═╡ 8bb5ec74-63c6-11eb-1c3e-d358b4ebc5c5
begin
	predd = DataFrame(Horsepower = 45:0.1:232);
	prd = GLM.predict(hyperbolic_model, predd, interval = :prediction, level = 0.95);
	plot(xlabel="Horsepower", ylabel="MPG", legend=:best)
	plot!(auto.Horsepower, auto.MPG, label="data", seriestype=:scatter)
	plot!(predd.Horsepower, prd.prediction, label="model", linewidth=3,
        ribbon = (prd.prediction .- prd.lower, prd.upper .- prd.prediction))
end

# ╔═╡ 01019b9e-63dc-11eb-3148-9b08890b53c9
md"Теперь перейдем к множественной регрессии.

Мы видим, что есть линейная связь между **MPG** и **Horsepower** и **MPG** и **Weight**"

# ╔═╡ 0b27def8-63dc-11eb-003c-454b21b38132
model_GLM = lm(@formula(MPG ~ Weight + Horsepower), auto)

# ╔═╡ 38766bac-63dc-11eb-2dbc-978b5fb3a637
md"
Видим достаточно знакомую и родную табличку! Но фанаты **R** и **Stata** могут справедливо заметить, что как-то маловато всего представлено в таблице. Например, не видно коэффициента детерминации. Не унывать! **Julia** - за минимализм. Да и вывести всё, что Вас интересует, невероятно просто. Достаточно лишь указать интересующий метод и применить его к нашей модели.
"

# ╔═╡ 550be2bc-63dc-11eb-331b-d71a3071ec66
R2 = GLM.r2(model_GLM)

# ╔═╡ 6cd53b1e-63dc-11eb-074a-41c43b6dfaaf
coefs = GLM.coef(model_GLM)

# ╔═╡ 7df28c9e-63dc-11eb-3c8e-5732e36a7e5c
md"
Что ещё можно посмотреть, применив метод к модели?

*dof_residual*: степени свободы остатков (если имеет смысл)

*stderror*: стандартные ошибки коэффициентов

*vcov*: ковариационная матрица оценок коэффициентов

И много всего другого! Интересующиеся могут обратиться к документации `GLM`."

# ╔═╡ ce4f3232-63dc-11eb-374d-412a97f8d2c8
md"
А теперь давайте предсказать значение, задав заранее значения **Horsepower** и **Weight**.
"

# ╔═╡ de284b94-63dc-11eb-2d26-1391517ed17d
new_data = DataFrame(Horsepower = [200, 250], Weight = [4000, 3800])

# ╔═╡ 09221186-63dd-11eb-3c1a-a7fa2ef46e3e
prediction_GLM = DataFrame(MPG_pred = GLM.predict(model_GLM, new_data))

# ╔═╡ 3a4f1000-63e3-11eb-3f33-d1590d15b611
md"
Можем также посмотреть на имеющиеся наблюдения и сравнить с предсказанными значениями
"

# ╔═╡ 88873fa4-63e3-11eb-2d3a-e9d30a385d93
md"Посмотрим на распределение ошибок"

# ╔═╡ 76daf6ba-63e3-11eb-3e19-ff67bdb86a0e
begin
	performance = DataFrame(y_actual = auto.MPG, y_predicted =      GLM.predict(model_GLM));
	error = performance[:y_actual] - performance[:y_predicted];
	histogram(error, bins = 50, title = "Error Analysis", ylabel = "Frequency", 		xlabel = "Error",legend = false)
end

# ╔═╡ 92918540-63e3-11eb-1b57-3d49dfa75b09
md"Ошибки распределены нормально, судя по гистограмме"

# ╔═╡ c080265a-63e3-11eb-1aad-8f85c3d18f3e
md"
Также рассмотрим пакет `Econometrics`, который позволяет строить линейные регрессии, но немного отличается от `GLM`. Конечно, основные различия можно увидеть при более глубоком анализе, нежели тот, что мы приводим в нашем туториале, поэтому заинтересовавшиеся могут обратиться к документации пакета.
"

# ╔═╡ e66b840e-63e3-11eb-12f7-8d4155f58990
model_Econometrics = fit(EconometricModel, @formula(MPG ~ Weight + Horsepower), auto)

# ╔═╡ 6c1f9266-63e4-11eb-08bf-cbabc486dd30
md"
Фанаты **R** и **Stata** могут ликовать! Теперь мы видим табличку, в которой есть и коэффициент детерминации, и LR-тест, и много всего интересного для благородных Донов!
"

# ╔═╡ 0255de96-63e4-11eb-3336-fb9a753ae588
md"Сравним с аналогичной моделью из `GLM`"

# ╔═╡ 13df2904-63e4-11eb-0600-d1c1ebfe9d84
model_GLM

# ╔═╡ 2d6038be-63e4-11eb-2320-d7251a567785
md"Как мы видим, все статистические характеристики коэффициентов идентичны, только в `Econometrics` выводится еще разная полезная информация. В целом, пакет `Econometrics` умеет делать то же самое, что и пакет `GLM`"

# ╔═╡ da6bb162-665c-11eb-20f1-fb018a64646b
md"Может случиться такое, что мы построили модель, а коэффициент **Intercept** оказался незначимым, и хочется от него избавиться. Тогда в формуле надо просто записать ноль как слагаемое"

# ╔═╡ 055e78f2-665d-11eb-188e-8f88fb707c74
model_without_constant = fit(EconometricModel, @formula(MPG ~ 0 + Weight + Horsepower), auto)

# ╔═╡ 2afd1c5a-663b-11eb-26bf-5fbcc7b8b14a
md"### Работа с категориальными признаками"

# ╔═╡ 3e865ba6-663b-11eb-31e0-b57ea1e08bac
md"Если мы присмотримся к нашему датасету, то увидим, что у нас есть такой признак, как **Cylinders**, который мы можем рассмотреть, как категориальный. Убедимся в этом с помощью графика"

# ╔═╡ 6d6f598a-63da-11eb-0de2-fb70c755a0ef
bar!(auto.Cylinders, auto.MPG, xaxis="Cylinders", yaxis="Count", legend = false)

# ╔═╡ 0d648614-663c-11eb-0dc7-4de628f3f80d
md"Для того, чтобы рассматривать эту переменную как отдельные категории, создадим новый столбец с данными уже типа *categorical*"

# ╔═╡ deed8534-6635-11eb-3a64-37f4db496443
auto.Cat_cyl = categorical(auto.Cylinders)

# ╔═╡ 673274bc-663c-11eb-3766-03d3d4699611
md"Теперь построим регрессию с этим категориальным признаком"

# ╔═╡ e7fb9120-6635-11eb-2c9f-bdad924f344e
model_with_dummies = fit(EconometricModel, @formula(MPG ~ 1 + Cat_cyl), auto)

# ╔═╡ ab477a9e-663c-11eb-3bdb-7bbcc20dcaf8
md"Может смутить, что на графике четко видно, что **Cylinders** принимает 5 уникальных значений, а в выходе мы видим только четыре значения. Это связано с тем, что мы работает с дамми-переменными, и, чтобы избежать мультиколлинеарности, одну переменную не включают явно в уравнение регрессии (это верно только для модели с константой)"

# ╔═╡ ed569b2c-66be-11eb-3a87-bbc80441bde7
dummies_without_constant = fit(EconometricModel, @formula(MPG ~ 0 + Cat_cyl), auto)

# ╔═╡ 08585c6c-66bf-11eb-07e7-ab4d21eb9de2
md"Тут уже мы видим коэффициенты при всех категориях"

# ╔═╡ 8fca184a-6460-11eb-3685-9d8f282862fb
md"### Доверительные интервалы для коэффициентов"

# ╔═╡ 9b3e3364-6657-11eb-266d-29b8805c4a4c
md"Если вам нужен стандартный доверительный интервал (то есть доверительный интервал с уровнем доверия 95%), то достаточно просто посмотреть на выход модели, или вызвать следующую функцию"

# ╔═╡ 2df89a88-646f-11eb-353d-f9eec629f6e0
GLM.coeftable(model_GLM)

# ╔═╡ c17a0f14-6657-11eb-261a-837f9ca8bf81
md"Если нужен другой уровень доверия, то можно воспользоваться библиотекой `StatsBase`. Например, возьмем уровень доверия 80%"

# ╔═╡ 74743d1e-646f-11eb-3900-1789409e3f54
StatsBase.confint(model_GLM, level = 0.8)

# ╔═╡ e8f4608c-6657-11eb-2760-bbb016fc985e
StatsBase.confint(model_Econometrics, level = 0.8)

# ╔═╡ 20550cdc-6658-11eb-3359-b19b970c50de
md"Как мы видим, данная функция уживается и с `GLM`, и c `Econometrics`, и результаты не различаются, что не может не радовать"

# ╔═╡ 5b7c64c2-665d-11eb-0ea2-fb284873a9ed
md"Если нужен доверительный интервал с уровнем доверия 85% только для константы, тогда можно записать следующим образом"

# ╔═╡ 6ed3411c-665d-11eb-0e16-91bab8d97e52
StatsBase.confint(model_GLM, level = 0.85)[1, :]

# ╔═╡ 835a22cc-665d-11eb-12d1-19401e383788
md"Или для коэффициента при **Weight** при уровне доверия 92%"

# ╔═╡ 8156834e-665d-11eb-0c71-dba7bf9d54cb
StatsBase.confint(model_Econometrics, level = 0.92)[2, :]

# ╔═╡ d8a330c2-63db-11eb-2ad9-6b6a7126d393
md"### F-тест и LR-тест для сравнения моделей"

# ╔═╡ d3efe6bc-63b1-11eb-2ad2-f1df05dd812d
md"Рассмотрим две модели:

$$\widehat{MPG}_1 = \widehat{\beta_0} + \widehat{\beta_1} \cdot Weight + \widehat{\beta_2} \cdot Cylinders$$

$$\widehat{MPG}_2 = \widehat{\gamma_0} + \widehat{\gamma_1} \cdot Weight + \widehat{\gamma_2} \cdot Cylinders + \widehat{\gamma_3} \cdot Horsepower$$

F-тест позволит нам сравнить данные модели, и решить, является ли коэффициент $$\gamma_3$$ статистически значимым. Для большей ясности запишем гипотезы:

$$H_0: \gamma_3 = 0$$

$$H_1: \gamma_3 \neq 0$$

С помощью библиотеки `GLM` проведем F-тест
"

# ╔═╡ 1b8b1ad2-63b2-11eb-0d2e-75df725c1846
unrestricted_model = lm(@formula(MPG ~ Weight + Cylinders + Horsepower), auto)

# ╔═╡ 4617f266-63b2-11eb-257f-63906ea45f97
restricted_model = lm(@formula(MPG ~ Weight + Cylinders), auto)

# ╔═╡ 5a5e4770-63b2-11eb-3228-5dc5ea457b7b
GLM.ftest(unrestricted_model.model, restricted_model.model)

# ╔═╡ ee475c76-63b5-11eb-01ea-4121230e1b6f
md"Как мы видим, p-value, или p(>F), равен 0.0003, то есть он меньше любого разумного уровня значимости, отсюда $$\gamma_3 \neq 0$$ при любом разумном уровне значимости"

# ╔═╡ 52fca18e-665a-11eb-158f-69a0618e9b37
md"Так же мы можем провести LR-тест для сравнения наших моделей. Все те же условия, что и ранее"

# ╔═╡ 3637d41a-665a-11eb-27b7-5960cdc86bf4
GLM.lrtest(unrestricted_model.model, restricted_model.model)

# ╔═╡ 37b0229c-665b-11eb-3e81-4f52cf9babaf
md" Видим, что p-value, или p(>Chisq), очень маленький, то есть при любом разумном уровне значимости лучше выбирать неограниченную модель"

# ╔═╡ d65fc5c2-6215-11eb-3f82-2de3ce633b93
md"### Регуляризация"

# ╔═╡ b7f68a00-63b7-11eb-2789-1dce3762f0a6
md"Так же в **Julia** можно построить модели линейной регрессии с L1-регуляризацией (то есть **LASSO**) и L2-регуляризацией (то есть **Ridge**). Для этого воспользуемся библиотекой `GLMNet`. Заодно введем квадраты четырех переменных, поэлементно возводя их в квадрат (просто потому что)"

# ╔═╡ 17994b30-6216-11eb-05de-d1e8fef75c42
begin
	auto["Displacement_sq"] = auto["Displacement"].^2;
	auto["Horsepower_sq"] = auto["Horsepower"].^2;
	auto["Weight_sq"] = auto["Weight"].^2;
	auto["Acceleration_sq"] = auto["Acceleration"].^2;
end

# ╔═╡ 9ebe9d56-63b8-11eb-204f-fb4007c3392c
md"Важное уточнение: `GLMNet` работает не с датафреймами, а с массивами, поэтому выделим нужные признаки и передевем их в массивы"

# ╔═╡ b77d7afc-61e1-11eb-0e0f-55c1b299680c
X = hcat(Array(auto[:, 3:7]), Array(auto[:, 10:14]))

# ╔═╡ c2e747e2-61e1-11eb-3ac5-d7bc26e4b98e
y = Array(auto[:, 1])

# ╔═╡ e42b986c-63b8-11eb-32c7-45cbb84738bc
md"Сначала инициализируем модель **LASSO**, потом **Ridge**.

Подбираем оптимальные $$\lambda$$ с помощью кросс-валидации"

# ╔═╡ c6253aa6-61e4-11eb-214c-412220911078
lasso = glmnetcv(X, y, alpha = 1)

# ╔═╡ 17f4ee00-62d8-11eb-03f6-e78344c2183b
ridge = glmnetcv(X, y, alpha = 0)

# ╔═╡ fdd471ee-63b8-11eb-1231-f707db1c55d8
md"Хочется сравнить эти модели, напишем простую функцию для рассчета коэффциента детерминации. Тут нам нужна функция mean() из пакета `Statistics`"

# ╔═╡ d35f5116-61ee-11eb-1d2d-a70a56e4f5c4
function r2(model, X, y)
	RSS = sum((GLMNet.predict(model, X) - y).^2);
	TSS = sum((y .- Statistics.mean(y)).^2);
	(TSS - RSS)/TSS;
end

# ╔═╡ fdb7edfa-6212-11eb-07e4-19667961840b
r2(lasso, X, y)

# ╔═╡ 1fe3177c-62d8-11eb-0193-cb9420f3f4c7
r2(ridge, X, y)

# ╔═╡ 1f5c2190-63b9-11eb-302e-75f90c39c114
md"Тут модель **LASSO** показала себя лучше"

# ╔═╡ Cell order:
# ╟─8f1873f2-63b0-11eb-289f-4d880cee7e5a
# ╠═a6648952-5a6d-11eb-2c81-99d1514983c0
# ╟─b842ca4a-63b1-11eb-2190-0b6c834f85e0
# ╟─2f38687c-63b3-11eb-384e-c1847038f4c5
# ╠═093d42ea-5a72-11eb-19b1-6b0abb122099
# ╟─ad79eb3c-63b5-11eb-3d99-1f2b7aab9a18
# ╠═22840156-5a72-11eb-0e57-ff3178e5ee6f
# ╟─c24ee52a-63d4-11eb-1815-0d184c506af2
# ╟─61db537c-63d9-11eb-013c-9f74989da4d7
# ╠═fa02d9c0-63d9-11eb-1867-2f9780b312fe
# ╠═757e36ae-63d9-11eb-277a-9b2ac8f3b483
# ╠═d46d4030-63d9-11eb-32d0-1780207fae50
# ╠═e2e70718-63d9-11eb-02e5-0dff07a9c6c7
# ╟─b26f1aac-63da-11eb-2e66-27cd5597e043
# ╠═d987355c-63da-11eb-1170-910084cc4274
# ╟─f0845d38-63db-11eb-0be2-eb5217bf6170
# ╟─17e60466-63c4-11eb-1f97-072f5539cb38
# ╠═8b294048-63bd-11eb-3df5-959b55069357
# ╟─e86909f8-63c4-11eb-3f85-9dd2450e8c5f
# ╠═1fc60af8-63c1-11eb-0814-054e7b9be6f3
# ╟─f42c722a-63c4-11eb-262c-514bf90aa295
# ╠═4161b2ca-63c1-11eb-1f47-3f66d74863c3
# ╟─35163060-63c7-11eb-2013-d1878e8e2c0b
# ╠═65664ad0-63c6-11eb-3e2f-7bf9841d5265
# ╠═7582b77a-63c6-11eb-3ac2-615c4e294aab
# ╟─5ba889e6-63c7-11eb-14a1-216dee432abe
# ╠═8bb5ec74-63c6-11eb-1c3e-d358b4ebc5c5
# ╟─01019b9e-63dc-11eb-3148-9b08890b53c9
# ╠═0b27def8-63dc-11eb-003c-454b21b38132
# ╟─38766bac-63dc-11eb-2dbc-978b5fb3a637
# ╠═550be2bc-63dc-11eb-331b-d71a3071ec66
# ╠═6cd53b1e-63dc-11eb-074a-41c43b6dfaaf
# ╟─7df28c9e-63dc-11eb-3c8e-5732e36a7e5c
# ╟─ce4f3232-63dc-11eb-374d-412a97f8d2c8
# ╠═de284b94-63dc-11eb-2d26-1391517ed17d
# ╠═09221186-63dd-11eb-3c1a-a7fa2ef46e3e
# ╟─3a4f1000-63e3-11eb-3f33-d1590d15b611
# ╟─88873fa4-63e3-11eb-2d3a-e9d30a385d93
# ╠═76daf6ba-63e3-11eb-3e19-ff67bdb86a0e
# ╟─92918540-63e3-11eb-1b57-3d49dfa75b09
# ╟─c080265a-63e3-11eb-1aad-8f85c3d18f3e
# ╠═e66b840e-63e3-11eb-12f7-8d4155f58990
# ╟─6c1f9266-63e4-11eb-08bf-cbabc486dd30
# ╟─0255de96-63e4-11eb-3336-fb9a753ae588
# ╠═13df2904-63e4-11eb-0600-d1c1ebfe9d84
# ╟─2d6038be-63e4-11eb-2320-d7251a567785
# ╟─da6bb162-665c-11eb-20f1-fb018a64646b
# ╠═055e78f2-665d-11eb-188e-8f88fb707c74
# ╟─2afd1c5a-663b-11eb-26bf-5fbcc7b8b14a
# ╟─3e865ba6-663b-11eb-31e0-b57ea1e08bac
# ╠═6d6f598a-63da-11eb-0de2-fb70c755a0ef
# ╟─0d648614-663c-11eb-0dc7-4de628f3f80d
# ╠═deed8534-6635-11eb-3a64-37f4db496443
# ╟─673274bc-663c-11eb-3766-03d3d4699611
# ╠═e7fb9120-6635-11eb-2c9f-bdad924f344e
# ╟─ab477a9e-663c-11eb-3bdb-7bbcc20dcaf8
# ╠═ed569b2c-66be-11eb-3a87-bbc80441bde7
# ╟─08585c6c-66bf-11eb-07e7-ab4d21eb9de2
# ╟─8fca184a-6460-11eb-3685-9d8f282862fb
# ╟─9b3e3364-6657-11eb-266d-29b8805c4a4c
# ╠═2df89a88-646f-11eb-353d-f9eec629f6e0
# ╟─c17a0f14-6657-11eb-261a-837f9ca8bf81
# ╠═74743d1e-646f-11eb-3900-1789409e3f54
# ╠═e8f4608c-6657-11eb-2760-bbb016fc985e
# ╟─20550cdc-6658-11eb-3359-b19b970c50de
# ╟─5b7c64c2-665d-11eb-0ea2-fb284873a9ed
# ╠═6ed3411c-665d-11eb-0e16-91bab8d97e52
# ╟─835a22cc-665d-11eb-12d1-19401e383788
# ╠═8156834e-665d-11eb-0c71-dba7bf9d54cb
# ╟─d8a330c2-63db-11eb-2ad9-6b6a7126d393
# ╟─d3efe6bc-63b1-11eb-2ad2-f1df05dd812d
# ╠═1b8b1ad2-63b2-11eb-0d2e-75df725c1846
# ╠═4617f266-63b2-11eb-257f-63906ea45f97
# ╠═5a5e4770-63b2-11eb-3228-5dc5ea457b7b
# ╟─ee475c76-63b5-11eb-01ea-4121230e1b6f
# ╟─52fca18e-665a-11eb-158f-69a0618e9b37
# ╠═3637d41a-665a-11eb-27b7-5960cdc86bf4
# ╟─37b0229c-665b-11eb-3e81-4f52cf9babaf
# ╟─d65fc5c2-6215-11eb-3f82-2de3ce633b93
# ╟─b7f68a00-63b7-11eb-2789-1dce3762f0a6
# ╠═17994b30-6216-11eb-05de-d1e8fef75c42
# ╟─9ebe9d56-63b8-11eb-204f-fb4007c3392c
# ╠═b77d7afc-61e1-11eb-0e0f-55c1b299680c
# ╠═c2e747e2-61e1-11eb-3ac5-d7bc26e4b98e
# ╟─e42b986c-63b8-11eb-32c7-45cbb84738bc
# ╠═c6253aa6-61e4-11eb-214c-412220911078
# ╠═17f4ee00-62d8-11eb-03f6-e78344c2183b
# ╟─fdd471ee-63b8-11eb-1231-f707db1c55d8
# ╠═d35f5116-61ee-11eb-1d2d-a70a56e4f5c4
# ╠═fdb7edfa-6212-11eb-07e4-19667961840b
# ╠═1fe3177c-62d8-11eb-0193-cb9420f3f4c7
# ╟─1f5c2190-63b9-11eb-302e-75f90c39c114

max_stars_repo_path: Regression with Julia/regression.jl

max_issues_repo_path: Regression with Julia/regression.jl

max_forks_repo_path: Regression with Julia/regression.jl


text:
include(joinpath(pwd(), "src", "Heuristics.jl"))
include(joinpath(pwd(), "src", "DFS_BB.jl"))
include(joinpath(pwd(), "src", "VertexEliminationOrder.jl"))

import Test
import LightGraphs as LG
import GraphPlot as GP
import BenchmarkTools as BT
import ProfileView as PV

# Testing Heuristics
Test.@testset "VertexEliminationOrder.jl" begin
    # Set up
    graph_file = "circuit_graphs/qflex_line_graph_files_decomposed_true_hyper_true/test.gr"
    G = VertexEliminationOrder.graph_from_gr(graph_file)

    # min fill test
    g = copy(G)
    @BT.time mf = min_fill(g)
    Test.@test LG.nv(G) == length(mf[1])

    # min width test
    g = copy(G)
    @BT.time mw = min_width(g)
    Test.@test LG.nv(G) == length(mw[1])
end

# Testing DFS
Test.@testset "DFS" begin
    # Set up
    # graph_file = "circuit_graphs/qflex_line_graph_files_decomposed_true_hyper_true/rectangular_4x4_1-16-1_0.gr"
    # graph_file = "circuit_graphs/qflex_line_graph_files_decomposed_true_hyper_true/test.gr"
    # graph_file = "circuit_graphs/qflex_line_graph_files_decomposed_true_hyper_true/sycamore_53_20_0.gr" # very big
    graph_file = "circuit_graphs/qflex_line_graph_files_decomposed_true_hyper_true/sycamore_53_8_0.gr" # very big
    G = VertexEliminationOrder.graph_from_gr(graph_file)

    # min fill test
    seconds = 5.0
    println("Max alg duration: ", seconds, " seconds")
    @time r = branch_bound(G, seconds, 10)
    # @PV.profview r = branch_bound(G, seconds)

    println("Upper bound: ", r.ub)
    println("Best order: ", r.best_order)
    println("Timed out: ", r.timeout)
    println("Number of nodes visited: ", r.numNodes)
    println("Graph: ", r.g)

    # Checking to see if the function returns something
    Test.@test r !== nothing

    # Checking to see if the best order has unique values
    Test.@test length(unique(r.best_order)) == length(r.best_order)
end

max_stars_repo_name: OleeL/treewidth-heuristics

max_issues_repo_name: OleeL/treewidth-heuristics

max_forks_repo_name: OleeL/treewidth-heuristics


text:
module algT
    using GeometryTypes

    export Node, Edge, Obstacle, Room, AlgParameters, queueTmp, Graph, Pt2D, Line2D, Point

    Pt2D  = Point{2, Float64}
    UNDEF = -999999 #undefined
    Line2D = LineSegment{ Point{2, Float64}}

    mutable struct Node
        id::Int64
        state::Pt2D
        parentID::Int64
        cost::Float64

        Node(id::Int64, state::Pt2D, parentID::Int64 = -999, cost::Float64 = 9999999.) =
            new(id, state, parentID,  cost)
        end

    struct Edge
        startNode::algT.Node
        endNode::algT.Node
        #  startID::Int64
        #  endID::Int64
        #  startpt::Pt2D
        #  endpt::Pt2D
        #  line::Line2D
    end

    struct Obstacle
        id::Int64
        rect::HyperRectangle{2, Vec}
    end

    struct Room
        width::Int64
        height::Int64
        walls::Vector{LineSegment}
        obstacles::Vector{HyperRectangle}
    end

    struct AlgParameters
        numSamples::Int64
        connectRadius::Int64
    end

    struct queueTmp
        node::Node
        statesList::Vector{Node}
        cost::Int64
    end

    struct Graph
        startstate::Pt2D
        goalstate::Pt2D
        nodeslist::Vector{Node}
        edgeslist::Vector{Edge}
    end

    Base.isless(q1::queueTmp, q2::queueTmp) = q1.cost < q2.cost

    import Base.show
    function Base.show(io::IO, n::Node)
        print(io, "($(n.id), $(n.parentID), $(n.state))\n")
    end
    function Base.show(io::IO, e::Edge)
        print(io, "<$(e.startNode.id), $(e.startNode.state); $(e.endNode.id), $(e.endNode.state)>\n")
    end
end


text:
#######################################################################
# ModelBlock
#######################################################################

abstract type ModelTerm{F,S} end

struct Onsite{F,S} <: ModelTerm{F,S}
    o::F
    sublats::S              # SL === Missing means any sublats
end

struct Hopping{F,S,D,R} <: ModelTerm{F,S}
    h::F
    sublats::S              # S === Missing means any sublats
    ndists::D               # D  === Missing means any ndist
    range::R
end

isonsite(term::Onsite) = true
isonsite(term::Hopping) = false

(o::Onsite{F})(r, dr) where {F<:Function} = ensureSMatrix((o.o(r) + o.o(r)')/2)
(o::Onsite{S})(r, dr) where {S} = ensureSMatrix((o.o + o.o')/2)

(h::Hopping{F})(r, dr) where {F<:Function} = ensureSMatrix(h.h(r, dr))
(h::Hopping{S})(r, dr) where {S}  = ensureSMatrix(h.h)

Base.show(io::IO, o::Onsite) = print(io, "ModelTerm: onsite term on sublats = $(o.sublats)")
Base.show(io::IO, h::Hopping) =
    print(io, "ModelTerm: hopping term between sublat pairs = $(h.sublats), Bravais ndists = $(h.ndists) and maximum range = $(round(h.range, digits = 6))")

#######################################################################
# Model
#######################################################################
"""
    Model(terms::ModelTerm...)
    Model{Tv}(terms::ModelTerm...)

Define a `Model` of numeric type `Tv` (`Complex{Float64}` by default) from a collection of
`terms`, of type `Onsite` or `Hopping`. See the latter for more information.

# See also
    `Onsite`, `Hopping`
"""
struct Model{Tv,N,MB<:NTuple{N,ModelTerm}}
    terms::MB
end

Model{Tv}(hs::Vararg{ModelTerm,N}) where {Tv,N} = Model{Tv,N,typeof(hs)}(hs)
Model(hs...) = Model{Complex{Float64}}(hs...)

nonsites(m::Model) = count(isonsite, m.terms)

nhoppings(m::Model) = count(!isonsite, m.terms)

function Base.show(io::IO, m::Model{Tv,N}) where {Tv,N}
    print(io, "Model{$Tv,$N}: Model with $N terms ($(nonsites(m)) onsite and $(nhoppings(m)) hopping)")
end

"""
    Onsite(o; sublats = missing)

Define a `Onsite<:ModelTerm` that models an onsite energy `o` at sites on specific
`sublats`. Wrap several `ModelBlock`s in a `Model` to apply to a system.

 `o` can be a `o::SMatrix{N,N,Tv}` of dimension `NxN`, a number  `o::Tv`, or a function of
 position (`r -> o(r)`). The dimensions `N` of `o` will be checked for consistency between
 different terms when building a `System` with the model.

 `sublats` can be a single or a tuple of sublattice numbers or sublattice names. If
 `missing`, all sublattices are included.
"""
Onsite(o; sublats = missing) = _onsite(o, _normaliseSL(sublats))

_onsite(o, sublats) = Onsite(o, sublats)
_onsite(o::SMatrix{N,N}, sublats) where {N} = Onsite(o, sublats)
_onsite(o::SMatrix{N,M}, sublats) where {N,M} = throw(
    DimensionMismatch("Onsite energy must be a scalar or a square matrix"))
_onsite(o::AbstractArray, sublats) = Onsite(ensureSMatrix(o), sublats)

"""
    Hopping(h; sublats = missing, ndists = missing, range = 1)

Define a `Hopping<:ModelTerm` that models a hopping `t` between sites on specific `sublats`
and at a specific `range` and/or `ndists`. Wrap several `ModelTerms`s in a `Model` to apply
to a system.

`h` can be a `h::SMatrix{N,M,Tv}` of dimension `NxM`, a number `h::Tv`, or a function of two
positions (`(r,dr) -> h(r)`). Here, `r = (r1 + r2)/2` is the center of a link between two
sites (at `r1`, `r2`), while `dr = r2 - r1` is the link vector. The dimensions `N,M` of `h`
will be checked for consistency between different terms when building a `System` with the
model.

`sublats` can be a single or pair of sublattice numbers, sublattice
name or a tuple of pairs. If `missing`, all sublattices are included.

`ndists::Union{NTuple{D,NTuple{L,Int}},Missing}` specifies a tuple of
`D` unit cell distances (each a tuple of `L` integers) at which to
apply the hopping. If `missing` all `ndist`s will be included.

`range::Union{Real,Missing}` is a maximum Euclidean distance of sites to
link, or any distance if `missing`.
"""
Hopping(h; sublats = missing, ndists = missing, range = 1) =
    Hopping(ensureSMatrix(h), _normaliseSLpairs(sublats), _normaliseND(ndists),
    Float64(range) + sqrt(eps(Float64)))

_normaliseND(::Missing) = missing
_normaliseND(n::NTuple{L,Int}) where L = [SVector{L,Int}(n)]
_normaliseND(n::NTuple{N,NTuple{L,Int}}) where {N,L} = [SVector{L,Int}.(n)...]
_normaliseND(n::Vector{NTuple{L,Int}}) where {L} = SVector{L,Int}.(n)
_normaliseND(n) = throw(ErrorException(
    "`ndists` must be either `missing`, a `NTuple{L,Int}` or a list of them"))

_normaliseSLpairs(::Missing) = missing
_normaliseSLpairs(s::Union{Int,NameType}) = ((s,s),)
_normaliseSLpairs(ss::NTuple{N, Any}) where N =
    ntuple(n -> first(_normaliseSLpairs(ss[n])), Val(N))
_normaliseSLpairs(ss::Tuple{Union{Int,NameType}, Union{Int,NameType}}) = (ss,)
_normaliseSLpairs(n) = throw(ErrorException(
    "`sublats` must be either `missing` or a tuple of `s::Union{S, NTuple{2,S}}`, where `S<:Union{Int,$(NameType)}` specifies a sublattice number or name"))

_normaliseSL(::Missing) = missing
_normaliseSL(s::Union{Int,NameType}) = ((s,s),)
_normaliseSL(ss::NTuple{N, Any}) where N = ntuple(n -> first(_normaliseSL(ss[n])), Val(N))
_normaliseSL(n) = throw(ErrorException("`sublats` must be either `missing` or a tuple of `s::Union{S, NTuple{2,S}}`, where `S<:Union{Int,$(NameType)}` specifies a sublattice number or name"))


text:
module TestHP

using Distributions
using Test
using TreeParzen
import TreeParzen: Delayed

# Test that weights sum to 1

@test_throws ArgumentError HP.PChoice(:a, [
    Prob(0.1, 0),
    Prob(0.2, 1),
])

# Test that the space is constructed correctly

space = HP.PChoice(:a, [
    Prob(0.1, 0),
    Prob(0.2, 1),
    Prob(0.7, 2),
])
@test isa(Multinomial(1, space.choice.obj.probabilities), Multinomial)
@test isa(space, Delayed.AbstractSwitch)
@test length(space.options) == 3
@test isa(space.choice, Delayed.AbstractParam)
@test space.choice.label == :a
@test isa(space.choice.obj, Delayed.CategoricalIndex)
@test space.options == [0, 1, 2]

space = HP.PChoice(:b, [Prob(0.1, false), Prob(0.9, true)])
@test isa(space, Delayed.AbstractSwitch)
@test length(space.options) == 2
@test isa(space.choice, Delayed.AbstractParam)
@test space.choice.label == :b
@test isa(space.choice.obj, Delayed.CategoricalIndex)
@test space.options == [false, true]

@test_throws ArgumentError HP.PChoice(:c, Prob[])
@test_throws ArgumentError HP.Choice(:d, [])

end
true


text:
#!/usr/bin/env julia

# A Julia script that reproduces expected results of the test module
# 'test/ratmatTest.cpp'.
#
# From a shell, run the script as:
#   julia /path/to/ratmat.jl
#
# From Julia, the script can be run as:
#   include("path/to/ratmat.jl")

function ratmatTest()
    a = [ -1//2  2//5; 3//4  -1//3 ]
    println("a:")
    println(a);  println()

    invv = inv(a)
    println("inverse of a")
    println(invv);  println()

    a += invv
    println("a += invv:")
    println(a);  println()

    invv *= a
    println("determinant of inv * (a+inv): ", det(invv));  println()

    r = 3//4;
    b = a + r
    println("a + 3/4 :");
    println(b);  println()
    b = r + a
    println("3/4 + a :");
    println(b);  println()

    b = a - r
    println("a - 3/4 :");
    println(b);  println()
    b = r - a
    println("3/4 - a :");
    println(b);  println()

    r = 1//2;
    b = a;
    b += r
    println("a + 1/2 :");
    println(b);  println()
    b = a;
    b -= r
    println("a - 1/2 :");
    println(b);  println()
    b = a / r
    println("a / 1/2 :");
    println(b);  println()
    b /= r
    println("a / (1/2)^2 :");
    println(b);  println()

    b = a .* (a+r)
    println("a .* (a+1/2) :");
    println(b);  println()
    b = a;
    b .*= a
    println("a .* a :");
    println(b);  println()
    b = (a-r) ./ a
    println("(a-1/2) ./ a :");
    println(b);  println()
    b = a - r;
    b ./= (a + r)
    println("(a-1/2) ./ (a+1/2) :");
    println(b);  println()
end

ratmatTest()

max_stars_repo_name: jkovacic/math

max_issues_repo_name: jkovacic/math

max_forks_repo_name: jkovacic/math


text:
"""
$(SIGNATURES)

Find the only vertex in `graph` satifsying `f`. Function `f` must take arguments `graph`
and `vertex` and return a boolean.

Return the vertex.
"""
function find_vertex(f, graph)
    idx = findonly(v -> f(graph, v), collect(LightGraphs.vertices(graph)))
    return idx
end

"""
$(SIGNATURES)

Find the only edge in `graph` satifsying `f`. Function `f` must take arguments `graph`
and `edge` and return a boolean.

Return the edge.
"""
function find_edge(f, graph)
    edges = collect(LightGraphs.edges(graph))
    idx = findonly(v -> f(graph, v), edges)
    return edges[idx]
end

# Why do these not exist in LightGraphs.jl??
inedges(g, v) = (LightGraphs.edgetype(g)(u, v) for u in LightGraphs.inneighbors(g, v))
outedges(g, v) = (LightGraphs.edgetype(g)(v, u) for u in LightGraphs.outneighbors(g, v))

"""
Custom graph type that allows associating metadata with graph edges and vertices.

$(METHODLIST)
"""
struct MetaGraph{T,E,V} <: LightGraphs.AbstractGraph{T}
    graph::LightGraphs.SimpleDiGraph{T}
    edge_meta::E
    vertex_meta::V

    """
    $(SIGNATURES)

    Construct a `MetaGraph` from `graph` with edge metadata type `E` and vertes metadata
    type `V`.
    """
    function MetaGraph(graph::LightGraphs.SimpleDiGraph{T}, ::Type{E}, ::Type{V}) where {T,E,V}
        edge_meta = Dict{LightGraphs.edgetype(graph),E}()
        vertex_meta = Dict{eltype(graph),V}()
        return new{T,typeof(edge_meta),typeof(vertex_meta)}(graph, edge_meta, vertex_meta)
    end
end

const LIGHTGRAPHS_INTERFACE = (
    :(Base.reverse),
    :(LightGraphs.dst),
    :(LightGraphs.edges),
    :(LightGraphs.edgetype),
    :(LightGraphs.has_edge),
    :(LightGraphs.has_vertex),
    :(LightGraphs.inneighbors),
    :(LightGraphs.is_directed),
    :(LightGraphs.ne),
    :(LightGraphs.nv),
    :(LightGraphs.outneighbors),
    :(LightGraphs.src),
    :(LightGraphs.vertices),
)
for f in LIGHTGRAPHS_INTERFACE
    eval(:($f(M::MetaGraph, x...) = $f(M.graph, x...)))
end

LightGraphs.add_edge!(M::MetaGraph, src, dst) = LightGraphs.add_edge!(M.graph, src, dst)
function LightGraphs.add_edge!(M::MetaGraph, src, dst, metadata)
    success = LightGraphs.add_edge!(M.graph, src, dst)
    success || error()

    M.edge_meta[LightGraphs.edgetype(M.graph)(src, dst)] = metadata
    return true
end

LightGraphs.add_vertex!(M::MetaGraph) = LightGraphs.add_vertex!(M.graph)
function LightGraphs.add_vertex!(M::MetaGraph, metadata)
    success = LightGraphs.add_vertex!(M.graph)
    success || error()

    M.vertex_meta[LightGraphs.nv(M)] = metadata
    return true
end

getmeta(M::MetaGraph, v::Integer) = M.vertex_meta[v]
getmeta(M::MetaGraph, e) = M.edge_meta[e]

function rem_edge!(M::MetaGraph, e)
    success = LightGraphs.rem_edge!(M.graph, e)
    if success
        delete!(M.edge_meta, e)
    end
    return success
end

max_stars_repo_name: darchr/AutoTM

max_issues_repo_name: darchr/AutoTM

max_forks_repo_name: darchr/AutoTM


text:
#== # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
#
#   Tests related to orbit functions.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# References
#
#   [1] Vallado, D. A (2013). Fundamentals of Astrodynamics and Applications.
#       Microcosm Press, Hawthorn, CA, USA.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # ==#

# File: ./src/orbit/anomalies.jl
# ==============================

# Function M_to_E
# ---------------

################################################################################
#                                 Test Results
################################################################################
#
# Scenario 01
# ===========
#
# Example 2-1: Using Kepler's Equation [1, p. 66-67].
#
#   M = 235.4°, e = 0.4 ===> E = 220.512074767522º
#
# Using this function, the following result was obtained:
#
#   julia> M_to_E(0.4,235.4*pi/180)*180/pi
#   220.51207476752163
#
################################################################################

@testset "Function M_to_E" begin
    E = M_to_E(0.4, 235.4*pi/180)*180/pi
    @test E ≈ 220.512074767522 atol=1e-12
end

# File: ./src/transformation/gmst.jl
# ==================================

# Function JDtoGMST
# -----------------

################################################################################
#                                 Test Results
################################################################################
#
# Scenario 01
# ===========
#
# Example 3-5: Finding GMST and LST (Method 1) [1, p. 188].
#
# Considering the Julian Day [UT1] 2448855.009722, the Greenwich Mean Sideral
# Time was computed as 152.578787810°.
#
# Using SatelliteToolbox, the following was obtained:
#
#   julia> JDtoGMST(2448855.009722)*180/pi
#   152.57870762832462
#
# NOTE: The difference was also found by replicating the algorithm in MATLAB.
#
################################################################################

@testset "Function JDtoGMST" begin
    θ_GMST = JDtoGMST(2448855.009722)*180/pi
    @test θ_GMST ≈ 152.578787810 atol = 1e-4
end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule libpolymake_julia_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("libpolymake_julia")
JLLWrappers.@generate_main_file("libpolymake_julia", UUID("4d8266f6-2b3b-57e3-ad7a-d431eaaac945"))
end  # module libpolymake_julia_jll


text:
using Test
using PowerSpectra
using Healpix
using StableRNGs

@testset "binning matrix" begin
    b = binning_matrix([0,1,2], [0,1,3], ℓ->1; lmax=3)
    ref = [
        1.0  0.0  0.0  0.0
        0.0  1.0  0.0  0.0
        0.0  0.0  0.5  0.5
    ]
    @test all(ref .≈ b)
end

@testset "commented header reader" begin
    t = read_commented_header("data/commented_header.txt")
    @test all(t[!, "a"] .≈ [1.0, 2.0])
    @test all(t[!, "b"] .≈ [2.0, 3.0])
end

##

@testset "synalm" begin
    rng = StableRNG(123)
    nside = 16
    nsims = 100
    C0 = [3.  2.;  2.  5.]
    Cl = repeat(C0, 1, 1, 3nside);
    cls11 = Array{Float64, 2}(undef, (3nside, nsims))
    cls12 = Array{Float64, 2}(undef, (3nside, nsims))
    cls22 = Array{Float64, 2}(undef, (3nside, nsims))
    for i in 1:nsims
        alms = synalm(rng, Cl, nside)
        cls11[:,i] .= alm2cl(alms[1], alms[1])
        cls12[:,i] .= alm2cl(alms[1], alms[2])
        cls22[:,i] .= alm2cl(alms[2], alms[2])
    end
    @test isapprox(sum(cls11) ./ (nsims * (3nside)), 3.0, rtol=0.05)
    @test isapprox(sum(cls12) ./ (nsims * (3nside)), 2.0, rtol=0.05)
    @test isapprox(sum(cls22) ./ (nsims * (3nside)), 5.0, rtol=0.05)
end

##
@testset "fit_dipole" begin
    𝓂, 𝒹 = (0.5913756583917124,
        (0.0010145299203208928, 0.0009034021919689896, -0.07481236768343041))
    mask = readMapFromFITS("data/example_mask_1.fits", 1, Float64)
    𝓂′, 𝒹′ = fitdipole(mask)
    @test 𝓂 ≈ 𝓂′
    @test collect(𝒹) ≈ collect(𝒹′)
end

##
@testset "subtract monopole dipole" begin
    mask = readMapFromFITS("data/example_mask_1.fits", 1, Float64)
    𝓂, 𝒹 = fitdipole(mask)
    subtract_monopole_dipole!(mask, 𝓂, 𝒹)
    𝓂′, 𝒹′ = fitdipole(mask)
    @test 𝓂′ ≈ 0.0 atol=eps()
    @test collect(𝒹′) ≈ [0.0, 0.0, 0.0] atol=eps()
end

max_stars_repo_name: xzackli/PlanckCov

max_issues_repo_name: xzackli/PlanckCov


text:
using BDF, Base.Test

speedModes = [0,4,5,6,7,9]

for i in speedModes
    dats, evtTab, trigs, statusChan = readBDF(string("MK2_speedmode", i, "_CMS_not_in_range_battery_charged.bdf"))
    statusChanInfo = decodeStatusChannel(statusChan)
    @test isequal(unique(statusChanInfo["CMSInRange"]) , [false])
    @test isequal(unique(statusChanInfo["speedMode"]) , [i])
    @test isequal(unique(statusChanInfo["batteryLow"]) , [false])
    @test isequal(unique(statusChanInfo["isMK2"]) , [true])
end