text:
using FFTW
using StaticArrays
using Statistics

function msd_for_corr_time(traj_msd, corr_time, unwrapped = false)
    a = 0
    for k in 1:size(traj_msd)[3]
        for j in 1:size(traj_msd)[2]
            for i in 1:size(traj_msd)[1]-corr_time
                a += (traj_msd[i+corr_time, j, k] - traj_msd[i, j, k])^2
            end
        end
    end
    return a /= (size(traj_msd)[1] - corr_time) * size(traj_msd)[3]
end

function msd_direct(traj::Trajectory, targetatom, resolution, timerange)
    if !traj.unwrapped
        println("Trajectory is not unwrapped -> MSD might be faulty.")
    end

    return msd_direct(traj.coords, traj.timestep_in_fs, traj.atomlabels, targetatom, resolution, timerange)
end

function msd_direct(coords, timestep_fs, atom, targetatom, resolution, timerange)
    timesteps = size(coords)[3]
    corr_time_list = range(0, length = resolution, step = round(Int, timesteps * timerange / resolution))
    # Permuting the trajectory to shape (timesteps, dimensions, atoms) is faster, because it allows the innermost loop over the timesteps to be along the first dimension
    traj_msd = permutedims(coords[:, atom .== targetatom, :], [3, 1, 2])
    msd_result = Array{Float64, 1}(undef, resolution)
    Threads.@threads for i in 1:resolution
        msd_result[i] = msd_for_corr_time(traj_msd, corr_time_list[i])
    end
    # correlation time in ps, MSD in A^2
    return corr_time_list * timestep_fs / 1000, msd_result
end

function msd_fft(traj, timestep_fs, atom, targetatom, timerange)
    traj_msd = @view traj[:, atom .== targetatom, :]
    traj_msd = permutedims(traj_msd, [3, 1, 2])

    # compute S2
    fft_result = rfft(cat(dims = 1, traj_msd, zeros(Float64, size(traj_msd))), 1)
    fft_result = abs2.(fft_result) * size(traj_msd)[1]
    fft_result = irfft(fft_result, size(traj_msd)[1] * 2, 1)
end

function rdf(traj::Trajectory, atomtype1, atomtype2, d_min, d_max, bins)
    coord, atom, pbc = traj.coords, traj.atomlabels, traj.mdbox.pbc_matrix
    # Each bin i contains distances from distance_limits[i] to distance_limits[i+1]; distnace_mean[i] will be the average of both values
    distance_limits = range(d_min, d_max, length = bins + 1)
    distance_mean = (distance_limits[1:end-1] + distance_limits[2:end]) / 2
    # Shell volume of each bin i and overall ideal density used to normalize RDF
    shell_volumes = 4/3 * pi * (distance_limits[2:end].^3 - distance_limits[1:end-1].^3)
    cell_volume = det(pbc) #pbc[1,:] ⋅ (pbc[2,:] × pbc[3,:])
    ideal_density = count(atom .== atomtype1) * count(atom .== atomtype2) / cell_volume
    # Create distance histogram; delegate calculation of all distances and binning to create_dist_histogram_multi
    dist_histogram = create_dist_histogram_multi(traj, atomtype1, atomtype2, distance_limits)
    rdf = dist_histogram ./ (shell_volumes .* ideal_density .* size(coord)[3])
    return distance_mean, rdf
end

function create_dist_histogram_multi(traj::Trajectory, atomtype1, atomtype2, distance_limits)
    coord, atom, mdbox = traj.coords, traj.atomlabels, traj.mdbox
    dist_histogram = zeros(Int, distance_limits.len - 1)
    # Copy trajectories of specified atom types for later use; using @view here causes a lot more memory usage
    coord1 = coord[:, atom .== atomtype1, :]
    coord2 = coord[:, atom .== atomtype2, :]
    d_min, d_max, d_step = distance_limits[1], distance_limits[end], convert(Float64, distance_limits.step)
    # Split length of trajecgory $size(coord)[3] into roughly equal chunks, one for each available thread
    thread_limits = convert.(Int, floor.(range(1, size(coord)[3] + 1, length = Threads.nthreads() + 1)))
    # Lock to prevent two threads from editing dist_histogram
    Rlock = ReentrantLock()
    # Each thread will process one iteration of this loop; each will produce its own dist_histogram_thread and add it to the overall dist_histogram at the end
    Threads.@threads for i_thread in 1:Threads.nthreads()
        dist_histogram_thread = zeros(Int, distance_limits.len - 1)
        mdbox_thread = deepcopy(mdbox)
        i_limits = thread_limits[Threads.threadid()]:(thread_limits[Threads.threadid() + 1] - 1)
        # Iterate through chunk of timesteps set aside for this thread
        @views for i in i_limits
            # Iterate through atoms of both atom types
            for x in eachcol(coord1[:, :, i])
                for y in eachcol(coord2[:, :, i])
                    distance = pbc_dist(x, y, mdbox_thread)
                    bin!(dist_histogram_thread, distance, d_min, d_max, d_step)
                end
            end
        end
        # Add thread-local dist_histogram_thread to overall dist_histogram; lock it in the meantime in order to prevent threads from interfering with each other
        lock(Rlock)
        dist_histogram .+= dist_histogram_thread
        unlock(Rlock)
    end
    return dist_histogram
end

function bin!(dist_histogram, distance, d_min, d_max, d_step)
    # each bin i ∈ [1, bins] covers distances from (d_min + d_step * (i-1)) to (d_min + d_step * i)
    if d_min < distance < d_max
        index = floor(Int, (distance - d_min) / d_step + 1)
        dist_histogram[index] += 1
    end
end

function autocorrFFT(inputarray)
    arraylength = size(inputarray)[end]
    paddedlength = arraylength * 2 + 1
    paddedarray = zeros(eltype(inputarray), size(inputarray)[1:end-1]..., paddedlength)
    @views paddedarray[:, :, 1:arraylength] = inputarray

    FTresult = rfft(paddedarray, 3)
    FTresult .= FTresult .* conj(FTresult)
    autocorr = irfft(FTresult, paddedlength, 3)[:, :, 1:arraylength]
    for i in 1:arraylength
        @views autocorr[:, :, i] ./= (arraylength - i + 1)
    end
    return autocorr
end

function oacf(traj::Trajectory, atomtype1, atomtype2)
    atomno1 = count(traj.atomlabels .== atomtype1)
    atomno2 = count(traj.atomlabels .== atomtype2)
    frameno = size(traj.coords)[3]

    neighborsfirstframe = zeros(Int64, atomno1)
    frame1type2 = traj.coords[:, traj.atomlabels .== atomtype2, 1]
    @views for (i, pointtype1) in enumerate(eachcol(traj.coords[:, traj.atomlabels .== atomtype1, 1]))
        neighborsfirstframe[i] = TrajTools.next_neighbor(pointtype1, frame1type2, traj.mdbox)[1]
    end
    coordstype1 = @view traj.coords[:, traj.atomlabels .== atomtype1, :]
    coordstype2 = @view traj.coords[:, traj.atomlabels .== atomtype2, :]
    bondvectors = zeros(Float64, 3, atomno1, frameno)
    @views for frame in 1:size(traj.coords)[3]
        for atom1 in 1:count(traj.atomlabels .== atomtype1)
            neighborindex = neighborsfirstframe[atom1]
            coord1 = coordstype1[:, atom1, frame]
            coord2 = coordstype2[:, neighborindex, frame]
            bondvectors[:, atom1, frame] .= TrajTools.minimum_image_vector(coord1, coord2, traj.mdbox)
            normalize!(bondvectors[:, atom1, frame])
        end
    end
    auto1 = autocorrFFT(bondvectors)
    # orientational auto correlation
    oacf = vec(mean(sum(auto1, dims=1), dims = 2))
    # time in ps
    time_oacf = (1:frameno) * traj.timestep_in_fs / 1E3
    return time_oacf, oacf
end


text:
stderr_copy = Base.STDERR

rserr, wrerr = redirect_stderr()

using Weave
# For some reason this is needed on Windows, without it we see lots of errors
info("Ignore")

redirect_stderr(stderr_copy)
close(rserr)
close(wrerr)

input_file = readline()
input_file = chomp(input_file)

output_file = readline()
output_file = chomp(output_file)

doctype = readline()
doctype = chomp(doctype)

if doctype=="PREVIEW"
    template_path = joinpath(dirname(@__FILE__), "preview.tpl")

    Weave.weave(input_file, out_path=output_file, template=template_path)
else
    Weave.weave(input_file, out_path=:doc, doctype=doctype)
end

ext: jl

max_stars_repo_name: ChristianKurz/julia-vscode

max_issues_repo_name: ChristianKurz/julia-vscode

max_forks_repo_name: ChristianKurz/julia-vscode


text:
function list_notebooks(file, basedir = "")
    sp = splitpath(file)
    path = joinpath(@__DIR__, "..", "notebooks")
    filename = split(sp[end], "#")[1]
    list = join(["1. " * (nb == filename ?
                            name * " (本章)" :
                            "[$name](" * _linkname(path, nb, basedir) * ")")
                 for (name, nb) in NOTEBOOKS], "\n")
    Markdown.parse("""## 课程总目录 \n $list """)
end

function _linkname(path, nb, basedir)
    if haskey(ENV, "html_export") && ENV["html_export"] == "true"
        joinpath(basedir, "$(splitext(nb)[1]).html")
    else
        "open?path=" * joinpath(path, nb)
    end
end


text:
import Patchwork.div
import Base: size

export inset,
       empty,
       container,
       offset,
       width,
       height,
       minwidth,
       minheight,
       maxwidth,
       maxheight,
       size,
       topleft,
       midtop,
       topright,
       midleft,
       middle,
       midright,
       bottomleft,
       midbottom,
       bottomright,
       vertical,
       horizontal,
       depth,
       left,
       right,
       top,
       bottom,
       inward,
       outward,
       flow,
       hbox,
       vbox,
       vskip,
       hskip,
       grow,
       shrink,
       flex,
       wrap,
       axisstart,
       axisend,
       center,
       spacebetween,
       spacearound,
       stretch,
       packitems,
       packacross,
       packlines,
       space,
       pad,
       padcontent

# 0. Width and height

@api width => Width <: Tile begin
    typedarg(prefix::String="")
    arg(w::Length)
    curry(tile::Tile)
end

render(t::Width) =
    render(t.tile) & [:style => [(t.prefix == "" ? "width" : t.prefix * "Width")=> t.w]]

@api height => Height <: Tile begin
    typedarg(prefix::String="")
    arg(h::Length)
    curry(tile::Tile)
end

render(t::Height) =
    render(t.tile) & [:style => [(t.prefix == "" ? "height" : t.prefix * "Height")=> t.h]]

minwidth(w, x...) = width("min", w, x...)
minheight(h, x...) = height("min", h, x...)

maxwidth(w, x...) = width("max", w, x...)
maxheight(h, x...) = height("max", h, x...)

size(w::Length, h::Length, t) =
    width(w, height(h, t))
size(w::Length, h::Length) =
    t -> size(w, h, t)

container(w, h) =
    empty |> size(w, h)

# 1. Placing a Tile inside another
abstract Position
abstract Corner <: Position

@terms Corner begin
    topleft => TopLeft
    midtop => MidTop
    topright => TopRight
    midleft => MidLeft
    middle => Middle
    midright => MidRight
    bottomleft => BottomLeft
    midbottom => MidBottom
    bottomright => BottomRight
end

render{C <: Corner}(::C) = string(C)

immutable Relative{T <: Corner} <: Position
    x::Length
    y::Length
    # z::Length
end

@api inset => Inset <: Tile begin
    typedarg(position::Position=topleft)
    arg(containing::Tile)
    curry(contained::Tile)
end

render_position(p::TopLeft, x, y) =
    [:top => y, :left => x]
render_position(p::MidTop, x, y) =
    [:left =>  50cent, :top => y,
     :transform => "translate(-50%)",
     :marginLeft => x]
render_position(p::TopRight, x, y) =
    [:top => x, :right => y]
render_position(p::MidLeft, x, y) =
    [:top => 50cent, :left => x,
     :marginTop => y,
     :transform => "translate(0, -50%)"]
render_position(p::Middle, x, y) =
    [:top => 50cent, :left=>50cent,
     :marginLeft => x, :marginTop => y,
     :transform => "translate(-50%, -50%)"]
render_position(p::MidRight, x, y) =
    [:top => 50cent,
    :transform => "translate(0, -50%)",
    :marginTop => y, :right => x]
render_position(p::BottomLeft, x, y) =
    [:bottom => y, :left => x]
render_position(p::MidBottom, x, y) =
    [:left => 50cent, :bottom => y,
     :marginLeft => x,
     :transform => "translate(-50%)"]
render_position(p::BottomRight, x, y) =
    [:bottom => y, :right => x]

render_position(c::Corner) = [:style => render_position(c, 0, 0)]
render_position{C <: Corner}(p::Relative{C}) =
    [:style => render_position(C(), p.x, p.y)]

function render(tile::Inset)
    outer = render(tile.containing)
    inner = render(tile.contained)

    outer &= [:style => [:position => :relative]]
    inner &= [:style => [:position => :absolute]]

    outer << (inner & render_position(tile.position))
end

# 2. Axes, Directions and Flow

abstract Axis

abstract FixedAxis <: Axis
@terms FixedAxis begin
    horizontal => Horizontal
    vertical => Vertical
    depth => Depth
end

abstract  Side{Perpendicular <: Axis}
@terms Side{Horizontal} begin
    left => Left
    right => Right
end

@terms Side{Vertical} begin
    top => TopSide
    bottom => Bottom
end

name(s::Left) = "Left"
name(s::Right) = "Right"
name(s::TopSide) = "Top"
name(s::Bottom) = "Bottom"

@terms Side{Depth} begin
    inward => Inward
    outward => Outward
end

# TODO: render inward, outward flow

abstract FlowRelativeAxis <: Axis
@terms FlowRelativeAxis begin
    mainaxis => MainAxis
    crossaxis => CrossAxis
end

@terms Side{MainAxis} begin
    mainstart => MainStart
    mainend => MainEnd
end

@terms Side{CrossAxis} begin
    crossstart => CrossStart
    crossend => CrossEnd
end

abstract FlexContainer <: Tile

render(f::FlexContainer) =
    addclasses(render(f.tile), classes(f))

@api flow => Flow{A <: FixedAxis} <: FlexContainer begin
    typedarg(axis::A)
    curry(tiles::TileList)
    kwarg(reverse::Bool=false)
end
Flow{T}(x::T, y, z) = Flow{T}(x, y, z) # Julia issue 10641

flow(axis::FixedAxis, tiles...; reverse=false) =
    flow(axis, [t for t in tiles]; reverse=reverse)

Base.reverse(flow::Flow) =
    Flow(flow.axis, flow.tiles, !flow.reverse)

classes(f::Flow{Horizontal}) =
    f.reverse ? "flow flow-reverse horizontal" : "flow horizontal"
classes(f::Flow{Vertical}) =
    f.reverse ? "flow flow-reverse vertical" : "flow vertical"

render(f::Flow) =
    addclasses(render(f.tiles, :div), classes(f))

hbox(args...) = flow(horizontal, args...)
vbox(args...) = flow(vertical, args...)

hbox(arg) = flow(horizontal, [arg])
vbox(arg) = flow(vertical, [arg])

vskip(y) = size(0px, y, empty)
hskip(x) = size(x, 0px, empty)

@api ordering => FlowOrder <: FlexContainer begin
    arg(ordering::AbstractArray)
    curry(flow::FlexContainer)
end

# TODO: render ordering

@api wrap => Wrap <: FlexContainer begin
    arg(tile::FlexContainer)
    kwarg(reverse::Bool=false)
end
wrapreverse(t) = wrap(t, reverse=true)

classes(f::Wrap) =
    f.reverse ? "flex-wrap-reverse" : "flex-wrap"

# 3. Flexing and alignment

@api float => FloatingTile{T <: Side{Horizontal}} <: Tile begin
    typedarg(side::T)
    curry(tile::Tile)
end

render(f::FloatingTile) =
    render(tile) & [:style => [:float => lowercase(name(f.side))]]

@api grow => Grow <: Tile begin
    arg(factor::Float64)
    curry(tile::Tile)
end
grow(t::Tile) = grow(1.0, t)
grow(t::AbstractVector) = map(grow, t)

render(t::Grow) =
    render(t.tile) & [:style => [:flexGrow => t.factor]]

@api shrink => Shrink <: Tile begin
    arg(factor::Float64)
    curry(tile::Tile)
end
shrink(t::Tile) = shrink(1.0, t)
shrink(t::AbstractVector) = map(shrink, t)

render(t::Shrink) =
    render(t.tile) & [:style => [:flexShrink => t.factor]]

@api flexbasis => FlexBasis <: Tile begin
    arg(basis::Union(Length, Symbol))
    curry(tile::Tile)
end

render(t::FlexBasis) =
    render(t.tile) & [:style => [:flexBasis => t.basis]]

# Flex ignores the width and distributes forcefully
flex(factor::Real, t) =
    flexbasis(0mm, grow(factor, t))
flex(t) = flex(1.0, t)
flex() = flex(empty)
flex{T <: Real}(factor::AbstractVector{T}, t) =
    map(flex, factor, t)
flex(t::AbstractVector) =
    map(flex, t)

flex(factor::Real) = t -> flex(factor, t)
flex{T <: Real}(factor::AbstractVector{T}) = t -> flex(factor, t)

# Flow alignment
abstract Packing

@terms Packing begin
    axisstart => AxisStart
    axisend => AxisEnd
    center => AxisCenter
    stretch => Stretch
    baseline => Baseline
    spacebetween => SpaceBetween
    spacearound => SpaceAround
end

@api packlines => PackedLines{T <: Packing} <: FlexContainer begin
    typedarg(packing::T)
    curry(tile::FlexContainer)
end

classes(t::PackedLines{AxisStart}) = "pack-lines-start"
classes(t::PackedLines{AxisEnd}) = "pack-lines-end"
classes(t::PackedLines{AxisCenter}) = "pack-lines-center"
classes(t::PackedLines{Stretch}) = "pack-lines-stretch"
classes(t::PackedLines{SpaceBetween}) = "pack-lines-space-between"
classes(t::PackedLines{SpaceAround}) = "pack-lines-space-around"

@api packitems => PackedItems{T <: Packing} <: FlexContainer begin
    typedarg(packing::T)
    curry(tile::FlexContainer)
end

classes(t::PackedItems{AxisStart}) = "pack-start"
classes(t::PackedItems{AxisEnd}) = "pack-end"
classes(t::PackedItems{AxisCenter}) = "pack-center"
classes(t::PackedItems{SpaceBetween}) = "pack-space-between"
classes(t::PackedItems{SpaceAround}) = "pack-space-around"

@api packacross => PackedAcross{T <: Packing} <: FlexContainer begin
    typedarg(packing::T)
    curry(tile::FlexContainer)
end

classes(t::PackedAcross{AxisStart}) = "pack-across-start"
classes(t::PackedAcross{AxisEnd}) = "pack-across-end"
classes(t::PackedAcross{AxisCenter}) = "pack-across-center"
classes(t::PackedAcross{Stretch}) = "pack-across-stretch"
classes(t::PackedAcross{Baseline}) = "pack-across-baseline"

# 4. Padding

immutable Container <: Tile
    tile::Tile
end

render(cont::Container) = Elem(:div, render(cont.tile))

immutable Padded <: Tile
    sides::AbstractVector
    length::Length
    tile::Tile
end

render(padded::Padded) =
    render(padded.tile) &
        (isempty(padded.sides) ? # Apply padding to all sides if none specified
                [:style => ["padding" => padded.length]] :
                [:style => ["padding" * name(p) => padded.length for p=padded.sides]])

padcontent(len::Length, tile) =
    Padded(Side[], len, tile)

padcontent(sides::AbstractVector, len::Length, tile) =
    Padded(sides, len, tile)

pad(len::Length, tile) =
    padcontent(len, Container(tile))

pad(sides::AbstractVector, len::Length, tile) =
    padcontent(len, Container(tile))

pad(len::Length) =
    t -> pad(len, t)

pad(d::AbstractVector{Side}, len::Length) =
    t -> pad(d, len, t)


text:
type CrawlCorpus
    cloc::AbstractString
    debug::Bool
    function CrawlCorpus(cache_location::AbstractString, debug::Bool=false)
        isempty(cache_location) && error("cache location must be set to a valid directory")
        new(cache_location, debug)
    end
end

#TODO:
# server type
# character encoding
# server OS
# top level domain
# public suffix of domain
# second level domain
type ArchiveEntry
    uri::AbstractString
    mime::AbstractString
    len::Int
    http_status::AbstractString
    http_hdrs::Dict{ASCIIString,AbstractString}
    data::Array
end

function clear_cache(cc::CrawlCorpus)
    cc.debug && println("clearing cached files...")
    for f in readdir(cc.cloc)
        cc.debug && println("\t$f")
        rm(joinpath(cc.cloc, f))
    end
end

function clear_cache(cc::CrawlCorpus, archive::URI)
    fname = basename(archive.path)
    docsfile = joinpath(cc.cloc, fname)
    cc.debug && println("clearing cached file $fname")
    isfile(docsfile) && rm(docsfile)
    nothing
end

function segments(cc::CrawlCorpus)
    file = joinpath(cc.cloc, "valid_segments.txt")
    if !isfile(file)
        cc.debug && println("fetching valid segments...")
        t1 = time()
        os = open(file, "w")
        ho = HTTPClient.HTTPC.RequestOptions(ostream=os)
        get("http://aws-publicdatasets.s3.amazonaws.com/common-crawl/parse-output/valid_segments.txt", ho)
        close(os)
        cc.debug && println("\tfetched in $(time()-t1)secs")
    end
    segnames = AbstractString[]
    open(file) do f
        for str in readlines(f)
            push!(segnames, chomp(str))
        end
    end
    segnames
end

function archives(cc::CrawlCorpus, segment::AbstractString)
    file = joinpath(cc.cloc, string("segment_list_",segment,".txt"))
    arcnames = URI[]
    if !isfile(file)
        cc.debug && println("listing segment $segment")
        t1 = time()
        env = AWSEnv(timeout=60.0)
        segname = string("common-crawl/parse-output/segment/", segment)
        os = open(file, "w")
        opts = GetBucketOptions(prefix=segname)
        while true
            resp = S3.get_bkt(env, "aws-publicdatasets", options=opts)
            for elem in resp.obj.contents
                if endswith(elem.key, ".arc.gz")
                    uri_str = string("http://aws-publicdatasets.s3.amazonaws.com/", elem.key)
                    push!(arcnames, URI(uri_str))
                    println(os, uri_str)
                end
                opts.marker = elem.key
            end
            !resp.obj.isTruncated && break
        end
        close(os)
        cc.debug && println("\tfetched in $(time()-t1)secs")
    else
        cc.debug && println("opening cached file [$(file)]")
        open(file) do f
            for str in readlines(f)
                push!(arcnames, URI(chomp(str)))
            end
        end
    end
    cc.debug && println("$(length(arcnames)) archives in segment $(segment)")
    arcnames
end

function archives(cc::CrawlCorpus, count::Int=0)
    arcs = URI[]
    for seg in segments(cc)
        arcs_in_seg = archives(cc, seg)
        append!(arcs, arcs_in_seg)
        (count > 0) && (length(arcs) >= count) && break
    end
    (count == 0) ? arcs : arcs[1:count]
end

function open(cc::CrawlCorpus, archive::URI)
    fname = basename(archive.path)
    docsfile = joinpath(cc.cloc, fname)
    cc.debug && println("opening $archive. ($docsfile)")
    if !isfile(docsfile)
        cc.debug && println("\tdownloading $archive to $docsfile")
        t1 = time()
        os = open(docsfile, "w")
        ho = HTTPClient.HTTPC.RequestOptions(ostream=os)
        get(string(archive), ho)
        close(os)
        cc.debug && println("\tdownloaded in $(time()-t1)secs")
    end
    GZip.open(docsfile, "r", 8192*10)
end

function read_entry(cc::CrawlCorpus, f::IO, mime_part::AbstractString="", metadata_only::Bool=false)
    while true
        l = readline(f)
        while !eof(f) && isempty(l)
            l = readline(f)
        end
        eof(f) && isempty(l) && return ArchiveEntry("","",0,"",Dict{ASCIIString,AbstractString}(),[])
        vs = split(l)

        uri = vs[1]
        mime = vs[4]
        len = parse(Int, vs[5])

        if !isempty(mime_part) && !startswith(mime, mime_part)
            skip(f, len)
            continue
        end

        # read the http header
        hdrs = Dict{ASCIIString,AbstractString}()
        http_status = ""
        hdrlen = 0
        while !eof(f)
            l = readline(f)
            hdrlen += length(l)
            l = strip(l)
            isempty(l) && break
            nv = split(l, ':'; limit=2)
            if length(nv) == 2
                hdrs[lowercase(strip(nv[1]))] = strip(nv[2])
            elseif startswith(l, "HTTP")
                http_status = l
            end
        end

        len -= hdrlen
        # read or skip the http data
        data = []
        if metadata_only
            skip(f, len)
        else
            data = read(f, Array(Uint8, len))
        end
        return ArchiveEntry(uri, mime, len, http_status, hdrs, data)
    end
end

function read_entries(cc::CrawlCorpus, f::IO, mime_part::AbstractString="", num_entries::Int=0, metadata_only::Bool=false)
    arcs = ArchiveEntry[]
    while !eof(f)
        (num_entries > 0) && (length(arcs) >= num_entries) && break
        arc = read_entry(cc, f, mime_part, metadata_only)
        (!metadata_only) && isempty(arc.data) && continue
        push!(arcs, arc)
    end
    arcs
end


text:
# Smallest multiple

using Primes

LIMIT = 20

function solve()
    return prod(n^floor(log(n, LIMIT)) for n in primes(LIMIT))
end

print(Int(solve()))

max_stars_repo_name: wephy/project-euler

max_issues_repo_name: wephy/project-euler

max_forks_repo_name: wephy/project-euler


text:
using LinearAlgebra
using Random: MersenneTwister
using Distributions: MvNormal, logpdf
using StatsFuns: logsumexp

using SMC.Models: lg_smc_model, lgc_smc_model, LinearGaussian
using SMC: smc

F0(t) = [0.0, 0.0, 0.0]
G0(t) = [0.0, 0.0, 0.0]
F(t) = diagm([1.0, 1.0, 1.0])
G(t) = diagm(([1.0, 1.0, 1.0]))
Σ(t) = 9.0 * diagm(([1.0, 1.0, 1.0]))
Tau(t) = diagm([1.0, 1.0, 1.0])
μ0 = vcat(0.0, 0.0, 0.0)
Tau0 = 100.0 * diagm([1.0, 1.0, 1.0])

y = [
    5.22896384735402,
    2.059848221254017,
    11.77696919584194,
    5.721925811242492,
    2.647894413351219,
    8.833040488128724,
    10.112474414603064,
    6.208953089526956,
    7.483218024927158,
    1.1450746725213778,
    2.3827281539488214,
    2.5154588929528154,
    7.6938119335857795,
    11.267170019825972,
    11.203900809640977,
    4.641402249798912,
    6.393560809006704,
    5.434277785253692,
    7.754107613915341,
    2.8018519282190573,
    8.844915959029999,
    6.179788030963477,
    8.986422350276268,
    8.129873884597485,
    10.464689898594969,
    7.43846896751153,
    6.195006607054435,
    7.684383628629019,
    6.432819127902091,
    9.784040586681458,
    3.9386281953008275,
    1.0524045623884053,
    5.2287370014479535,
    3.3793567409786296,
    9.109449328047496,
    2.4107833351064665,
    3.1442786924940815,
    8.694980615700167,
    9.34957815303634,
    1.9899759308134595,
    1.7676330607869204,
    6.164419273710474,
    7.378003270378701,
    3.730353213883039,
    8.15997706282842,
    3.912665621137947,
    2.8222087161027987,
    7.89108559137954,
    2.016432756238303,
    -0.7283326731242079,
    5.34191426139306,
    3.7422889485477566,
    6.879539434261266,
    1.697715484841611,
    5.7492808736295835,
    3.4683312150848638,
    7.897465524526299,
    9.358577329317843,
    9.01915967340802,
    5.098410508423983,
    10.161486865442141,
    12.905026331072838,
    11.193444598591357,
    5.284455585283422,
    5.403785641006497,
    2.984985087125228,
    -0.1108114821742765,
    3.8443055205297396,
    1.1335781598255243,
    7.313228382179326,
    1.1291030218984508,
    0.8188795644543472,
    13.044240581728545,
    4.551403331575773,
    5.182374381380807,
    4.064021944452551,
    13.385904690258663,
    2.7383489638189253,
    4.272321073129483,
    1.936449829003382,
    7.579915913135302,
    6.891186573209264,
    5.696714256363617,
    2.9401339639559136,
    4.045792031292042,
    -2.69331630340043,
    4.253305831682493,
    -6.440609695250053,
    1.2224659946193035,
    3.6232261847393623,
    5.212788618742227,
    -0.5932923334934714,
    6.036279124314067,
    6.862343503023183,
    -2.9726591794689665,
    4.290653945351619,
    -2.420184182764869,
    0.6612926367123761,
    0.6312700637171165,
    -0.041609972691161345,
    -0.24839219509891972,
];

Y = hcat(y, y, y)
ll_true = -275.019326809115 * 3
K = 20
T = size(Y, 1)
rng = MersenneTwister(34)
for md in [lg_smc_model, lgc_smc_model]
    smc_fns = md(K, F0, F, G0, G, Tau, Σ, μ0, Tau0, Y)

    Xs = @inferred smc_fns.rinit(rng)
    dprs = @inferred smc_fns.dpr(Xs)
    dm1 = @inferred smc_fns.dm(Xs, 1)
    dinits = @inferred smc_fns.dinit(Xs)
    dpres = @inferred smc_fns.dpre(Xs, 2)

    @test size(Xs) == (3, K)

    dpr_actual = logpdf(MvNormal(μ0, Tau0), Xs)
    @test dprs ≈ dpr_actual

    dm_actual = logpdf(MvNormal(Y[1, :], Σ(1)), Xs)
    @test dm1 ≈ dm_actual

    ## The adapted filter should have all weights equal to the total-data likelihood
    logweights = dprs + dm1 - dinits + dpres
    @test all(logweights .≈ ll_true)

    Xs_new = @inferred smc_fns.rp(rng, Xs, 2)
    dts = @inferred smc_fns.dt(Xs, Xs_new, 2)
    dm2 = @inferred smc_fns.dm(Xs_new, 2)
    dps = @inferred smc_fns.dp(Xs, Xs_new, 2)

    #= using Debugger =#

    #= @enter smc_fns.rp(rng, Xs, 2) =#

    #= @enter smc_fns.dt(Xs, Xs_new, 2) =#

    @test size(Xs_new) == (3, K)

    dt_actual = logpdf(MvNormal(Tau(2)), Xs_new - G(2) * Xs)
    @test dts ≈ dt_actual

    dm_actual = logpdf(MvNormal(Y[2, :], Σ(2)), Xs_new)
    @test dm2 ≈ dm_actual

    @time smc_out = @inferred smc(rng, T, smc_fns...; threshold = 0.5)

    ## There will still be some slight noise in the particle filter, so this approx may fail
    @test smc_out.loglik ≈ ll_true

    # Test out new interface
    smc_filter = SMC.SMCModel(
        record_history = true,
        rinit = smc_fns.rinit,
        rproposal = smc_fns.rp,
        dinit = smc_fns.dinit,
        dprior = smc_fns.dpr,
        dproposal = smc_fns.dp,
        dtransition = smc_fns.dt,
        dmeasure = smc_fns.dm,
        dpre = smc_fns.dpre,
        threshold = 1.0,
    )
    smc_out = @inferred smc_filter(rng, T)
    smc_out = @time smc_filter(rng, T)
    ## There will still be some slight noise in the particle filter, so this approx may fail
    @test smc_out.loglik ≈ ll_true

    fwd_weights =
        smc_out.logweight_history[:, 27] +
        smc_fns.dpre(smc_out.particle_history[:, :, 27], 28)
    fwd_w = fwd_weights .- logsumexp(fwd_weights)
    @test all(isapprox.(fwd_w, -log(K)))

    @inferred SMC.calc_filtered_mean(smc_out)

    @inferred SMC.simulate_backward(
        rng,
        smc_out.particle_history,
        smc_out.logweight_history,
        smc_fns.dt,
        5,
    )
    @time SMC.simulate_backward(
        rng,
        smc_out.particle_history,
        smc_out.logweight_history,
        smc_fns.dt,
        5,
    )
end


text:
function boundaries_to_sparse(boundary)
    I, J, V = Int[], Int[], Bool[]
    for faceindex in boundary
        cell, face = faceindex
        push!(I, face)
        push!(J, cell)
        push!(V, true)
    end
    return sparse(I, J, V)
end

"""
    generate_grid(celltype::Cell, nel::NTuple, [left::Vec, right::Vec)

Return a `Grid` for a rectangle in 1, 2 or 3 dimensions. `celltype` defined the type of cells,
e.g. `Triangle` or `Hexahedron`. `nel` is a tuple of the number of elements in each direction.
`left` and `right` are optional endpoints of the domain. Defaults to -1 and 1 in all directions.
"""
generate_grid

# Line
function generate_grid(::Type{Line}, nel::NTuple{1,Int}, left::Vec{1,T}=Vec{1}((-1.0,)), right::Vec{1,T}=Vec{1}((1.0,))) where {T}
    nel_x = nel[1]
    n_nodes = nel_x + 1

    # Generate nodes
    coords_x = collect(linspace(left[1], right[1], n_nodes))
    nodes = Node{1,T}[]
    for i in 1:n_nodes
        push!(nodes, Node((coords_x[i],)))
    end

    # Generate cells

    cells = Line[]
    for i in 1:nel_x
        push!(cells, Line((i, i+1)))
    end

    # Cell faces
    boundary = Vector([(1, 1),
                       (nel_x, 2)])

    boundary_matrix = boundaries_to_sparse(boundary)

    # Cell face sets
    facesets = Dict("left"  => Set{Tuple{Int,Int}}([boundary[1]]),
                    "right" => Set{Tuple{Int,Int}}([boundary[2]]))
    return Grid(cells, nodes, facesets=facesets, boundary_matrix=boundary_matrix)
end

# QuadraticLine
function generate_grid(::Type{QuadraticLine}, nel::NTuple{1,Int}, left::Vec{1,T}=Vec{1}((-1.0,)), right::Vec{1,T}=Vec{1}((1.0,))) where {T}
    nel_x = nel[1]
    n_nodes = 2*nel_x + 1

    # Generate nodes
    coords_x = collect(linspace(left[1], right[1], n_nodes))
    nodes = Node{1,T}[]
    for i in 1:n_nodes
        push!(nodes, Node((coords_x[i],)))
    end

    # Generate cells
    cells = QuadraticLine[]
    for i in 1:nel_x
        push!(cells, QuadraticLine((2*i-1, 2*i+1, 2*i)))
    end

    # Cell faces
    boundary = Tuple{Int,Int}[(1, 1),
                         (nel_x, 2)]

    boundary_matrix = boundaries_to_sparse(boundary)

    # Cell face sets
    facesets = Dict("left"  => Set{Tuple{Int,Int}}([boundary[1]]),
                    "right" => Set{Tuple{Int,Int}}([boundary[2]]))
    return Grid(cells, nodes, facesets=facesets, boundary_matrix=boundary_matrix)
end

function _generate_2d_nodes!(nodes, nx, ny, LL, LR, UR, UL)
      for i in 0:ny-1
        ratio_bounds = i / (ny-1)

        x0 = LL[1] * (1 - ratio_bounds) + ratio_bounds * UL[1]
        x1 = LR[1] * (1 - ratio_bounds) + ratio_bounds * UR[1]

        y0 = LL[2] * (1 - ratio_bounds) + ratio_bounds * UL[2]
        y1 = LR[2] * (1 - ratio_bounds) + ratio_bounds * UR[2]

        for j in 0:nx-1
            ratio = j / (nx-1)
            x = x0 * (1 - ratio) + ratio * x1
            y = y0 * (1 - ratio) + ratio * y1
            push!(nodes, Node((x, y)))
        end
    end
end

function generate_grid(C::Type{Cell{2,M,N}}, nel::NTuple{2,Int}, X::Vector{Vec{2,T}}) where {M,N,T}
    @assert length(X) == 4
    generate_grid(C, nel, X[1], X[2], X[3], X[4])
end

function generate_grid(C::Type{Cell{2,M,N}}, nel::NTuple{2,Int}, left::Vec{2,T}=Vec{2}((-1.0,-1.0)), right::Vec{2,T}=Vec{2}((1.0,1.0))) where {M,N,T}
    LL = left
    UR = right
    LR = Vec{2}((UR[1], LL[2]))
    UL = Vec{2}((LL[1], UR[2]))
    generate_grid(C, nel, LL, LR, UR, UL)
end

# Quadrilateral
function generate_grid(C::Type{Quadrilateral}, nel::NTuple{2,Int}, LL::Vec{2,T}, LR::Vec{2,T}, UR::Vec{2,T}, UL::Vec{2,T}) where {T}
    nel_x = nel[1]; nel_y = nel[2]; nel_tot = nel_x*nel_y
    n_nodes_x = nel_x + 1; n_nodes_y = nel_y + 1
    n_nodes = n_nodes_x * n_nodes_y

    # Generate nodes
    nodes = Node{2,T}[]
    _generate_2d_nodes!(nodes, n_nodes_x, n_nodes_y, LL, LR, UR, UL)

    # Generate cells
    node_array = reshape(collect(1:n_nodes), (n_nodes_x, n_nodes_y))
    cells = Quadrilateral[]
    for j in 1:nel_y, i in 1:nel_x
        push!(cells, Quadrilateral((node_array[i,j], node_array[i+1,j], node_array[i+1,j+1], node_array[i,j+1])))
    end

    # Cell faces
    cell_array = reshape(collect(1:nel_tot),(nel_x, nel_y))
    boundary = Tuple{Int,Int}[[(cl, 1) for cl in cell_array[:,1]];
                              [(cl, 2) for cl in cell_array[end,:]];
                              [(cl, 3) for cl in cell_array[:,end]];
                              [(cl, 4) for cl in cell_array[1,:]]]

    boundary_matrix = boundaries_to_sparse(boundary)

    # Cell face sets
    offset = 0
    facesets = Dict{String, Set{Tuple{Int,Int}}}()
    facesets["bottom"] = Set{Tuple{Int,Int}}(boundary[(1:length(cell_array[:,1]))   + offset]); offset += length(cell_array[:,1])
    facesets["right"]  = Set{Tuple{Int,Int}}(boundary[(1:length(cell_array[end,:])) + offset]); offset += length(cell_array[end,:])
    facesets["top"]    = Set{Tuple{Int,Int}}(boundary[(1:length(cell_array[:,end])) + offset]); offset += length(cell_array[:,end])
    facesets["left"]   = Set{Tuple{Int,Int}}(boundary[(1:length(cell_array[1,:]))   + offset]); offset += length(cell_array[1,:])

    return Grid(cells, nodes, facesets=facesets, boundary_matrix=boundary_matrix)
end

# QuadraticQuadrilateral
function generate_grid(::Type{QuadraticQuadrilateral}, nel::NTuple{2,Int}, LL::Vec{2,T}, LR::Vec{2,T}, UR::Vec{2,T}, UL::Vec{2,T}) where {T}
    nel_x = nel[1]; nel_y = nel[2]; nel_tot = nel_x*nel_y
    n_nodes_x = 2*nel_x + 1; n_nodes_y = 2*nel_y + 1
    n_nodes = n_nodes_x * n_nodes_y

    # Generate nodes
    nodes = Node{2,T}[]
    _generate_2d_nodes!(nodes, n_nodes_x, n_nodes_y, LL, LR, UR, UL)

    # Generate cells
    node_array = reshape(collect(1:n_nodes), (n_nodes_x, n_nodes_y))
    cells = QuadraticQuadrilateral[]
    for j in 1:nel_y, i in 1:nel_x
        push!(cells, QuadraticQuadrilateral((node_array[2*i-1,2*j-1],node_array[2*i+1,2*j-1],node_array[2*i+1,2*j+1],node_array[2*i-1,2*j+1],
                                             node_array[2*i,2*j-1],node_array[2*i+1,2*j],node_array[2*i,2*j+1],node_array[2*i-1,2*j],
                                             node_array[2*i,2*j])))
    end

    # Cell faces
    cell_array = reshape(collect(1:nel_tot),(nel_x, nel_y))
    boundary = Tuple{Int,Int}[[(cl, 1) for cl in cell_array[:,1]];
                              [(cl, 2) for cl in cell_array[end,:]];
                              [(cl, 3) for cl in cell_array[:,end]];
                              [(cl, 4) for cl in cell_array[1,:]]]

    boundary_matrix = boundaries_to_sparse(boundary)

    # Cell face sets
    offset = 0
    facesets = Dict{String, Set{Tuple{Int,Int}}}()
    facesets["bottom"] = Set{Tuple{Int,Int}}(boundary[(1:length(cell_array[:,1]))   + offset]); offset += length(cell_array[:,1])
    facesets["right"]  = Set{Tuple{Int,Int}}(boundary[(1:length(cell_array[end,:])) + offset]); offset += length(cell_array[end,:])
    facesets["top"]    = Set{Tuple{Int,Int}}(boundary[(1:length(cell_array[:,end])) + offset]); offset += length(cell_array[:,end])
    facesets["left"]   = Set{Tuple{Int,Int}}(boundary[(1:length(cell_array[1,:]))   + offset]); offset += length(cell_array[1,:])

    return Grid(cells, nodes, facesets=facesets, boundary_matrix=boundary_matrix)
end

# Hexahedron
function generate_grid(::Type{Hexahedron}, nel::NTuple{3,Int}, left::Vec{3,T}=Vec{3}((-1.0,-1.0,-1.0)), right::Vec{3,T}=Vec{3}((1.0,1.0,1.0))) where {T}
    nel_x = nel[1]; nel_y = nel[2]; nel_z = nel[3]; nel_tot = nel_x*nel_y*nel_z
    n_nodes_x = nel_x + 1; n_nodes_y = nel_y + 1; n_nodes_z = nel_z + 1
    n_nodes = n_nodes_x * n_nodes_y * n_nodes_z

    # Generate nodes
    coords_x = linspace(left[1], right[1], n_nodes_x)
    coords_y = linspace(left[2], right[2], n_nodes_y)
    coords_z = linspace(left[3], right[3], n_nodes_z)
    nodes = Node{3,T}[]
    for k in 1:n_nodes_z, j in 1:n_nodes_y, i in 1:n_nodes_x
        push!(nodes, Node((coords_x[i], coords_y[j], coords_z[k])))
    end

    # Generate cells
    node_array = reshape(collect(1:n_nodes), (n_nodes_x, n_nodes_y, n_nodes_z))
    cells = Hexahedron[]
    for k in 1:nel_z, j in 1:nel_y, i in 1:nel_x
        push!(cells, Hexahedron((node_array[i,j,k], node_array[i+1,j,k], node_array[i+1,j+1,k], node_array[i,j+1,k],
                                 node_array[i,j,k+1], node_array[i+1,j,k+1], node_array[i+1,j+1,k+1], node_array[i,j+1,k+1])))
    end

    # Cell faces
    cell_array = reshape(collect(1:nel_tot),(nel_x, nel_y, nel_z))
    boundary = Tuple{Int,Int}[[(cl, 1) for cl in cell_array[:,:,1][:]];
                              [(cl, 2) for cl in cell_array[:,1,:][:]];
                              [(cl, 3) for cl in cell_array[end,:,:][:]];
                              [(cl, 4) for cl in cell_array[:,end,:][:]];
                              [(cl, 5) for cl in cell_array[1,:,:][:]];
                              [(cl, 6) for cl in cell_array[:,:,end][:]]]

    boundary_matrix = boundaries_to_sparse(boundary)

    # Cell face sets
    offset = 0
    facesets = Dict{String,Set{Tuple{Int,Int}}}()
    facesets["bottom"] = Set{Tuple{Int,Int}}(boundary[(1:length(cell_array[:,:,1][:]))   + offset]); offset += length(cell_array[:,:,1][:])
    facesets["front"]  = Set{Tuple{Int,Int}}(boundary[(1:length(cell_array[:,1,:][:]))   + offset]); offset += length(cell_array[:,1,:][:])
    facesets["right"]  = Set{Tuple{Int,Int}}(boundary[(1:length(cell_array[end,:,:][:])) + offset]); offset += length(cell_array[end,:,:][:])
    facesets["back"]   = Set{Tuple{Int,Int}}(boundary[(1:length(cell_array[:,end,:][:])) + offset]); offset += length(cell_array[:,end,:][:])
    facesets["left"]   = Set{Tuple{Int,Int}}(boundary[(1:length(cell_array[1,:,:][:]))   + offset]); offset += length(cell_array[1,:,:][:])
    facesets["top"]    = Set{Tuple{Int,Int}}(boundary[(1:length(cell_array[:,:,end][:])) + offset]); offset += length(cell_array[:,:,end][:])

    return Grid(cells, nodes, facesets=facesets, boundary_matrix=boundary_matrix)
end

# Triangle
function generate_grid(::Type{Triangle}, nel::NTuple{2,Int}, LL::Vec{2,T}, LR::Vec{2,T}, UR::Vec{2,T}, UL::Vec{2,T}) where {T}
    nel_x = nel[1]; nel_y = nel[2]; nel_tot = 2*nel_x*nel_y
    n_nodes_x = nel_x + 1; n_nodes_y = nel_y + 1
    n_nodes = n_nodes_x * n_nodes_y

    # Generate nodes
    nodes = Node{2,T}[]
    _generate_2d_nodes!(nodes, n_nodes_x, n_nodes_y, LL, LR, UR, UL)

    # Generate cells
    node_array = reshape(collect(1:n_nodes), (n_nodes_x, n_nodes_y))
    cells = Triangle[]
    for j in 1:nel_y, i in 1:nel_x
        push!(cells, Triangle((node_array[i,j], node_array[i+1,j], node_array[i,j+1]))) # ◺
        push!(cells, Triangle((node_array[i+1,j], node_array[i+1,j+1], node_array[i,j+1]))) # ◹
    end

    # Cell faces
    cell_array = reshape(collect(1:nel_tot),(2, nel_x, nel_y))
    boundary = Tuple{Int,Int}[[(cl, 1) for cl in cell_array[1,:,1]];
                               [(cl, 1) for cl in cell_array[2,end,:]];
                               [(cl, 2) for cl in cell_array[2,:,end]];
                               [(cl, 3) for cl in cell_array[1,1,:]]]

    boundary_matrix = boundaries_to_sparse(boundary)

    # Cell face sets
    offset = 0
    facesets = Dict{String,Set{Tuple{Int,Int}}}()
    facesets["bottom"] = Set{Tuple{Int,Int}}(boundary[(1:length(cell_array[1,:,1]))   + offset]); offset += length(cell_array[1,:,1])
    facesets["right"]  = Set{Tuple{Int,Int}}(boundary[(1:length(cell_array[2,end,:])) + offset]); offset += length(cell_array[2,end,:])
    facesets["top"]    = Set{Tuple{Int,Int}}(boundary[(1:length(cell_array[2,:,end])) + offset]); offset += length(cell_array[2,:,end])
    facesets["left"]   = Set{Tuple{Int,Int}}(boundary[(1:length(cell_array[1,1,:]))   + offset]); offset += length(cell_array[1,1,:])

    return Grid(cells, nodes, facesets=facesets, boundary_matrix=boundary_matrix)
end

# QuadraticTriangle
function generate_grid(::Type{QuadraticTriangle}, nel::NTuple{2,Int}, LL::Vec{2,T}, LR::Vec{2,T}, UR::Vec{2,T}, UL::Vec{2,T}) where {T}
    nel_x = nel[1]; nel_y = nel[2]; nel_tot = 2*nel_x*nel_y
    n_nodes_x = 2*nel_x + 1; n_nodes_y = 2*nel_y + 1
    n_nodes = n_nodes_x * n_nodes_y

    # Generate nodes
    nodes = Node{2,T}[]
    _generate_2d_nodes!(nodes, n_nodes_x, n_nodes_y, LL, LR, UR, UL)

    # Generate cells
    node_array = reshape(collect(1:n_nodes), (n_nodes_x, n_nodes_y))
    cells = QuadraticTriangle[]
    for j in 1:nel_y, i in 1:nel_x
        push!(cells, QuadraticTriangle((node_array[2*i-1,2*j-1], node_array[2*i+1,2*j-1], node_array[2*i-1,2*j+1],
                                        node_array[2*i,2*j-1], node_array[2*i,2*j], node_array[2*i-1,2*j]))) # ◺
        push!(cells, QuadraticTriangle((node_array[2*i+1,2*j-1], node_array[2*i+1,2*j+1], node_array[2*i-1,2*j+1],
                                        node_array[2*i+1,2*j], node_array[2*i,2*j+1], node_array[2*i,2*j]))) # ◹
    end

    # Cell faces
    cell_array = reshape(collect(1:nel_tot),(2, nel_x, nel_y))
    boundary = Tuple{Int,Int}[[(cl, 1) for cl in cell_array[1,:,1]];
                              [(cl, 1) for cl in cell_array[2,end,:]];
                              [(cl, 2) for cl in cell_array[2,:,end]];
                              [(cl, 3) for cl in cell_array[1,1,:]]]

    boundary_matrix = boundaries_to_sparse(boundary)

    # Cell face sets
    offset = 0
    facesets = Dict{String,Set{Tuple{Int,Int}}}()
    facesets["bottom"] = Set{Tuple{Int,Int}}(boundary[(1:length(cell_array[1,:,1]))   + offset]); offset += length(cell_array[1,:,1])
    facesets["right"]  = Set{Tuple{Int,Int}}(boundary[(1:length(cell_array[2,end,:])) + offset]); offset += length(cell_array[2,end,:])
    facesets["top"]    = Set{Tuple{Int,Int}}(boundary[(1:length(cell_array[2,:,end])) + offset]); offset += length(cell_array[2,:,end])
    facesets["left"]   = Set{Tuple{Int,Int}}(boundary[(1:length(cell_array[1,1,:]))   + offset]); offset += length(cell_array[1,1,:])

    return Grid(cells, nodes, facesets=facesets, boundary_matrix=boundary_matrix)
end

# Tetrahedron
function generate_grid(::Type{Tetrahedron}, nel::NTuple{3,Int}, left::Vec{3,T}=Vec{3}((-1.0,-1.0,-1.0)), right::Vec{3,T}=Vec{3}((1.0,1.0,1.0))) where {T}
    nel_x = nel[1]; nel_y = nel[2]; nel_z = nel[3]; nel_tot = 5*nel_x*nel_y*nel_z
    n_nodes_x = nel_x + 1; n_nodes_y = nel_y + 1; n_nodes_z = nel_z + 1
    n_nodes = n_nodes_x * n_nodes_y * n_nodes_z

    # Generate nodes
    coords_x = linspace(left[1], right[1], n_nodes_x)
    coords_y = linspace(left[2], right[2], n_nodes_y)
    coords_z = linspace(left[3], right[3], n_nodes_z)
    nodes = Node{3,T}[]
    for k in 1:n_nodes_z, j in 1:n_nodes_y, i in 1:n_nodes_x
        push!(nodes, Node((coords_x[i], coords_y[j], coords_z[k])))
    end

    # Generate cells, case 13 from: http://www.baumanneduard.ch/Splitting%20a%20cube%20in%20tetrahedras2.htm
    node_array = reshape(collect(1:n_nodes), (n_nodes_x, n_nodes_y, n_nodes_z))
    cells = Tetrahedron[]
    for k in 1:nel_z, j in 1:nel_y, i in 1:nel_x
        tmp = (node_array[i,j,k], node_array[i+1,j,k], node_array[i+1,j+1,k], node_array[i,j+1,k],
               node_array[i,j,k+1], node_array[i+1,j,k+1], node_array[i+1,j+1,k+1], node_array[i,j+1,k+1])
        push!(cells, Tetrahedron((tmp[1], tmp[2], tmp[4], tmp[5])))
        push!(cells, Tetrahedron((tmp[2], tmp[3], tmp[4], tmp[7])))
        push!(cells, Tetrahedron((tmp[2], tmp[4], tmp[5], tmp[7])))
        push!(cells, Tetrahedron((tmp[2], tmp[5], tmp[6], tmp[7])))
        push!(cells, Tetrahedron((tmp[4], tmp[5], tmp[7], tmp[8])))
    end
    # Cell faces
    cell_array = reshape(collect(1:nel_tot),(5, nel_x, nel_y, nel_z))
    boundary = Tuple{Int,Int}[[(cl, 1) for cl in cell_array[1,:,:,1][:]];
                              [(cl, 1) for cl in cell_array[2,:,:,1][:]];
                              [(cl, 2) for cl in cell_array[1,:,1,:][:]];
                              [(cl, 1) for cl in cell_array[4,:,1,:][:]];
                              [(cl, 2) for cl in cell_array[2,end,:,:][:]];
                              [(cl, 4) for cl in cell_array[4,end,:,:][:]];
                              [(cl, 3) for cl in cell_array[2,:,end,:][:]];
                              [(cl, 4) for cl in cell_array[5,:,end,:][:]];
                              [(cl, 4) for cl in cell_array[1,1,:,:][:]];
                              [(cl, 2) for cl in cell_array[5,1,:,:][:]];
                              [(cl, 3) for cl in cell_array[4,:,:,end][:]];
                              [(cl, 3) for cl in cell_array[5,:,:,end][:]]]

    boundary_matrix = boundaries_to_sparse(boundary)

    # Cell face sets
    offset = 0
    facesets = Dict{String, Set{Tuple{Int,Int}}}()
    facesets["bottom"] = Set{Tuple{Int,Int}}(boundary[(1:length([cell_array[1,:,:,1][:];   cell_array[2,:,:,1][:]]))   + offset]); offset += length([cell_array[1,:,:,1][:];   cell_array[2,:,:,1][:]])
    facesets["front"]  = Set{Tuple{Int,Int}}(boundary[(1:length([cell_array[1,:,1,:][:];   cell_array[4,:,1,:][:]]))   + offset]); offset += length([cell_array[1,:,1,:][:];   cell_array[4,:,1,:][:]])
    facesets["right"]  = Set{Tuple{Int,Int}}(boundary[(1:length([cell_array[2,end,:,:][:]; cell_array[4,end,:,:][:]])) + offset]); offset += length([cell_array[2,end,:,:][:]; cell_array[4,end,:,:][:]])
    facesets["back"]   = Set{Tuple{Int,Int}}(boundary[(1:length([cell_array[2,:,end,:][:]; cell_array[5,:,end,:][:]])) + offset]); offset += length([cell_array[2,:,end,:][:]; cell_array[5,:,end,:][:]])
    facesets["left"]   = Set{Tuple{Int,Int}}(boundary[(1:length([cell_array[1,1,:,:][:];   cell_array[5,1,:,:][:]]))   + offset]); offset += length([cell_array[1,1,:,:][:];   cell_array[5,1,:,:][:]])
    facesets["top"]    = Set{Tuple{Int,Int}}(boundary[(1:length([cell_array[4,:,:,end][:]; cell_array[5,:,:,end][:]])) + offset]); offset += length([cell_array[4,:,:,end][:]; cell_array[5,:,:,end][:]])

    return Grid(cells, nodes, facesets=facesets, boundary_matrix=boundary_matrix)
end


text:
function rosenbrock2d(x)
  return (1.0 - x[1])^2 + 100.0 * (x[2] - x[1]^2)^2
end

facts("bboptimize smoketest") do

 context("test each method option in short runs, just to ensure they work") do

   for(m in keys(BlackBoxOptim.ValidMethods))
     b, f = bboptimize(rosenbrock2d; method = m,
      search_space = [(-5.0, 5.0), (-2.0, 2.0)], max_time = 0.3,
      parameters = {:ShowTrace => false})

     # println("Fitness for $m: $f")

     @fact size(b) => (1, 2)
     @fact typeof(f) => Float64
     @fact f < 100.0 => true # this can't be very tight since we give very little time for optimization...
   end

 end

end


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    signarea(A, B, C)

Compute signed area of triangle formed
by points `A`, `B` and `C`.
"""
function signarea(A::Point{2}, B::Point{2}, C::Point{2})
  ((B - A) × (C - A)) / 2
end

"""
    sideof(point, segment)

Determines on which side of the oriented `segment`
the `point` lies. Possible results are `:LEFT`,
`:RIGHT` or `:ON` the segment.
"""
function sideof(p::Point{2,T}, s::Segment{2,T}) where {T}
  a, b = vertices(s)
  area = signarea(p, a, b)
  ifelse(area > atol(T), :LEFT, ifelse(area < -atol(T), :RIGHT, :ON))
end

"""
    sideof(point, chain)

Determines on which side of the closed `chain` the
`point` lies. Possible results are `:INSIDE` or
`:OUTSIDE` the chain.
"""
function sideof(p::Point{2,T}, c::Chain{2,T}) where {T}
  w = windingnumber(p, c)
  ifelse(isapprox(w, zero(T), atol=atol(T)), :OUTSIDE, :INSIDE)
end

"""
    mahalanobis(radii, angles, convention)

Return the Mahalanobis distance corresponding to an ellipsoid
with `radii` rotated by given `angles` according to given `convention`.

- For 2D ellipses, there are two radii and one rotation angle.
- For 3D ellipsoids, there are three radii and three rotation angles.

The list of available conventions can be found with:

```julia
julia> subtypes(RotationConvention)
```
"""
function mahalanobis(radii, angles, convention)
  ndims, nangles = length(radii), length(angles)
  valid = (ndims == 3 && nangles == 3) || (ndims == 2 && nangles == 1)
  @assert valid "invalid number of radii/angles"

  # invert radii if necessary
  invert = mainaxis(convention) == :Y
  ranges = invert ? [radii[i] for i in reverse(1:ndims,1,2)] : vec(radii)

  # scaling matrix
  λ = one(eltype(ranges)) ./ ranges.^2
  Λ = Diagonal(SVector{ndims}(λ))

  # rotation matrix
  R = rotmat(angles, convention)

  # ellipsoid metric
  W = Symmetric(R'*Λ*R)
  Mahalanobis(W)
end


text:
using LinearAlgebra

"""
    a, g, b, c = mapABCD(ABCD, form="CRFB")

Compute the coefficients for the specified structure.

See realizeNTF for a list of supported structures.

stuffABCD is the inverse function.
"""
function mapABCD(ABCD, form="CRFB")
    ABCD = copy(ABCD)

    order = size(ABCD, 1) - 1
    odd = rem(order, 2)
    even = 1- odd
    diagonal = 1:(order+2):order*(order+1)
    subdiag = diagonal[1:order] .+ 1
    supdiag = diagonal[2+odd:2:order] .- 1

    # Do the mapping.
    if form == "CRFB" || form == "CIFB" || form == "CRFBD"
        c = ABCD[subdiag]
        g = -ABCD[supdiag]
        if form == "CRFB"
            dly = (2+odd):2:order;  # row numbers of delaying integrators.
            ABCD[dly,:] = ABCD[dly,:] - diagm(c[dly.-1]) * ABCD[dly.-1,:]
        elseif form == "CRFBD"
            dly = (1+odd):2:order;  # row numbers of delaying integrators.
            ABCD[dly,:] = ABCD[dly,:] + diagm(g) * ABCD[dly.+1,:]
            if order > 2
                coupl = 2+even:2:order
                ABCD[coupl,:] = ABCD[coupl,:] - diagm(c[coupl.-1]) * ABCD[coupl.-1,:]
            end
        end
        a = -ABCD[1:order,order+2]
        b = ABCD[:,order+1]

    elseif form == "CRFF"
        c = [-ABCD[1,order+2]; ABCD[subdiag[1:end-1]]]
        g = -ABCD[supdiag]
        if even == 1
            multg = 1:2:order   # Rows to have g*(following row) added.
            ABCD[multg,:] = ABCD[multg,:] + diagm(g) * ABCD[multg.+1,:]
        end
        multc = 3:2:order # Rows to have c*(preceding row) subtracted.
        ABCD[multc,:] = ABCD[multc,:] - diagm(c[multc]) * ABCD[multc.-1,:]
        a = zeros(order)
        a[2:2:order] = ABCD[order+1,2:2:order]
        for i = 2:2:order     # Recover a coeff.
            ABCD[order+1,:] = ABCD[order+1,:] - a[i] * ABCD[i,:]
        end
        a[1:2:order] = ABCD[order+1,1:2:order]
        b = ABCD[:,order+1]

    elseif form == "CRFFD"
        # CRFFD has an extra order. Correct the common variables
        order = order - 1
        odd = rem(order, 2)
        even = 1 - odd
        diagonal = diagonal[1:order]
        subdiag = diagonal[1:order-1] .+ 1
        supdiag = diagonal[2+odd:2:order] .- 1

        g = -ABCD[supdiag]
        c = [-ABCD[1,order+3]; ABCD[subdiag]]
        a = zeros(order)
        for i = 1:2:order	# Recover the a coefficients
            a[i] = ABCD[order+1,i]
            ABCD[order+1,:] = ABCD[order+1,:] - a[i] * ABCD[i,:]
        end
        a[2:2:order] = ABCD[order+1,2:2:order]
        b = ABCD[1:order+1,order+2]
        for i = 2:2:order	# Recover the b coefficients
            b[i] = b[i] - c[i] * b[i-1]
            if odd == 1
                b[i] = b[i] + g[i÷2] * b[i+1]
            end
        end
        yscale = ABCD[order+2,order+1]
        a = a * yscale;
        b[end] = b[end] * yscale

    elseif form == "CIFF" || form == "Stratos"
        a = ABCD[order+1,1:order]
        c = [-ABCD[1,order+2]; ABCD[subdiag[1:end-1]]]
        g = -ABCD[supdiag]
        b = ABCD[:,order+1]

    elseif form == "DSFB" || form == "PFF"
        error("Form $form is not yet supported")
    else
        error("Unsupported form $form")
    end

    return a, g, b, c
end

max_stars_repo_name: Mahmoud-Kharsa/DeltaSigma

max_issues_repo_name: Mahmoud-Kharsa/DeltaSigma

max_forks_repo_name: Mahmoud-Kharsa/DeltaSigma


text:
using PyPlot
using LinearAlgebra
using SparseArrays
using Arpack
using ForwardDiff
using BenchmarkTools
using Formatting
using FileIO

#readfile
println("opening file... \n")
a = open("../programing/efelkismos-simple.unv","r")
#include("readfile.jl") #returns array of: nodes, edge, elem, bc
include("example-mesh.jl")
include("elastic_functions.jl") #defines basic functions

#start building global array
n = length(elem[:,1])
dofs = 2*length(nodes[:,1])
global E = Ematrix(220e9,0.5,"stress")
global K = sparse([],[],Float64[],dofs,dofs);
global M = sparse([],[],Float64[],dofs,dofs);

println("calculating stiffness/mass matrix...\n")
for i =1:n
    i1 = elem[i,:]'
    xs = nodes[i1,1]
    ys = nodes[i1,2]
    i1 = [i1 (i1 .+ Int64(0.5*dofs))]
    kre = ke(E,xs,ys)
    mre = 7e3 .* me(xs,ys)
    K[i1,i1] = K[i1,i1][1,:,1,:] + kre
    M[i1,i1] = M[i1,i1][1,:,1,:] + mre
end

include("rhs.jl") #calculate F

#gib traction nodes where (u,v) are known
#reminder, ix contains boundary condition nodes INDEXES
println("solving for U....\n")
freedofs = zeros(dofs,1)
fixedofs = ix[1][:]'#take first's b.c. nodes

fixedofs = [fixedofs[1] (fixedofs[1] .+Int64(0.5*dofs))]
freedofs[fixedofs] .= 1
freedofs = [i for i=1:dofs if freedofs[i]==0]

U = zeros(dofs,1)
rhs = F
rhs = rhs[freedofs] - K[freedofs,fixedofs'][:,:,1]*U[fixedofs']
Kfree = K[freedofs,freedofs]
Mfree = M[freedofs,freedofs]

Ufree = U[freedofs]
U[freedofs] = Kfree\rhs

n = length(elem[:,1])

println("calculating element average stress....\n")
#calculate element stress
include("stress.jl")
#println(stress')
println("DONE! \n")

i1 = elem[2,:]'
xs = nodes[i1,1]
ys = nodes[i1,2]

println(E*LN([0 0 1],xs,ys)*U[[ i1 (i1 .+ Int64(0.5*dofs))]]')
display(stress[:,1:3:18])
xs = [1 0 0 0.5 0 0.5]
ys = [0 -.5 .5 -.25 0 .25]
l = [1 0 0]

#=
F = zeros(12,1)
F[1] = 1e6
freedofs = ones(1,12)
fixedofs = [5 ]
freedofs[fixedofs] .= 0
freedofs = freedofs .==1
freedofs = freedofs[:]
kre = ke(E,xs,ys)

U  = zeros(12,1)
Ufree = kre[freedofs,freedofs]\F[freedofs]
U[freedofs] .= Ufree
println(E*LN([2 1 1]./4,xs,ys)*U)
figure(11)
for l1 = 0.:0.01:1.
    l2 = 0.5*(1-l1)
    l3 = l2
    sigma = E*LN([l1 l2 l3],xs,ys)*U
    scatter([l1],[sigma[1]/1e6])
    #scatter!([l1],[sigma[2]/1e6],legend=false,label="sy")
    #scatter!([l1],[sigma[3]/1e6],legend=false,label="txy")
end

xs = nodes[elem[2,:],1]
ys = nodes[elem[2,:],2]
xs = xs' ; ys  = ys';
g(x) = f(x)[1]
l = [1 1 0]
print(BC(g,xs,ys,l)')
=#

max_stars_repo_name: billyon/first

max_issues_repo_name: billyon/first

max_forks_repo_name: billyon/first


text:
using DataFrames
using VectorAutoregression
using Distributions
using Plots
pyplot()
default(show=true, reuse=false)

# ------------------------------------------------------------------------------------------
# Import Data
# ------------------------------------------------------------------------------------------
qdata = readtable("importable_data/qdata_processed.csv")
mdata = readtable("importable_data/mdata_processed.csv")
# Convert from DataFrame to Array
qdatamat = Matrix(qdata[[:GDP]])
mdatamat = Matrix(mdata[[:IP, :INFL, :UNEMP]])
# Dimensions
Kl = size(qdatamat, 2)
Kh = size(mdatamat, 2)
m = 3
lags = 1
constant = true
trend = false
nsave = 5000
nburn = 1000
h = 20
# Create the mixed-frequency matrix of endogenous variables
Ymix = mixed_freq_matrix(qdatamat, mdatamat, m)

# Aggregate monthly data to quarterly by taking final month values
aggregatedata = mdata[collect(range(3,3,length(qdata[:GDP]))) , :]
aggregatedata[:GDP] = qdata[:GDP]
Yagg = Matrix(aggregatedata[[:IP, :INFL, :UNEMP, :GDP]])

# ------------------------------------------------------------------------------------------
# Specify prior for coefficient matrix : Assuming all variables are AR(1) only
# ------------------------------------------------------------------------------------------

# AR coefficients for high-freq variables
rho_ip = 0.0
rho_infl = 0.9
rho_unemp = 0.9
rhos = [rho_ip, rho_infl, rho_unemp]
# AR coefficient for low-freq variable(s)
a = 0.0
# Initialize coefficient matrix
Aprior = zeros(Kh*m+Kl, Kh*m+Kl)
# Fill in high-freq block
Aprior[1:3, 7:9] = diagm(rhos)
Aprior[4:6, 7:9] = diagm(rhos.^2)
Aprior[7:9, 7:9] = diagm(rhos.^3)
# Fill in low-freq block
Aprior[10:10] = a
# Add constant term
if constant
    Aprior = hcat(zeros(Kh*m + Kl), Aprior)
end

# ------------------------------------------------------------------------------------------
# Specify prior for variance of coefficient matrix
# ------------------------------------------------------------------------------------------

# Hyperparameters
lambda = 1.0    # regulates tightness of prior: larger = higher-variance prior
vartheta = 1.0  # regulates effect of low-freq variables on high-freq variables: larger = higher-variance prior
# Initialize matrix
Vprior = zeros(Kh*m+Kl, Kh*m+Kl)
# prior variance on time (t-1) lags of the high-freq variables
for a = 1:m:Kh*m
    for b = 1:m:(Kh*(m-1))
        denom = (m - (div(b, m) + 1) + (div(a, m) + 1))^2
        Vprior[a:a+m-1, b:b+m-1] = ((lambda^2)/denom)*ones(Kh,Kh)
    end
end
# prior variance on AR coefficients for high-freq variables
for a = 1:m:Kh*m
    b = Kh*(m-1)+1
    denom = (div(a, m) + 1)^2
    Vprior[a:a+m-1, b:b+m-1] = ((lambda^2)/denom)*ones(Kh,Kh)
end
# prior for variance of effect of low-freq variables on high-freq variables
# Note: missing a scaling ratio S_hl present Ghysels (2016) (equiv., assuming S_hl = 1 for all)
for a = 1:m:Kh*m
    b = Kh*m + 1
    denom = (div(a, m) + 1)^2
    Vprior[a:a+m-1, b] = vartheta*((lambda^2)/denom)*ones(Kh)
end
# prior variance for coeffecients in low-freq block:
# 1) effect of high-freq variables on low-freq varables
for b = 1:m:Kh*m
    a = Kh*m + 1
    denom = (m - (div(b, m) + 1) + 1)^2
    Vprior[a:end, b:b+m-1] = ((lambda^2)/denom)*ones(Kl, Kh)
end
# 2) variance of low-freq AR coefficient
Vprior[Kh*m + 1, Kh*m + 1] = (lambda^2)/(m^2)
if constant
    const_var = (lambda^2)/(m^2)*ones(Kh*m + Kl)
    Vprior = hcat(const_var, Vprior)
end
# Transpose so they're in the proper orientation
Aprior = Aprior'
Vprior = Vprior'
# ------------------------------------------------------------------------------------------
# End prior specification
# ------------------------------------------------------------------------------------------

# Estimate MF-BVAR with default Independent Normal-Wishart prior
# varout = BVAR(Ymix, lags, constant, trend, nsave, nburn, [""])

# Estimte MF-BVAR with specified priors
var_mix = BVAR(Ymix, lags, constant, trend, nsave, nburn, [""]; Aprior=Aprior, Vprior=Vprior)

# Estimate BVAR with aggregate data
var_agg = BVAR(Yagg, lags, constant, trend, nsave, nburn, [""])

# IP Shock
# ------------------------------------------------------------------------------------------
# Calculate IRFs for MF-BVAR
irf_mix1 = VectorAutoregression.irf_bvar(var_mix, h, 1, 10)
irf_mix2 = VectorAutoregression.irf_bvar(var_mix, h, 4, 10)
irf_mix3 = VectorAutoregression.irf_bvar(var_mix, h, 7, 10)
# Calculate IRFs for BVAR with aggregate data
irf_agg1 = VectorAutoregression.irf_bvar(var_agg, h, 1, 4)

plot(irf_mix1[:point],  linewidth=2, label="Q1 shock", title="Response of GDP to IP shock")
plot!(irf_mix2[:point], linewidth=2, label="Q2 shock")
plot!(irf_mix3[:point], linewidth=2, label="Q3 shock")
plot!(irf_agg1[:point], linewidth=2, label="Aggregate VAR", color=:black)
plot!(irf_agg1[:lower], linewidth=2, linestyle=:dash, color=:black, label="")
plot!(irf_agg1[:upper], linewidth=2, linestyle=:dash, color=:black, label="")
savefig("IPshock.pdf")

# INFL Shock
# ------------------------------------------------------------------------------------------
# Calculate IRFs for MF-BVAR
irf_mix1 = VectorAutoregression.irf_bvar(var_mix, h, 2, 10)
irf_mix2 = VectorAutoregression.irf_bvar(var_mix, h, 5, 10)
irf_mix3 = VectorAutoregression.irf_bvar(var_mix, h, 8, 10)
# Calculate IRFs for BVAR with aggregate data
irf_agg1 = VectorAutoregression.irf_bvar(var_agg, h, 2, 4)

plot(irf_mix1[:point],  linewidth=2, label="Q1 shock", title="Response of GDP to INFL shock")
plot!(irf_mix2[:point], linewidth=2, label="Q2 shock")
plot!(irf_mix3[:point], linewidth=2, label="Q3 shock")
plot!(irf_agg1[:point], linewidth=2, label="Aggregate VAR", color=:black)
plot!(irf_agg1[:lower], linewidth=2, linestyle=:dash, color=:black, label="")
plot!(irf_agg1[:upper], linewidth=2, linestyle=:dash, color=:black, label="")
savefig("INFLshock.pdf")

# UNEMP Shock
# ------------------------------------------------------------------------------------------
# Calculate IRFs for MF-BVAR
irf_mix1 = VectorAutoregression.irf_bvar(var_mix, h, 3, 10)
irf_mix2 = VectorAutoregression.irf_bvar(var_mix, h, 6, 10)
irf_mix3 = VectorAutoregression.irf_bvar(var_mix, h, 9, 10)
# Calculate IRFs for BVAR with aggregate data
irf_agg1 = VectorAutoregression.irf_bvar(var_agg, h, 3, 4)

plot(irf_mix1[:point],  linewidth=2, label="Q1 shock", title="Response of GDP to UNEMP shock")
plot!(irf_mix2[:point], linewidth=2, label="Q2 shock")
plot!(irf_mix3[:point], linewidth=2, label="Q3 shock")
plot!(irf_agg1[:point], linewidth=2, label="Aggregate VAR", color=:black)
plot!(irf_agg1[:lower], linewidth=2, linestyle=:dash, color=:black, label="")
plot!(irf_agg1[:upper], linewidth=2, linestyle=:dash, color=:black, label="")
savefig("UNEMPshock.pdf")

max_stars_repo_name: lnsongxf/MF-VAR

max_issues_repo_name: lnsongxf/MF-VAR

max_forks_repo_name: lnsongxf/MF-VAR


text:
export neuron_dropping, neuron_linear_regression, population_linear_regression

function neuron_dropping(spikes::Array{SpikeTrain,1},analog::Array{Float64,2},binsize::Float64,nlags::UnitRange{Int64},firstwin::UnitRange{Int64},secondwin::UnitRange{Int64},orderinds::Array{Int64,1})

    neuroncorr=neuron_linear_regression(spikes,analog,binsize,nlags,firstwin)

    # Y is first x minutes of analog data
    Y = analog[firstwin, :]
    mywin1=collect(linspace((firstwin[1]-nlags[end]-1)*binsize, (firstwin[end]+1)*binsize, length(firstwin)+length(nlags)))

    # Y2 is second x minutes of analog data
    Y2 = analog[secondwin, :]
    mywin2=collect(linspace((secondwin[1]-nlags[end]-1)*binsize, (secondwin[end]+1)*binsize, length(secondwin)+length(nlags)))

    #find order of neurons from highest correlation to lowest correlation
    myinds=sortperm(squeeze(maximum(neuroncorr[:,orderinds],2),2), rev=true)

    #setup master neuron firing rate matrix, which is ordered X1, X2, X3 like above (so columns = neuron * (nlags+1), with order defined by highest correlation
    #first column and last column are ones

    neurbin1=Array(Float64,size(Y,1),(length(nlags)*length(spikes))+2)
    neurbin2=Array(Float64,size(Y,1),(length(nlags)*length(spikes))+2)

    neurbin1[:,1]=ones(Float64,size(neurbin1,1))
    neurbin2[:,1]=ones(Float64,size(neurbin2,1))

    count=2

    for i in myinds

        ts1=hist(spikes[i].ts,mywin1)[2]
        ts2=hist(spikes[i].ts,mywin2)[2]

        for j in nlags
            neurbin1[:,count]=ts1[(nlags[end]-j+1):(end-j)]
            neurbin2[:,count]=ts2[(nlags[end]-j+1):(end-j)]
            count+=1
        end

    end

    neurbin1[:,end]=ones(Float64,size(neurbin1,1))
    neurbin2[:,end]=ones(Float64,size(neurbin2,1))

    mincorr=zeros(Float64,length(spikes),size(analog,2))

    for i=length(spikes):-1:1
        myrange=((i-1)*(length(nlags))+2):size(neurbin1,2)
        X=view(neurbin1,:,myrange)
        coeffs= X \ Y
        X2=view(neurbin2,:,myrange)
        Z= X2 * coeffs
        mincorr[length(spikes)-i+1,:]=abs(diag(cor(Z,Y2)))
    end

    #max corr iteration

    maxcorr=zeros(Float64,length(spikes),size(analog,2))

    for i=1:length(spikes)

        myrange=1:((i*length(nlags))+1)
        X=view(neurbin1,:,myrange)
        coeffs= X \ Y
        X2=view(neurbin2,:,myrange)
        Z = X2 * coeffs
        maxcorr[i,:]=abs(diag(cor(Z,Y2)))

    end

    hcat(mincorr,maxcorr)
end

function neuron_linear_regression(spikes::Array{SpikeTrain,1},analog::Array{Float64,2},binsize::Float64,nlags::UnitRange{Int64},firstwin::UnitRange{Int64})

    neuroncorr=Array(Float64,length(spikes),size(analog,2))

    # Y is first x minutes of analog data
    Y = analog[firstwin, :]
    mywin1=collect(linspace((firstwin[1]-nlags[end]-1)*binsize, (firstwin[end]+1)*binsize, length(firstwin)+length(nlags)))

    X=zeros(Float64,size(Y,1),length(nlags)+1)
    X[:,1]=ones(Float64,size(Y,1),1)

    for i=1:length(spikes)

        #bin time stamps for first part
        ts=hist(spikes[i].ts,mywin1)[2]
        count=2
        for j in nlags
            X[:,count]=ts[(nlags[end]-j+1):(end-j)]
            count+=1
        end

        coeffs = X \ Y

        Z = X * coeffs

        neuroncorr[i,:]=abs(diag(cor(Z,Y)))

    end

    neuroncorr

end

function population_linear_regression(spikes::Array{SpikeTrain,1},analog::Array{Float64,2},binsize::Float64,nlags::UnitRange{Int64},firstwin::UnitRange{Int64})

    # Y is first x minutes of analog data
    Y = analog[firstwin, :]
    mywin1=collect(linspace((firstwin[1]-nlags[end]-1)*binsize, (firstwin[end]+1)*binsize, length(firstwin)+length(nlags)))

    X=zeros(Float64,size(Y,1),(length(nlags)*length(spikes))+1)
    X[:,1]=ones(Float64,size(Y,1),1)

    count=2

    for i=1:length(spikes)

        ts1=hist(spikes[i].ts,mywin1)[2]

        for j in nlags
            X[:,count]=ts1[(nlags[end]-j+1):(end-j)]
            count+=1
        end

    end

    coeffs = X \ Y

    Z = X * coeffs

    neuroncorr=abs(diag(cor(Z,Y)))

end

function population_linear_regression_predict(spikes::Array{SpikeTrain,1},analog::Array{Float64,2},binsize::Float64,nlags::UnitRange{Int64},firstwin::UnitRange{Int64},secondwin::UnitRange{Int64})

        # Y is first x minutes of analog data
    Y = analog[firstwin, :]
    mywin1=collect(linspace((firstwin[1]-nlags[end]-1)*binsize, (firstwin[end]+1)*binsize, length(firstwin)+length(nlags)))
    Y2 = analog[secondwin,:]

    X=zeros(Float64,size(Y,1),(length(nlags)*length(spikes))+1)
    X[:,1]=ones(Float64,size(Y,1),1)

    count=2

    for i=1:length(spikes)

        ts1=hist(spikes[i].ts,mywin1)[2]

        for j in nlags
            X[:,count]=ts1[(nlags[end]-j+1):(end-j)]
            count+=1
        end

    end

    coeffs = X \ Y

    X2 = zeros(Float64,size(Y2,1),(length(nlags)*length(spikes))+1)
    X2[:,1]=ones(Float64,size(Y2,1),1)

    mywin2=collect(linspace((secondwin[1]-nlags[end]-1)*binsize, (secondwin[end]+1)*binsize, length(secondwin)+length(nlags)))

    count=2

    for i=1:length(spikes)

        ts2=hist(spikes[i].ts,mywin2)[2]

        for j in nlags
            X2[:,count]=ts2[(nlags[end]-j+1):(end-j)]
            count+=1
        end

    end

    Z = X2 * coeffs

    neuroncorr=abs(diag(cor(Z,Y2)))

end


text:
using Econometrics, StatsBase
using BSON: @load
using BSON: @save

function MakeData()
    data = 0.0
    datadesign = 0.0
    S = TrainingTestingSize
    # training and testing
    for s = 1:S
        θ = rand(size(lb,1)).*(ub-lb) + lb
        m = WileE_model(θ)
        if s == 1
            data = zeros(S, size(vcat(θ, m),1))
        end
        data[s,:] = vcat(θ, m)
    end
    # save needed items with standard format
    @save "raw_data.bson" data
    return nothing
end


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

using Markdown
using InteractiveUtils

# ╔═╡ 3cd07b90-9e21-4a74-8d99-a6375beb7d1e
begin
	using Revise
	using Pkg
	Pkg.activate("../")

	using BayesianModels
	using DelimitedFiles
	using ForwardDiff
	using LinearAlgebra
	using Plots
	using Zygote

	include("./plotting.jl")
end

# ╔═╡ 6958bde2-c50f-4afd-8d08-7dfbfd97169a
md"""
# Learning and Modeling the Languages Phonetics

*[Lucas Ondel](https://lucasondel.github.io/), 8 March 2022*
"""

# ╔═╡ b83311ab-ca4f-432e-8627-5207517a81db
md"""## Lingua Libre Formants data

We use the [Lingua Libre Formants](https://github.com/lucasondel/LinguaLibreFormants) data for our analysis.

### Getting the data

We download the data and store it in a local file.
"""

# ╔═╡ d1588a1c-babe-4223-bc6d-f353314ed2c0
localfile = "formants.tsv"

# ╔═╡ 262e0395-814f-4315-b25e-21c18a7a9a5f
begin
	# Database URL
	local dataurl = "https://raw.githubusercontent.com/lucasondel/LinguaLibreFormants/main/data/formants.tsv"

	# Download the data if it is not done already
	if ! isfile(localfile)
		run(`wget $dataurl`)
	end
end

# ╔═╡ 4bbdd3d2-eb62-4889-9c8a-20141c0b140a
md"""
### Data preparation

In a nutshell:
1. we extract the list of languages, speakers and phones
2. organized the data as a $3 \times N$ matrix where 3 is for 3 first formants and $N$ is the number of data points
3. for each data point we extract the annotation (language, speaker, phone)
"""

# ╔═╡ ee1d1c45-70f5-4fe7-a8e2-9134d45a2b95
begin
	local rawdata = readdlm(localfile)

	speakers = Set()
	languages = Set()
	phones = Set()
	phones_per_lang = Dict()
	labels = []
	data = []
	for (lang, spk, sex, ipaphone, f1, f2, f3, duration) in eachrow(rawdata)
		# Remove the duration symbol
		ipaphone = replace(ipaphone, "ː" => "")

		lang, spk, ipaphone = Symbol(lang), Symbol(spk), Symbol(ipaphone)
		push!(data, vcat(f1, f2, f3))

		phone = ipaphone
		spk = (lang, sex, spk)

		l = get(phones_per_lang, lang, Set())
		push!(l, phone)
		phones_per_lang[lang] = l

		push!(labels, (lang, spk, phone))
		push!(speakers, spk)
		push!(languages, lang)
		push!(phones, phone)
	end
	data = hcat(data...)
	labels = vcat(labels...)

	# Mean-variance normalization of the data
	local x̄ = sum(data, dims=2) / size(data, 2)
	local x̄² = sum(data .^ 2, dims=2) / size(data, 2)
	local var = x̄² - (x̄ .^ 2)
	data = (data .- x̄) ./ sqrt.(var)

	# Transform the languages, speakers and phone sets to ordered collections
	languages = collect(languages)
	speakers = collect(speakers)
	phones = collect(phones)

	# Build the element <-> id mappings
	lang2id = Dict(lang => i for (i, lang) in enumerate(languages))
	id2lang = Dict(i => lang for (i, lang) in enumerate(languages))
	spk2id = Dict(spk => i for (i, spk) in enumerate(speakers))
	id2spk = Dict(i => spk for (i, spk) in enumerate(speakers))
	phone2id = Dict(phone => i for (i, phone) in enumerate(phones))
	id2phone = Dict(i => phone for (i, phone) in enumerate(phones))

	# Transform the labels to use the id rather than symbols
	labels = [(lang2id[lang], spk2id[spk], phone2id[phone])
			  for (lang, spk, phone) in labels]

	# Group the data for each triplet (lang, speaker, phone)
	groups = Dict()
	for (label, sample) in zip(labels, eachcol(data))
		list = get(groups, label, [])
		push!(list, sample)
		groups[label] = list
	end
	groups = tuple((k => hcat(groups[k]...) for k in keys(groups))...)
end

# ╔═╡ 474abeb1-82f8-4a74-a03c-be1cf55df840
md"Some important constants:"

# ╔═╡ 31e56fed-e4f0-4c71-abb2-5c2eff5a9a47
D = size(data, 1)

# ╔═╡ 61cc1ab2-a615-4a9b-91f9-bc5909e29297
N = size(data, 2)

# ╔═╡ e08cc384-49dc-4585-a2ab-fbec9fcfe9d2
L = length(languages)

# ╔═╡ b8e7f92f-0065-46fe-ac6b-5adf29b3379f
S = length(speakers)

# ╔═╡ ef71661f-338e-4d57-86ed-fc2c308a744f
P = length(phones)

# ╔═╡ c01de0cd-43f8-49cc-a8c9-ffbbc182d9e9
begin
	md"""
	### Data statistics

	| \# languages | \# speakers  | \# phones |
	|:------------:|:------------:|:---------:|
	| $L           | $S           | $P        |

	"""
end

# ╔═╡ 8b535f17-7a55-42da-ba7a-c161bec9fd28
md"""
## Model

### Embeddings
1.  $\boldsymbol{\lambda} \sim p(\boldsymbol{\lambda}) \triangleq \mathcal{N}( \mathbf{0}, \mathbf{I})$, $\boldsymbol{\lambda}$ → language embedding of dimension $Q_\lambda$
2.  $\boldsymbol{\sigma} \sim p(\boldsymbol{\sigma}) \triangleq \mathcal{N}( \mathbf{0}, \mathbf{I})$, $\boldsymbol{\sigma}$ → speaker embedding of dimension $Q_\sigma$
3.  $\boldsymbol{\pi} \sim p(\boldsymbol{\pi}) \triangleq \mathcal{N}( \mathbf{0}, \mathbf{I})$ $\boldsymbol{\pi}$ → phone embedding of dimension $Q_\pi$
"""

# ╔═╡ 5175712e-7f27-459e-8e6a-9bed325d6b4e
md"""
### Subspace construction
1.  $\mathscr{L}_i \sim p(\mathscr{L}) \triangleq \mathcal{N}(\mathbf{0}, \mathbf{I})$, $\{\mathscr{L}_0, \dots, \mathscr{L}_{Q_\lambda} \}$ → Language hyper-subspace bases of dimension $(Q_\pi + 1) \times K$
2.  $\mathscr{S}_j \sim p(\mathscr{S}) \triangleq \mathcal{N}(\mathbf{0}, \mathbf{I})$, $\{\mathscr{S}_0, \dots, \mathscr{S}_{Q_\sigma} \}$ → Speaker hyper-subspace bases $(Q_\pi + 1) \times K$
3.  $\mathscr{P} = \sum_{i=0}^{Q_\lambda} \lambda_i \mathscr{L}_i + \sum_{j=0}^{Q_\sigma} \sigma_j \mathscr{S}_j$ where $\lambda_0 = \sigma_0 = 1$, → Phonetic subspace bases of dimension $(Q_\pi + 1) \times K$
"""

# ╔═╡ 763b5f95-db02-4ede-b632-848a646232d9
md"""
### Phone model
1.  $\boldsymbol{\psi} = \mathscr{P}^\top \begin{bmatrix} \boldsymbol{\pi} \\ 1 \end{bmatrix}$, $\boldsymbol{\psi} \in \mathbb{R}^K$
2.  $\boldsymbol{\eta} = f(\boldsymbol{\psi})$
3.  $\mathbf{x}_n \sim p(\mathbf{x} | \boldsymbol{\eta}) \triangleq \mathcal{N}(\boldsymbol{\mu}, \boldsymbol{\Lambda}^{-1})$, where $\boldsymbol{\eta}$ are the natural parameters of the multivariate Normal distribution
"""

# ╔═╡ 841a3f0e-8975-49d2-959e-afe832ed41e4
md"""
### Parameter function mapping

The function $f$ transforms the parameters $\boldsymbol{\psi}$ from the real space into the natural parameter space. It is defined as follows:

```math
	\begin{align}
	\boldsymbol{\psi} &= \begin{bmatrix}
		 \boldsymbol{\psi}^{(\boldsymbol{\Lambda\mu})} \\
		 \boldsymbol{\psi}^{(\mathbf{L})}_1 \\
		 \boldsymbol{\psi}^{(\mathbf{L})}_2
	\end{bmatrix} \\
	\boldsymbol{\Lambda\mu} &= \boldsymbol{\psi}^{(\boldsymbol{\Lambda\mu})} \\
	\text{diag}(\mathbf{L}) &= \exp \{ \boldsymbol{\psi}^{(\mathbf{L})}_1 \} \\
	\text{vech}^{\dagger}(\mathbf{L}) &= \boldsymbol{\psi}^{(\mathbf{L})}_2 \\
	\boldsymbol{\Lambda} &= \mathbf{L} \mathbf{L}^\top
	\end{align}
```
where $\text{vech}^{\dagger}$ is the half-vectorization of a matrix **not including the diagonal**
"""

# ╔═╡ b0761304-49df-4626-a677-488295dc786d
function f(ψ::AbstractVector)
		L = length(ψ)
		D = Int((-(3/2) + sqrt((9/4) + 2L)) )

		Λμ = ψ[1:D]
		A = diagm(exp.(1e-3 .+ ψ[D+1:2*D]))
		L = A + CompressedLowerTriangular(D, 1, ψ[2*D+1:end])
		L⁻¹ = inv(L)

		Λ = L * L'
		μᵀΛμ = ((L⁻¹' * L⁻¹) * Λμ)' * Λμ
		logdetΛ = 2*sum(ψ[D+1:2*D])

		Λμ, μᵀΛμ, Λ, logdetΛ
	end

# ╔═╡ 4154a267-326d-495c-aa60-156560aa5f72
Qλ = 2

# ╔═╡ 4c69229d-a272-49eb-9dc1-484c53bad338
pλ = Normal(Qλ)

# ╔═╡ eebb75b8-8ece-4e26-972b-ed962ce1dd42
Qσ = 10

# ╔═╡ 3a303018-4560-4553-96f6-1600b0e1d6c6
pσ = Normal(Qσ)

# ╔═╡ b629f310-dbf7-4b6e-9c3d-199a31a1bfb5
Qπ = 2

# ╔═╡ 51d10c40-5dce-47e2-aeba-958beadda7b8
pπ = Normal(Qπ)

# ╔═╡ c352df19-b653-43f3-b490-046447cea85a
phonemodel = NormalModel((pλ=pλ, pσ=pσ, pπ=pπ))

# ╔═╡ c47dcaf4-049a-45e1-853f-d589283b2038
# The dimension of the parameter space is determined
# by the dimension of the observed space
K = Int(2D + (D^2 - D) / 2)

# ╔═╡ add6cb90-47bd-40e8-bbbb-34938714a79b
pL = JointNormalFixedCov(zeros(Qπ + 1, K))

# ╔═╡ 5415acd0-5e5f-40db-97eb-91b659a1ac5b
pS = JointNormalFixedCov(zeros(Qπ + 1, K))

# ╔═╡ b69c63d7-3f45-45d2-8278-05bc2fb5f3f1
md"""
### Hyper-parameters

| name          | type            | value|  description                  |
|:--------------|:----------------|:----:|:------------------------------|
| ``Q_\lambda`` | hyper-parameter | $Qλ  | language embedding dimension  |
| ``Q_\sigma``  | hyper-parameter | $Qσ  | speaker embedding dimension   |
| ``Q_\pi``     | hyper-parameter | $Qπ  | phone embedding dimension     |
| ``D``         | constant        | $D   | observed space dimension      |
| ``K``         | constant        | $K   | phone model natural parameter dimension|
| ``L`` 		| constant 		  | $L   | number of languages in the data |
| ``S``         | constant        | $S   | number of speakers in  the data |
| ``P``         | constant        | $P   | number of phones in the data  |
"""

# ╔═╡ 946af84f-a5f4-4780-a8c5-13d4acbb9a80
md"""
### Model parameters summary

```math
\begin{align}

\boldsymbol{\theta} &= \Big\{ \{\mathbf{L}_i\}, \{\mathbf{S}_j \}, \{\boldsymbol{\lambda}_l \}, \{\boldsymbol{\sigma}_s \}, \{ \boldsymbol{\pi}_p \} \Big\} \\
\{\mathbf{L}_i \}, \; &\forall i \in \{1, \dots, Q_\lambda \} \\
\{\mathbf{S}_j \}, \;  &\forall j \in \{1, \dots, Q_\sigma \} \\
\{\mathbf{\boldsymbol{\lambda}}_l \}, \; &\forall l \in \{1, \dots, L\} \\
\{\mathbf{\boldsymbol{\sigma}}_s \}, \; &\forall s \in \{1, \dots, S\} \\
\{\mathbf{\boldsymbol{\pi}_p} \}, \; &\forall p \in \{1, \dots, P\}
\end{align}
```
"""

# ╔═╡ f2f5e843-3505-421e-9c9c-71a564387bb2
# Complete model: hierarchical generalized subspace model.
hgsm = GSMPhonetic(f, pL, pS, phonemodel)

# ╔═╡ 1dd9bfa7-b7c9-4075-9c6c-2f33871d5e0a
md"""
### Variational posterior

```math
\begin{align}
	q(\boldsymbol{\theta}) &= \prod_{m}q(\mathbf{L}_m) \prod_{n} q(\mathbf{S}_n) \prod_{l} q(\boldsymbol{\lambda}_l) \prod_{s} q(\boldsymbol{\sigma}_s) \prod_{p} q(\boldsymbol{\pi}_p) \\
	q(\mathbf{L}_m) &= \delta(\mathbf{L}_m - \boldsymbol{\mu}^{L}_m) \\
	q(\mathbf{S}_n) &=  \delta(\mathbf{S}_n - \boldsymbol{\mu}^{S}_n) \\
	q(\boldsymbol{\lambda}_l) &= \mathcal{N}(\boldsymbol{\mu}^{\lambda}_l, \boldsymbol{\Sigma}^{\lambda}_l) \\
	q(\boldsymbol{\sigma}_s) &= \mathcal{N}(\boldsymbol{\mu}^{\sigma}_s, \boldsymbol{\Sigma}^{\sigma}_s) \\
		q(\boldsymbol{\pi}_p) &= \mathcal{N}(\boldsymbol{\mu}^{\pi}_p, \boldsymbol{\Sigma}^{\pi}_p)
\end{align}
```
"""

# ╔═╡ 6f1225eb-8df1-4991-8292-c0b5f2c45c10
q₀L = [JointNormalFixedCov(i == Qλ + 1 ? zeros(Qπ + 1, K) : randn(Qπ + 1, K) / K)
	   for i = 1:(Qλ+1)]

# ╔═╡ 349ef71a-55f8-4c9a-8216-d6fbce8d42e2
q₀S = [JointNormalFixedCov(j == Qσ + 1 ? zeros(Qπ + 1, K) : randn(Qπ + 1, K) / K)
	   for j = 1:(Qσ+1)]

# ╔═╡ 27fa667c-4a76-4fed-9afc-7482dc1339a0
q₀λ = [Normal(zeros(Qλ), Matrix(I, Qλ, Qλ)) for lang in languages]

# ╔═╡ 3a5cadd7-d15b-4da3-8a0b-9e7678514ff0
q₀σ = [Normal(zeros(Qσ), Matrix(I, Qσ, Qσ)) for spk in speakers]

# ╔═╡ e343d7f3-9a50-4952-b429-40b982b2b49f
q₀π = [Normal(zeros(Qπ), Matrix(I, Qπ, Qπ)) for phone in phones]

# ╔═╡ 42ccdf0d-8a5f-404d-8965-cccc7f85d164
q₀θ = (q₀L, q₀S, q₀λ, q₀σ, q₀π)

# ╔═╡ ffe22379-2f91-4b54-8f74-413ac6d2ce08
md"""### Minimum Divergence"""

# ╔═╡ 4d374801-703c-4fb0-b77d-bbb41fbfd54b
function renorm_hyperbases(qM, qh, ph)
	μh = map(μ, qh)
	stdh = stdparams.(qh, μh)

	μh₀ = sum(μh) / length(qh)
	m₀, Σ₀ = stdparams(ph, μh₀)
	L₀ = cholesky(Σ₀).L
	L₀⁻¹ = inv(L₀)

	new_qh = [Normal(L₀⁻¹ * (m - m₀), Symmetric(L₀⁻¹ * Σ * L₀⁻¹')) for (m, Σ) in stdh]

	B = vec.(getproperty.(qM, (:M)))
	dims = size(qM[1].M)
	M, b = hcat(B[1:end-1]...), B[end]
	new_M = reshape.(eachcol(M * L₀), dims...)
	new_b = reshape(b + M * m₀, dims...)

	new_qM = [JointNormalFixedCov(V) for V in [new_M..., new_b]]
	new_qM, new_qh
end

# ╔═╡ 546fd344-42e2-4881-9607-da99db000c33
function renorm_bases(qL, qS, qh, ph)
	μh = map(μ, qh)
	stdh = stdparams.(qh, μh)

	μh₀ = sum(μh) / length(qh)
	m₀, Σ₀ = stdparams(ph, μh₀)
	L₀ = cholesky(Σ₀).L
	L₀⁻¹ = inv(L₀)

	new_qh = [Normal(L₀⁻¹ * (m - m₀), Symmetric(L₀⁻¹ * Σ * L₀⁻¹')) for (m, Σ) in stdh]

	Ms = getproperty.(qL, (:M))
	new_Ms = map(Ms) do M
		W, b = M[1:end-1,:], M[end:end, :]
		vcat(L₀' * W, b + m₀' * W)
	end
	new_qL = [JointNormalFixedCov(V) for V in new_Ms]

	Ms = getproperty.(qS, (:M))
	new_Ms = map(Ms) do M
		W, b = M[1:end-1,:], M[end:end, :]
		vcat(L₀' * W, b + m₀' * W)
	end
	new_qS = [JointNormalFixedCov(V) for V in new_Ms]

	new_qL, new_qS, new_qh
end

# ╔═╡ 8b939634-e3f2-4a29-91de-4ed56787bc78
function minimum_divergence(hgsm, qθ)
	qL, qS, qλ, qσ, qπ = qθ

	qL2, qλ2 = renorm_hyperbases(qL, qλ, hgsm.model.pθ.pλ)
	qS2, qσ2 = renorm_hyperbases(qS, qσ, hgsm.model.pθ.pσ)
	qL2, qS2, qπ2 = renorm_bases(qL2, qS2, qπ, hgsm.model.pθ.pπ)

	qL2, qS2, qλ2, qσ2, qπ2
end

# ╔═╡ 85e89c27-0371-42a9-8b3c-4fc677d1f7a7
md"""
### Training

"""

# ╔═╡ 2822d23b-197d-4dcc-85d0-3df653fcf3b6
begin
	qθ_md = q₀θ

	kls_md = []
	loss_md = []
	for t in 1:100000
		μ = getstats(hgsm, qθ_md)
		(l, (∇μ,)) = withgradient(μ -> elbo(hgsm, groups, qθ_md, μ) / N, μ)

		qθ_md = newposterior(hgsm, qθ_md, ∇μ; lrate=1e-2)

		if t % 100 == 0
			qθ_md = minimum_divergence(hgsm, qθ_md)
		end

		push!(kls_md, kldiv(hgsm, qθ_md, getstats(hgsm, qθ_md)))
		push!(loss_md, l)
	end
	loss_md
end

# ╔═╡ 5162a02c-b803-4be5-bff8-801ec24001d3
plot(loss_md[1000:1:end], legend=false, linealpha=0.5)

# ╔═╡ 05c5717a-c0fa-47a4-9de9-c747b6fb3f93
md"""### Visualization"""

# ╔═╡ 4274c7f4-38d5-4bf0-87ae-417d715d9b3b
begin
	local qλ = qθ_md[3]
	local p = plot(legend=false)
	for (i, q) in enumerate(qλ)
		x, diagxxᵀ, hxxᵀ = unpack(q, μ(q))
		m = x
		Σ = diagm(diagxxᵀ) + CompressedSymmetric(Qλ, 1, hxxᵀ) - x * x'
		m = [m[1], m[2]]
		Σ = Σ[1:2,1:2]
		annotate!(m[1], m[2], id2lang[i])
		plotnormal!(m, Σ, σ=2, fillalpha=0.4, linealpha=0, label=false)
	end
	p
end

# ╔═╡ aab42b7b-b269-4e6f-af74-f09ec9c6b862
begin
	local qσ = qθ_md[4]
	local p = plot(size=(800,800), legend=false, xlims=(-3, 3), ylims=(-3, 3))
	for (i, q) in enumerate(qσ)
		x, diagxxᵀ, hxxᵀ = unpack(q, μ(q))
		m = x
		Σ = diagm(diagxxᵀ) + CompressedSymmetric(Qσ, 1, hxxᵀ) - x * x'
		m = [m[1], m[2]]
		Σ = [Σ[1,1] Σ[1,2]; Σ[2, 1] Σ[2,2]] #Σ[1:2,1:2]
		#annotate!(m[1], m[2], speakers[i][1])
		color = speakers[i][2] == "male" ? :blue : :green
		plotnormal!(m, Σ, σ=2, fillcolor=color, fillalpha=0.4, linealpha=0, label=false)
	end
	p
end

# ╔═╡ b9bba539-ba31-420e-a718-a952db80d4c2
begin
	local qπ = qθ_md[5]
	local d1, d2 = 1, 2
	local p = plot(size=(600, 500), legend=false, xlabel="π₁", ylabel="π₂")
	for (i, q) in enumerate(qπ)
		x, diagxxᵀ, hxxᵀ = unpack(q, μ(q))
		m = x
		Σ = diagm(diagxxᵀ) + CompressedSymmetric(Qπ, 1, hxxᵀ) - x * x'
		m = [m[d1], m[d2]]
		Σ = [Σ[d1, d1] Σ[d1, d2]; Σ[d2, d1] Σ[d2, d2]]
		annotate!(m[1], m[2], id2phone[i])
		plotnormal!(m, Σ, σ=2, fillcolor=:blue, fillalpha=0.1, linealpha=0, label=false)
	end
	p
end

# ╔═╡ f7d8e930-3d7b-4135-bb64-a8b17367f683
qL, qS, qλ, qσ, qπ = qθ_md

# ╔═╡ cd347af2-ad3e-4e7d-b471-54b668511e69
function hgsm_η(hgsm::GSMPhonetic, qL, qS, λ, σ, π)
	L = getproperty.(qL, (:M,))
	S = getproperty.(qS, (:M,))
	P = sum(L .* vcat(λ, 1)) + sum(S .* vcat(σ, 1))
	ψ = P' * vcat(π, 1)
	Λμ, μᵀΛμ, Λ, logdetΛ = hgsm.f(ψ)
	vcat(Λμ, -(1/2)*diag(Λ), -vech(Λ, 1))
end

# ╔═╡ 7d457821-67c0-47e0-9c00-7bd2a6560cc8
hgsm_η(hgsm, qL, qS, qλ[1].μ, qσ[1].μ, qπ[1].μ)

# ╔═╡ bffb14d2-1173-4511-80b2-0676af0cd3f8
function hgsm_μ(hgsm::GSMPhonetic, η)
	L = length(η)
	D = Int((-(3/2) + sqrt((9/4) + 2L)))

	Λμ = η[1:D]
	diagΛ = -2 * η[D+1:2*D]
	vechΛ = - η[2*D+1:end]
	Λ = diagm(diagΛ) + CompressedSymmetric(D, 1, vechΛ)

	Σ = inv(Symmetric(Λ))
	μ = Σ * Λμ
	xxᵀ = Σ + μ * μ'
	vcat(μ, diag(xxᵀ), vech(xxᵀ, 1))
end

# ╔═╡ 3a4921ab-d8e4-463b-ae41-d9ac4b034c41
hgsm_μ(hgsm, [0., 0., -1, -1, 0.5])

# ╔═╡ c3f6fc47-2843-498d-8d86-0c99e8638ff4
hgsm_μ(hgsm, hgsm_η(hgsm, qL, qS, qλ[1].μ, qσ[1].μ, qπ[1].μ))

# ╔═╡ 5d54df87-0d8a-4371-93e9-cc20fd659974
function hgsm_A(hgsm::GSMPhonetic, η)
	L = length(η)
	D = Int((-(3/2) + sqrt((9/4) + 2L)))

	Λμ = η[1:D]
	diagΛ = -2 * η[D+1:2*D]
	vechΛ = - η[2*D+1:end]
	Λ = diagm(diagΛ) + CompressedSymmetric(D, 1, vechΛ)
	Σ = inv(Λ)

	-(1/2) * logdet(Λ) + (1/2) * Λμ' * (Σ * Λμ)
end

# ╔═╡ f0d90bdc-aa48-4edf-a673-e6dc6d27e2e0
function trajectory_λ(hgsm, qL, qS, λ₁, λ₂, σ, π; step=1e-2)
	traj = [λ₁]
	dist = 0
	H = η -> ForwardDiff.hessian(η -> hgsm_A(hgsm, η), η)
	J = λ -> ForwardDiff.jacobian(λ -> hgsm_η(hgsm, qL, qS, λ, σ, π), λ)

	η₂ = hgsm_η(hgsm, qL, qS, λ₂, σ, π)

	while norm(traj[end] - λ₂) > 1e-1
		λ = traj[end]
		ηₜ = hgsm_η(hgsm, qL, qS, λ, σ, π)
		μₜ = hgsm_μ(hgsm, ηₜ)
		Jλ = J(λ)
		∇λ = (1/2) * Jλ' * (η₂ - ηₜ)
		#∇λ = -Jλ' * (ηₜ - (η₂ + μₜ))
		#∇λ = (λ₂ - λ)
		∇λ = ∇λ ./ norm(∇λ)

		G = Jλ' * Symmetric(H(ηₜ)) * Jλ
		v = ∇λ #inv(G) * ∇λ
		v = step * v ./ norm(v)
		dist += sqrt((step * ∇λ)' * G * (step * ∇λ))
		#dist += sqrt((step * ∇λ)' * G * (step * ∇λ))

		push!(traj, λ + v)
	end
	hcat(traj...), dist
end

# ╔═╡ 0d8c54e2-3677-4d5d-adc1-b6b3f482e409
t1, dist1 = trajectory_λ(hgsm, qL, qS, qλ[1].μ, qλ[6].μ, qσ[1].μ, qπ[1].μ)

# ╔═╡ b4c44e20-ffbc-4376-8cc6-3d9b1f64af9a
t2, dist2 = trajectory_λ(hgsm, qL, qS, qλ[6].μ, qλ[1].μ, qσ[1].μ, qπ[1].μ)

# ╔═╡ b317dbf6-2a60-400f-96d5-addd740be016
begin
	# Manual sorting so that languages from the same family appear together
	sorted_languages = [:spa, :ita, :fra, :rus, :pol, :deu]

	dist_lang = zeros(L, L)
	local idmap = Dict(i => findall(l -> l == lang, languages)[1] for (i, lang) in enumerate(sorted_languages))
	for i = 1:L, j = 1:L
		x, y = idmap[i], idmap[j]
		for phone in phones_per_lang[languages[x]]
			p = phone2id[phone]
			println(i, ": ", languages[x], " ", j, ": ", languages[y], " ", phones[p])
			traj, Δ = trajectory_λ(hgsm, qL, qS, qλ[x].μ, qλ[y].μ, zeros(Qσ), qπ[p].μ)
			dist_lang[i, j] += Δ / length(phones_per_lang[languages[x]])
		end
	end
end

# ╔═╡ 8d416ab7-3ae7-4984-a7e6-1f56692d1b7b
begin
	mdstr = "| language | # speakers | phones |\n"
	mdstr *= "|:---------|:-----------:|:-------|\n"
	for lang in sorted_languages
		mdstr *= "| $lang | "

		ns = 0
		for (l, spk) in speakers
			if l == lang ns += 1end
		end
		mdstr *= "$ns | "

		for phone in sort(collect(phones_per_lang[lang]))
			mdstr *= "$phone "
		end
		mdstr *= "|\n"
	end
	Markdown.parse(mdstr)
end

# ╔═╡ 69777c50-240b-4be7-a824-b5dd0eb5fb7c
languages

# ╔═╡ 00d323cd-b928-49cd-a590-0e8852ffb264
begin
	heatmap(
		dist_lang,
		xticks=(1:6, sorted_languages),
		yticks=(1:6, sorted_languages),
		tickfont = (14,),
		xmirror=true,
		yflip=true,
		c=:grays
	)
end

# ╔═╡ 4d44afe7-b0dc-40b1-a00d-73ad346e0576
dist_lang

# ╔═╡ cea584da-2d0e-4c4f-a964-2f7254bf74aa
languages

# ╔═╡ 295da023-8cb9-4b13-99f8-fa7690c20b5c
begin
	local d1, d2 = 1, 2
	local qλ = qθ_md[3]
	local p = plot(legend=false)
	for (i, q) in enumerate(qλ)
		x, diagxxᵀ, hxxᵀ = unpack(q, μ(q))
		m = x
		Σ = diagm(diagxxᵀ) + CompressedSymmetric(Qλ, 1, hxxᵀ) - x * x'
		m = [m[d1], m[d2]]
		Σ = [Σ[d1, d1] Σ[d1, d2]; Σ[d2, d1] Σ[d2, d2]]
		annotate!(m[1], m[2], id2lang[i])
		plot!(t1[d1, :], t1[d2, :], linecolor=:black)
		plot!(t2[d1, :], t2[d2, :], linecolor=:black)
		plotnormal!(m, Σ, σ=2, fillalpha=0.4, linealpha=0, label=false)
	end
	p
end

# ╔═╡ d15d2ae4-01b4-4312-b921-2fac2d943036
function hgsm_llh(hgsm::GSMPhonetic, X, qL, qS, λ, σ, π)
	L = getproperty.(qL, (:M,))
	S = getproperty.(qS, (:M,))
	P = sum(L .* vcat(λ, 1)) + sum(S .* vcat(σ, 1))
	ψ = P' * vcat(π, 1)
	#Λμ, μᵀΛμ, Λ, logdetΛ = hgsm.f(ψ)
	#vcat(Λμ, -(1/2)*diag(Λ), -vech(Λ, 1))
	loglikelihood(hgsm.model, X, hgsm.f(ψ))
end

# ╔═╡ fa30a558-5bed-4ce1-8be2-ed55a3d7a9f8
function hgsm_stdparams(hgsm::GSMPhonetic, qL, qS, λ, σ, π)
	L = getproperty.(qL, (:M,))
	S = getproperty.(qS, (:M,))
	P = sum(L .* vcat(λ, 1)) + sum(S .* vcat(σ, 1))
	ψ = P' * vcat(π, 1)
	Λμ, μᵀΛμ, Λ, logdetΛ = hgsm.f(ψ)
	Σ = inv(Λ)
	μ = Σ * Λμ
	Σ, μ
end

# ╔═╡ a36c7fd9-6e47-4386-a9e0-cb9a9dc6a644
id2spk

# ╔═╡ 0f26901f-e592-4768-ac67-ee526a640ede
phone2id

# ╔═╡ 8509deeb-f156-4350-a21a-d1d70949091d
lang2id

# ╔═╡ 696d1d7d-68df-45c3-b5d9-0cae311e553a
repeat([2], 10)

# ╔═╡ 6836772d-67c1-475f-8d99-e11fc7f752a0
begin
	local rawdata = readdlm(localfile)
	local data = []

	#pdataraw = Dict("œ" => [], "ɛ" => [])
	#pdataraw = Dict("ã" => [], "õ" => [])
	#pdataraw = Dict("ẽ" => [], "õ" => [])
	pdataraw = Dict("$p" => [] for p in phones)
	testphones = collect(keys(pdataraw))
	for (lang, spk, sex, ipaphone, f1, f2, f3, duration) in eachrow(rawdata)
		# Remove the duration symbol
		ipaphone = replace(ipaphone, "ː" => "")

		lang, spk, ipaphone = lang, spk, ipaphone

		#if lang == "fra" && spk == "spk8" && sex == "male" && ipaphone in keys(pdataraw)
		if ipaphone in keys(pdataraw)
			push!(pdataraw[ipaphone], vcat(f1, f2, f3))
		end
		push!(data, vcat(f1, f2, f3))
	end

	pdata = Dict()
	for (k, v) in pdataraw
		pdata[k] = hcat(v...)
	end
	data = hcat(data...)

	local x̄ = sum(data, dims=2) / size(data, 2)
	local x̄² = sum(data .^ 2, dims=2) / size(data, 2)
	local var = x̄² - (x̄ .^ 2)
	for (k, v) in pdata
		pdata[k] = (v .- x̄) ./ sqrt.(var)
	end
	pdata
end

# ╔═╡ 9c9d2747-b6db-4957-84be-3c327c4dc1df
begin
	local spk = 3
	local lang = 4
	#local testphones = ["ẽ", "õ"] #["œ", "ɛ"]
	p = plot(xlims=(-3, 3), ylims=(-3, 3))
	for phone in testphones
		i = phone2id[Symbol(phone)]
		Σ, m = hgsm_stdparams(hgsm, qL, qS, qλ[lang].μ, qσ[spk].μ, qπ[i].μ)
		plotnormal!(m[1:2], Σ[1:2,1:2], σ=2, fillalpha=0.4, linealpha=0, label="$phone")
	end

	for phone in testphones
		scatter!(pdata["$phone"][1,1:10], pdata["$phone"][2,:], label="$phone")
	end
	p
end

# ╔═╡ d711ad75-3fc6-4726-a506-962aaccb2de7
begin
	local p = plot()
	for phone in testphones
		scatter!(pdata["$phone"][1,1:10], pdata["$phone"][2,1:10], label="$phone")
	end
	p
end

# ╔═╡ c80d9947-e02e-4fb0-bcb3-e29a4702688e
begin
	local spk = 3
	local i = phone2id[Symbol(testphones[1])]
	local j = phone2id[Symbol(testphones[2])]
	local errors = zeros(L)

	for lang in 1:L

		llh1 = vcat(
			hgsm_llh(hgsm, pdata[testphones[1]], qL, qS, qλ[lang].μ, qσ[spk].μ, qπ[i].μ),
			hgsm_llh(hgsm, pdata[testphones[1]], qL, qS, qλ[lang].μ, qσ[spk].μ, qπ[j].μ)
		)

		llh2 = vcat(
			hgsm_llh(hgsm, pdata[testphones[2]], qL, qS, qλ[lang].μ, qσ[spk].μ, qπ[i].μ),
			hgsm_llh(hgsm, pdata[testphones[2]], qL, qS, qλ[lang].μ, qσ[spk].μ, qπ[j].μ)
		)
		errors[lang] = ( sum(llh1[1,:] .> llh1[2,:]) + sum(llh2[1, :] .< llh2[2, :]) ) / (size(llh1, 2) + size(llh2, 2))
	end
	errors
end

# ╔═╡ 43692f77-f402-41c7-b22e-0b8a42bf8451
speakers

# ╔═╡ 6d1e2c27-4c08-4324-af85-45ed1c4950e0
local rawdata = readdlm(localfile)

# ╔═╡ Cell order:
# ╟─6958bde2-c50f-4afd-8d08-7dfbfd97169a
# ╠═3cd07b90-9e21-4a74-8d99-a6375beb7d1e
# ╟─b83311ab-ca4f-432e-8627-5207517a81db
# ╠═d1588a1c-babe-4223-bc6d-f353314ed2c0
# ╠═262e0395-814f-4315-b25e-21c18a7a9a5f
# ╟─4bbdd3d2-eb62-4889-9c8a-20141c0b140a
# ╠═ee1d1c45-70f5-4fe7-a8e2-9134d45a2b95
# ╟─474abeb1-82f8-4a74-a03c-be1cf55df840
# ╠═31e56fed-e4f0-4c71-abb2-5c2eff5a9a47
# ╠═61cc1ab2-a615-4a9b-91f9-bc5909e29297
# ╠═e08cc384-49dc-4585-a2ab-fbec9fcfe9d2
# ╠═b8e7f92f-0065-46fe-ac6b-5adf29b3379f
# ╠═ef71661f-338e-4d57-86ed-fc2c308a744f
# ╟─c01de0cd-43f8-49cc-a8c9-ffbbc182d9e9
# ╟─8d416ab7-3ae7-4984-a7e6-1f56692d1b7b
# ╟─8b535f17-7a55-42da-ba7a-c161bec9fd28
# ╠═4c69229d-a272-49eb-9dc1-484c53bad338
# ╠═3a303018-4560-4553-96f6-1600b0e1d6c6
# ╠═51d10c40-5dce-47e2-aeba-958beadda7b8
# ╟─5175712e-7f27-459e-8e6a-9bed325d6b4e
# ╠═add6cb90-47bd-40e8-bbbb-34938714a79b
# ╠═5415acd0-5e5f-40db-97eb-91b659a1ac5b
# ╟─763b5f95-db02-4ede-b632-848a646232d9
# ╠═c352df19-b653-43f3-b490-046447cea85a
# ╟─841a3f0e-8975-49d2-959e-afe832ed41e4
# ╠═b0761304-49df-4626-a677-488295dc786d
# ╟─b69c63d7-3f45-45d2-8278-05bc2fb5f3f1
# ╠═4154a267-326d-495c-aa60-156560aa5f72
# ╠═eebb75b8-8ece-4e26-972b-ed962ce1dd42
# ╠═b629f310-dbf7-4b6e-9c3d-199a31a1bfb5
# ╠═c47dcaf4-049a-45e1-853f-d589283b2038
# ╟─946af84f-a5f4-4780-a8c5-13d4acbb9a80
# ╠═f2f5e843-3505-421e-9c9c-71a564387bb2
# ╟─1dd9bfa7-b7c9-4075-9c6c-2f33871d5e0a
# ╠═6f1225eb-8df1-4991-8292-c0b5f2c45c10
# ╠═349ef71a-55f8-4c9a-8216-d6fbce8d42e2
# ╠═27fa667c-4a76-4fed-9afc-7482dc1339a0
# ╠═3a5cadd7-d15b-4da3-8a0b-9e7678514ff0
# ╠═e343d7f3-9a50-4952-b429-40b982b2b49f
# ╠═42ccdf0d-8a5f-404d-8965-cccc7f85d164
# ╟─ffe22379-2f91-4b54-8f74-413ac6d2ce08
# ╠═4d374801-703c-4fb0-b77d-bbb41fbfd54b
# ╠═546fd344-42e2-4881-9607-da99db000c33
# ╠═8b939634-e3f2-4a29-91de-4ed56787bc78
# ╟─85e89c27-0371-42a9-8b3c-4fc677d1f7a7
# ╠═2822d23b-197d-4dcc-85d0-3df653fcf3b6
# ╠═5162a02c-b803-4be5-bff8-801ec24001d3
# ╟─05c5717a-c0fa-47a4-9de9-c747b6fb3f93
# ╠═4274c7f4-38d5-4bf0-87ae-417d715d9b3b
# ╠═aab42b7b-b269-4e6f-af74-f09ec9c6b862
# ╠═b9bba539-ba31-420e-a718-a952db80d4c2
# ╠═f7d8e930-3d7b-4135-bb64-a8b17367f683
# ╠═cd347af2-ad3e-4e7d-b471-54b668511e69
# ╠═7d457821-67c0-47e0-9c00-7bd2a6560cc8
# ╠═bffb14d2-1173-4511-80b2-0676af0cd3f8
# ╠═3a4921ab-d8e4-463b-ae41-d9ac4b034c41
# ╠═c3f6fc47-2843-498d-8d86-0c99e8638ff4
# ╠═5d54df87-0d8a-4371-93e9-cc20fd659974
# ╠═f0d90bdc-aa48-4edf-a673-e6dc6d27e2e0
# ╠═0d8c54e2-3677-4d5d-adc1-b6b3f482e409
# ╠═b4c44e20-ffbc-4376-8cc6-3d9b1f64af9a
# ╠═b317dbf6-2a60-400f-96d5-addd740be016
# ╠═69777c50-240b-4be7-a824-b5dd0eb5fb7c
# ╠═00d323cd-b928-49cd-a590-0e8852ffb264
# ╠═4d44afe7-b0dc-40b1-a00d-73ad346e0576
# ╠═cea584da-2d0e-4c4f-a964-2f7254bf74aa
# ╠═295da023-8cb9-4b13-99f8-fa7690c20b5c
# ╠═d15d2ae4-01b4-4312-b921-2fac2d943036
# ╠═fa30a558-5bed-4ce1-8be2-ed55a3d7a9f8
# ╠═9c9d2747-b6db-4957-84be-3c327c4dc1df
# ╠═d711ad75-3fc6-4726-a506-962aaccb2de7
# ╠═c80d9947-e02e-4fb0-bcb3-e29a4702688e
# ╠═a36c7fd9-6e47-4386-a9e0-cb9a9dc6a644
# ╠═0f26901f-e592-4768-ac67-ee526a640ede
# ╠═8509deeb-f156-4350-a21a-d1d70949091d
# ╠═696d1d7d-68df-45c3-b5d9-0cae311e553a
# ╠═6836772d-67c1-475f-8d99-e11fc7f752a0
# ╠═43692f77-f402-41c7-b22e-0b8a42bf8451
# ╠═6d1e2c27-4c08-4324-af85-45ed1c4950e0

max_stars_repo_path: examples/Learning and Modeling Languages Phonetics v2.jl

max_stars_repo_name: BUTSpeechFIT/PPCA

max_issues_repo_name: BUTSpeechFIT/PPCA

max_forks_repo_name: BUTSpeechFIT/PPCA


text:
module O_structs

export ms_struct,mp_struct,pp_struct,ps_struct,hc_struct,sd_struct,al_params

mutable struct ms_struct{UI<:UnitRange{Int64}}
      P::UI
     Ii::UI
     Io::UI
end

mutable struct mp_struct{AF1<:Array{Float64,1},AF2<:Array{Float64,2}}
      xh::AF2
      xm::AF1
      ci::AF2
      cf::AF1
     map::Dict{Tuple{Int64,Int64},Array{Int64,1}}
end

mutable struct pp_struct{AF1<:Array{Float64,1},AF2<:Array{Float64,2},AF3<:Array{Float64,3},F<:Float64}
   cvg::AF1
   cfg::AF1
    eg::AF1
    ηg::AF1
    rg::AF1
    ηb::AF1
    pb::AF1
    pr::AF3
    pd::AF2
    cs::F
end

mutable struct ps_struct{UI<:UnitRange{Int64},I<:Int64}
    S::UI
    H::UI
    P::UI
    G::UI
    B::UI
    R::UI
   b0::I
   r0::I
end

mutable struct hc_struct{I<:Int64,AF2<:Array{Float64,2}}
    h::AF2
    c::AF2
   nx::I
   nc::I
  nx0::I
   ny::I
end

mutable struct al_params
   ϵ::Float64
   J::Int64
end

mutable struct sd_struct{AI1<:Array{Int64,1}}
   yrs::AI1
   yri::AI1
    bi::AI1
 techs::Array{String,1}
end

end

max_stars_repo_name: nimazzi/PlanModelDet

max_issues_repo_name: nimazzi/PlanModelDet

max_forks_repo_name: nimazzi/PlanModelDet


text:
include("euler/euler.jl")
using .Primes: prime_numbers

function compute()::Int
    for prime ∈ reverse(prime_numbers(7_654_321))
        str_prime = string(prime)
        if sort(parse.(Int, ∪(str_prime))) == collect(1:length(str_prime))
            return prime
        end
    end
end

max_stars_repo_name: Dynortice/Project-Euler

max_issues_repo_name: Dynortice/Project-Euler

max_forks_repo_name: Dynortice/Project-Euler