text:
##############################################################################
##
## Define a type that iterates through q-grams of a string
##
############################################################################
struct QGramIterator{S <: AbstractString}
	s::S # string
	q::Int # Length of Qgram
end

function Base.iterate(qgram::QGramIterator,
	state = (1, nextind(qgram.s, 0, qgram.q)))
	istart, iend = state
	iend > ncodeunits(qgram.s) && return nothing
	element = SubString(qgram.s, istart, iend)
	nextstate = nextind(qgram.s, istart), nextind(qgram.s, iend)
	element, nextstate
end
Base.length(qgram::QGramIterator) = max(length(qgram.s) - qgram.q + 1, 0)
Base.eltype(qgram::QGramIterator{SubString{S}}) where {S} = SubString{S}
Base.eltype(qgram::QGramIterator{S}) where {S} = SubString{S}

"""
Return an iterator that iterates on the QGram of the string

### Arguments
* `s::AbstractString`
* `q::Integer`: length of qgram

## Examples
```julia
using StringDistances
for x in qgram("hello", 2)
	@show x
end
```
"""
function qgram(s::AbstractString, q::Integer)
	QGramIterator{typeof(s)}(s, q)
end

##############################################################################
##
##
##
##
##
##############################################################################
"""
   count_map(x1, x2)

For two iterators `x1` and `x2`, `count_map(x1, x2)`  returns an dictionary
that returns,  for each element in `x1` or `x2`, a tuple with the numbers of
times it appears in `x1` and the number of times it appears in `x2`
"""
function count_map(s1, s2)
	K = promote_type(eltype(s1), eltype(s2))
	d = Dict{K, Tuple{Int, Int}}()
	# I use a faster way to change a dictionary key
	# see setindex! in https://github.com/JuliaLang/julia/blob/master/base/dict.jl#L380
	sizehint!(d, length(s1) + length(s2))
	for x1 in s1
		index = Base.ht_keyindex2!(d, x1)
		if index > 0
			d.age += 1
			@inbounds d.keys[index] = x1
			@inbounds d.vals[index] = (d.vals[index][1] + 1, 0)
		else
			@inbounds Base._setindex!(d, (1, 0), x1, -index)
		end
	end
	for x2 in s2
		index = Base.ht_keyindex2!(d, x2)
		if index > 0
			d.age += 1
			@inbounds d.keys[index] = x2
			@inbounds d.vals[index] = (d.vals[index][1], d.vals[index][2] + 1)
		else
			@inbounds Base._setindex!(d, (0, 1), x2, -index)
		end
	end
	return d
end

#= Trie
function count_map(s1, s2)
	d = Trie{Tuple{Int, Int}}()
	for ch1 in s1
		node = d
		for char in ch1
			if !haskey(node.children, char)
	            node.children[char] = Trie{Tuple{Int, Int}}()
	        end
	        node = node.children[char]
	    end
	    node.value = node.is_key ? (node.value[1] + 1, 0) : (1, 0)
	    node.is_key = true
	end
	for ch2 in s2
		node = d
		for char in ch2
			if !haskey(node.children, char)
	            node.children[char] = Trie{Tuple{Int, Int}}()
	        end
	        node = node.children[char]
	    end
	    node.value = node.is_key ? (node.value[1], node.value[2]+ 1) : (0, 1)
	    node.is_key = true
	end
	return iterator(d)
end
function iterator(t::Trie, found = Tuple{Int, Int}[])
    if t.is_key
    	t.is_key = false
    	push!(found, t.value)
    else
    	for k in values(t.children)
    		iterator(k, found)
    	end
    end
    return found
end
=#

##############################################################################
##
## Distance on strings is computed by set distance on qgram sets
##
##############################################################################
abstract type AbstractQGramDistance <: SemiMetric end

function evaluate(dist::AbstractQGramDistance, s1::AbstractString, s2::AbstractString)
	x = count_map(qgram(s1, dist.q), qgram(s2, dist.q))
	evaluate(dist, x)
end

##############################################################################
##
## q-gram
##
##############################################################################
"""
For an AbstractString s, denote v(s) the vector on the space of q-grams of length N, that contains the number of times a q-gram appears in s
The q-gram distance is ||v(s1) - v(s2)||
"""

"""
	QGram(q::Int)

Creates a QGram metric.

The distance corresponds to

``||v(s1, q) - v(s2, q)||``

where ``v(s, q)`` denotes the vector on the space of q-grams of length q, that contains the number of times a q-gram appears for the string s
"""
struct QGram <: AbstractQGramDistance
	q::Int
end

function evaluate(dist::QGram, count_dict)
	n = 0
	for (n1, n2) in values(count_dict)
		n += abs(n1 - n2)
	end
	n
end

##############################################################################
##
## cosine
##
##
##############################################################################
"""
	Cosine(q::Int)

Creates a Cosine metric.

The distance corresponds to

`` 1 - v(s1, q).v(s2, q)  / ||v(s1, q)|| * ||v(s2, q)||``

where ``v(s, q)`` denotes the vector on the space of q-grams of length q, that contains the number of times a q-gram appears for the string s
"""
struct Cosine <: AbstractQGramDistance
	q::Int
end

function evaluate(dist::Cosine, count_dict)
	norm1, norm2, prodnorm = 0, 0, 0
	for (n1, n2) in values(count_dict)
		norm1 += n1^2
		norm2 += n2^2
		prodnorm += n1 * n2
	end
	1.0 - prodnorm / (sqrt(norm1) * sqrt(norm2))
end

##############################################################################
##
## Jaccard
##
##############################################################################
"""
	Jaccard(q::Int)

Creates a Jaccard metric.

The distance corresponds to

``1 - |Q(s1, q) ∩ Q(s2, q)| / |Q(s1, q) ∪ Q(s2, q))|``

where ``Q(s, q)``  denotes the set of q-grams of length n for the string s
"""
struct Jaccard <: AbstractQGramDistance
	q::Int
end

function evaluate(dist::Jaccard, count_dict)
	ndistinct1, ndistinct2, nintersect = 0, 0, 0
	for (n1, n2) in values(count_dict)
		ndistinct1 += n1 > 0
		ndistinct2 += n2 > 0
		nintersect += (n1 > 0) & (n2 > 0)
	end
	1.0 - nintersect / (ndistinct1 + ndistinct2 - nintersect)
end

##############################################################################
##
## SorensenDice
##
##############################################################################
"""
	SorensenDice(q::Int)

Creates a SorensenDice metric

The distance corresponds to

``1 - 2 * |Q(s1, q) ∩ Q(s2, q)|  / (|Q(s1, q)| + |Q(s2, q))|)``

where ``Q(s, q)``  denotes the set of q-grams of length n for the string s
"""
struct SorensenDice <: AbstractQGramDistance
	q::Int
end

function evaluate(dist::SorensenDice, count_dict)
	ndistinct1, ndistinct2, nintersect = 0, 0, 0
	for (n1, n2) in values(count_dict)
		ndistinct1 += n1 > 0
		ndistinct2 += n2 > 0
		nintersect += (n1 > 0) & (n2 > 0)
	end
	1.0 - 2.0 * nintersect / (ndistinct1 + ndistinct2)
end

##############################################################################
##
## overlap
##
## 1 -  |intersect(Q(s1, q), Q(s2, q))| / min(|Q(s1, q)|, |Q(s2, q)))
##############################################################################
"""
	Overlap(q::Int)

Creates a Overlap metric

The distance corresponds to

``1 - |Q(s1, q) ∩ Q(s2, q)|  / min(|Q(s1, q)|, |Q(s2, q)|)``

where ``Q(s, q)``  denotes the set of q-grams of length n for the string s
"""
struct Overlap <: AbstractQGramDistance
	q::Int
end

function evaluate(dist::Overlap, count_dict)
	ndistinct1, ndistinct2, nintersect = 0, 0, 0
	for (n1, n2) in values(count_dict)
		ndistinct1 += n1 > 0
		ndistinct2 += n2 > 0
		nintersect += (n1 > 0) & (n2 > 0)
	end
	1.0 - nintersect / min(ndistinct1, ndistinct2)
end


text:
#https://www.mathworks.com/help/matlab/ref/geoplot.html
function geoplot(lat,lon;options=Dict())
    MATLAB.@mput lat
    MATLAB.@mput lon
    optionsString = namevaluepairs(options)
    MATLAB.eval_string("geoplot(lat,lon"*optionsString*")")
end

function geoplot(lat,lon,lineSpec;options=Dict())
    MATLAB.@mput lat
    MATLAB.@mput lon
    optionsString = namevaluepairs(options)
    MATLAB.eval_string("geoplot(lat,lon"*",\"$lineSpec\""*optionsString*")")
end

ext: jl


text:
using QuantumComputing

psi = Ket(0,0,0,0)
n = length(psi)

println("input:")
println(psi)

psi = Hadamard(n,4)*GCNot(n,1,2,3)*Not(n,3)*Not(n,2)*psi

println("output:")
println(psi)


text:
struct Buoyancy{M, G}
                   model :: M
    vertical_unit_vector :: G
end

struct ZDirection end

"""
    Buoyancy(; model, vertical_unit_vector=ZDirection())

Uses a given buoyancy `model` to create buoyancy in a model. The optional keyword argument
`vertical_unit_vector` can be used to specify the direction opposite to the gravitational
acceleration (which we take here to mean the "vertical" direction).

Example
=======

```julia
using Oceananigans

grid = RegularRectilinearGrid(size=(1, 8, 8), extent=(1, 1000, 100))
θ = 45 # degrees
g̃ = (0, sind(θ), cosd(θ))

buoyancy = Buoyancy(model=BuoyancyTracer(), vertical_unit_vector=g̃)

model = NonhydrostaticModel(grid=grid, buoyancy=buoyancy, tracers=:b)
```
"""
function Buoyancy(; model, vertical_unit_vector=ZDirection())
    vertical_unit_vector = validate_vertical_unit_vector(vertical_unit_vector)
    return Buoyancy(model, vertical_unit_vector)
end

validate_vertical_unit_vector(ĝ::ZDirection) = ĝ

function validate_vertical_unit_vector(ĝ)
    length(ĝ) == 3 || throw(ArgumentError("`vertical_unit_vector` must have length 3"))

    gx, gy, gz = ĝ

    gx^2 + gy^2 + gz^2 ≈ 1 ||
        throw(ArgumentError("`vertical_unit_vector` must be a unit vector with g[1]² + g[2]² + g[3]² ≈ 1"))

    return tuple(ĝ...)
end

@inline ĝ_x(buoyancy) = @inbounds buoyancy.vertical_unit_vector[1]
@inline ĝ_y(buoyancy) = @inbounds buoyancy.vertical_unit_vector[2]
@inline ĝ_z(buoyancy) = @inbounds buoyancy.vertical_unit_vector[3]

@inline ĝ_x(::Buoyancy{M, ZDirection}) where M = 0
@inline ĝ_y(::Buoyancy{M, ZDirection}) where M = 0
@inline ĝ_z(::Buoyancy{M, ZDirection}) where M = 1

#####
##### For convenience
#####

@inline required_tracers(bm::Buoyancy) = required_tracers(bm.model)

@inline get_temperature_and_salinity(bm::Buoyancy, C) = get_temperature_and_salinity(bm.model, C)

@inline ∂x_b(i, j, k, grid, b::Buoyancy, C) = ∂x_b(i, j, k, grid, b.model, C)
@inline ∂y_b(i, j, k, grid, b::Buoyancy, C) = ∂y_b(i, j, k, grid, b.model, C)
@inline ∂z_b(i, j, k, grid, b::Buoyancy, C) = ∂z_b(i, j, k, grid, b.model, C)

@inline top_buoyancy_flux(i, j, grid, b::Buoyancy, args...) = top_buoyancy_flux(i, j, grid, b.model, args...)

regularize_buoyancy(b) = b
regularize_buoyancy(b::AbstractBuoyancyModel) = Buoyancy(model=b)


text:
src = "../src"
if !in(src,LOAD_PATH)
    push!(LOAD_PATH, src)
end
using HTSLIB
using Base.Test
using Logging
@Logging.configure(level=DEBUG)

bios = open("data/100.bam","rb","bam")
handle = bios.handle
@debug("handle is $handle")
phts = convert(Ptr{HTSLIB.HTSFile},handle)
@debug("phts is $phts")
htsfl  = unsafe_load(phts)
@show htsfl
@debug("htsFile fn is $(htsfl.fn)")
@show bytestring(htsfl.fn)

#=
for i=0:150
    pkstr = convert(Ptr{HTSLIB.KString}, phts+i)
    kstr = unsafe_load(pkstr)
    @show i,kstr,Int128(kstr.l)
    if i >= 33
        segm = [39,40,41,42,43,55,56,60,64,72,78,79,80,81,100,101,102,103,104]
        if i in segm
            println("$(string(segm)) segment fault")
            continue
        end
        str = try
            HTSLIB.strptr(kstr)
           catch ex
            @show ex
            continue
        end
        info("Great success: $str")
        break
    end
end
#@show HTSLIB.strptr(kstr)
#@show htsfl
=#


text:
# author: aliyyousuf@gmail.com

# assignment of week_03.

data_to_fit = readdlm("Week3_PR_Data.dat",',')

# Converting data into Float64

data_to_fit_2 = zeros(Float64, 15,2)
for r = 1:15
    for c = 1
        data_to_fit_2[r,c] = parse(Float64,split(data_to_fit[r,c])[1])
        data_to_fit_2[r,2] = parse(Float64,split(data_to_fit[r,c])[2])
    end
end

# "For loop" to show all the values of data:

for r = 1:15
    for c = 1:2
        #println(data_to_fit_2[r,c])
    end
end

# assign the values from first column to variable x
# and second column to variable y.

x = data_to_fit_2[:,1]
y = data_to_fit_2[:,2]

using Plots
pyplot()
plot(x,y, title = "Model versus Data",label = "Actual values from given data",
xlabel = " x values",
ylabel = "y values",
line =:scatter)
gui()

function parab_fit(x)
    return a * (x)^2 + b * x + c
end

a = -0.789
b = 2.988
c = -1.5
plot!(parab_fit,title = "Model versus Data",-10,10,
label = "Model values")
gui()

max_stars_repo_name: aliyyousuff/MOOC

max_issues_repo_name: aliyyousuff/MOOC

max_forks_repo_name: aliyyousuff/MOOC


text:
#####
#
# Example script for TintiNet.jl
# Model IGBT (InceptGOR-Bert): 3-state SS Regressor
#
# To run training, revise your code and uncomment the last line.
# On slower GPUs, an epoch may take several minutes with the default configuration.
# Adjust model size using the `hyperpar` and `config` dictionaries.
#
#####

#####
# 0. Loading Packages and defining Hyperparameters
#####

using CUDA, Zygote, Flux, Transformers, JSON, BSON, Random, Dates, Printf
using IterTools: ncycle;
using Transformers.Basic, Transformers.Pretrain, Transformers.Datasets, Transformers.Stacks, Transformers.BidirectionalEncoder
using BSON: @load, @save
using Flux: @functor, onehot, onecold, mse, shuffle, sigmoid, gradient, unsqueeze
using Flux.Optimise: update!
using Statistics: mean, std

using TintiNet
dformat = Dates.DateFormat("yyyy-m-d-HH-MM-SS")
Random.seed!(102528)

##
# Algorithm Hyperparameters

hyperpar = Dict(
    "max_seq_len"       => 128,     # Maximum sequence length on the model (default 192)
    "batch_size"        => 128,     # Batch size for training (default 128)
    "number_of_epochs"  => 200,     # Number of epochs for training (10 for the example file)
    "learning_rate"     => 1e-3,    # Learning rate for optimizer (default 0.001)
    "smooth_epsilon"    => 1e-6     # Small value to smooth onehot result before logcrossentropy (default 1e-6)
)

##
# Model Configuration

const config = Dict(
  "input_vocab_size"                => 22,      # Size of the SEQUENCE/FASTA vocabulary, usually 20 aminoacids + "X" + "-"
  "output_size"                     => 5,       # Output size. Usually 5 due to sin/cos phi, sin/cos psi and acc.
  "embed_dimension"                 => 128,     # Size of the individual Embed vectors
  "bert_num_hidden_layers"          => 2,       # Number of consecutive Transformer layers on the BERT encoder
  "bert_transformer_layer_heads"    => 8,       # Number of heads in each Transformer layer
  "bert_transformer_layer_hidden"   => 128,     # Hidden size of the positionwise linear layer in each Transformer layer
  "bert_transformer_dropout_prob"   => 0.1,     # Dropout probability for the Transformer layers
  "bert_transformer_activation"     => gelu,    # Activation function for each Transformer layer (default: gelu)
  "detection_num_filters"           => 32       # Number of filters in each convolutional decoder layer
)

#####
# 2. Obtaining and shaping the data
#####

## 2.1. Training sequences, structures and lengths

remapnan(x) = x == nothing ? missing : x

training_set_datapath = "../data/fold01training_dataset.json"

open(training_set_datapath, "r") do file

    entries = JSON.parse(file)
    num_entries = length(entries)

    parsed_sequences = Vector{Vector{String}}(undef, num_entries)
    parsed_phis = Vector{Vector{Union{Missing, Float64}}}(undef, num_entries)
    parsed_psis = Vector{Vector{Union{Missing, Float64}}}(undef, num_entries)
    parsed_accs = Vector{Vector{Union{Missing, Float64}}}(undef, num_entries)
    parsed_lengths = Vector{Int}(undef, num_entries)

    for (i,j) in enumerate(entries)
        parsed_sequences[i] = string.(j["fasta_seq"])
        parsed_phis[i] = map(remapnan, j["dssp_phi"])
        parsed_psis[i] = map(remapnan, j["dssp_psi"])
        parsed_accs[i] = map(remapnan, j["dssp_acc"])
        parsed_lengths[i] = length(parsed_sequences[i])
    end

    order_idx = sortperm(parsed_lengths)

    global train_sequences_ordered = parsed_sequences[order_idx]
    global train_phis_ordered = parsed_phis[order_idx]
    global train_psis_ordered = parsed_psis[order_idx]
    global train_accs_ordered = parsed_accs[order_idx]
    global train_lengths_ordered = parsed_lengths[order_idx]

end

## 2.2. Testing sequences, structures and lengths

testing_set_datapath = "../data/fold01testing_dataset.json"

open(testing_set_datapath, "r") do file

    entries = JSON.parse(file)
    num_entries = length(entries)

    parsed_sequences = Vector{Vector{String}}(undef, num_entries)
    parsed_phis = Vector{Vector{Union{Missing, Float64}}}(undef, num_entries)
    parsed_psis = Vector{Vector{Union{Missing, Float64}}}(undef, num_entries)
    parsed_accs = Vector{Vector{Union{Missing, Float64}}}(undef, num_entries)
    parsed_lengths = Vector{Int}(undef, num_entries)

    for (i,j) in enumerate(entries)
        parsed_sequences[i] = string.(j["fasta_seq"])
        parsed_phis[i] = map(remapnan, j["dssp_phi"])
        parsed_psis[i] = map(remapnan, j["dssp_psi"])
        parsed_accs[i] = map(remapnan, j["dssp_acc"])
        parsed_lengths[i] = length(parsed_sequences[i])
    end

    order_idx = sortperm(parsed_lengths)

    global test_sequences_ordered = parsed_sequences[order_idx]
    global test_phis_ordered = parsed_phis[order_idx]
    global test_psis_ordered = parsed_psis[order_idx]
    global test_accs_ordered = parsed_accs[order_idx]
    global test_lengths_ordered = parsed_lengths[order_idx]

end

## 2.3. Lazy DataLoaders

training_dataset = Flux.Data.DataLoader((train_sequences_ordered, train_phis_ordered, train_psis_ordered, train_accs_ordered, train_lengths_ordered); batchsize=hyperpar["batch_size"], shuffle=true)
eval_train_dataset = Flux.Data.DataLoader((train_sequences_ordered, train_phis_ordered, train_psis_ordered, train_accs_ordered, train_lengths_ordered); batchsize=64,shuffle=false)
eval_test_dataset = Flux.Data.DataLoader((test_sequences_ordered, test_phis_ordered, test_psis_ordered, test_accs_ordered, test_lengths_ordered); batchsize=64,shuffle=false)

#####
# 3. Network Structure
#####

model = Stack(
        @nntopo(x → enc → enc → enc → inc → inc → inc → inc → inc → inc → inc → inc → inc → inc → pe:(inc, pe) → bert_input → bert_output → bert_output → bert_output → y → y → y → y → y),
         # Embedding Layers
        Embed(config["embed_dimension"], config["input_vocab_size"]),
        unsqueeze(3),
        x -> permutedims(x, [3,2,1,4]),
        # Processing Layers
        Inception8(config["embed_dimension"]),
        Dropout(config["bert_transformer_dropout_prob"]),
        Inception8(config["embed_dimension"]),
        Dropout(config["bert_transformer_dropout_prob"]),
        Inception8(config["embed_dimension"]),
        Dropout(config["bert_transformer_dropout_prob"]),
        Inception8(config["embed_dimension"]),
        Dropout(config["bert_transformer_dropout_prob"]),
        x -> permutedims(x, [3,2,1,4]),
        x -> reshape(x, size(x,1), size(x,2), size(x,4)),
        PositionEmbedding(config["embed_dimension"], hyperpar["max_seq_len"]; trainable = false),
        (e, pe) -> e .+ pe,
        # Processing Layers
        Bert(
            config["embed_dimension"],
            config["bert_transformer_layer_heads"],
            config["bert_transformer_layer_hidden"],
            config["bert_num_hidden_layers"];
            act = config["bert_transformer_activation"],
            pdrop = config["bert_transformer_dropout_prob"],
            attn_pdrop = config["bert_transformer_dropout_prob"]
        ),
        # Detection Layers
        unsqueeze(3),
        bert_output -> permutedims(bert_output, [3,2,1,4]),
        Conv((1,7), config["embed_dimension"] => config["detection_num_filters"], relu; pad = (0,3)),
        Conv((1,7), config["detection_num_filters"] => config["detection_num_filters"], relu; pad = (0,3)),
        Conv((1,7), config["detection_num_filters"] => config["output_size"], tanh; pad = (0,3)),
        x -> permutedims(x, [3,2,1,4]),
        x -> reshape(x, size(x,1), size(x,2), size(x,4)),
    ) |> gpu

#####
# 4. Model training
#####

opt = Flux.Optimise.ADAM(hyperpar["learning_rate"])     # Initializes the ADAM optimizer
ps = params(model)      # Collects model parameters

function regressor_train_loss(model, x, y, mask)
    sum( ( (model(x) .- y) .^ 2 ) .* mask ) / sum(mask)
end

function regressor_eval_loss(predictions, y, mask)
    sum( ( (predictions .- y) .^ 2 ) .* mask ) / sum(mask)
end

function build_regressor_target_array(
    phis::Vector{Vector{Union{Missing, Float64}}},
    psis::Vector{Vector{Union{Missing, Float64}}},
    accs::Vector{Vector{Union{Missing, Float64}}},
    seqlen::Int64)

    y = zeros(5, seqlen, length(phis))

    for seq_i in 1:length(phis)
        for aa_j in 1:length(phis[seq_i])
            if !ismissing(phis[seq_i][aa_j])
                y[1, aa_j, seq_i] = sind(phis[seq_i][aa_j])
                y[2, aa_j, seq_i] = cosd(phis[seq_i][aa_j])
            end

            if !ismissing(psis[seq_i][aa_j])
                y[3, aa_j, seq_i] = sind(psis[seq_i][aa_j])
                y[4, aa_j, seq_i] = cosd(psis[seq_i][aa_j])
            end

            if !ismissing(accs[seq_i][aa_j])
                y[5, aa_j, seq_i] = (accs[seq_i][aa_j] - 100) / 100
            end

        end
    end

    return(y)

end

function prepare_regressor_mask(batchofseqs::Vector{Vector{String}},
                                phis::Vector{Vector{Union{Missing, Float64}}},
                                psis::Vector{Vector{Union{Missing, Float64}}},
                                accs::Vector{Vector{Union{Missing, Float64}}},
                                seqlen::Int64)

    mask = getmask(preprocess(batchofseqs, seqlen)) # Uses getmask() from Transformers.Basic

    for i in 1:length(batchofseqs)
        for j in 1:length(batchofseqs[i])
            if ( ismissing(phis[i][j]) | ismissing(psis[i][j]) | ismissing(accs[i][j]) )
                mask[1, j, i] = 0.0
            end
        end
    end

    return(mask)

end

#define training loop
function train!(hyperpar, training_dataset,
                model, ps, opt,
                train_loss=regressor_train_loss, eval_loss=regressor_eval_loss,
                eval_train_dataset=eval_train_dataset, eval_test_dataset=eval_test_dataset,
                fastaVocab=fastaVocab)
    @info "start training"
    for epoch in 1:hyperpar["number_of_epochs"]
    i = 0
        for batch in training_dataset
            i += 1
            x, phi, psi, acc, len = batch

            train_x_processed = fastaVocab(preprocess(x, hyperpar["max_seq_len"])) |> gpu
            train_y_processed = build_regressor_target_array(phi, psi, acc, hyperpar["max_seq_len"]) |> gpu
            train_mask_processed = prepare_regressor_mask(x, phi, psi, acc, hyperpar["max_seq_len"]) |> gpu

#            grad = gradient(()->train_loss(model, train_x_processed, train_y_processed, train_mask_processed), ps)

            l, back = Flux.pullback(ps) do
                train_loss(model, train_x_processed, train_y_processed, train_mask_processed)
            end
            grad = back(Flux.Zygote.sensitivity(l))

            update!(opt, ps, grad)

            if i%64 == 0

                Flux.testmode!(model)

                let
                batch_loss_ideal_train_set = 0.0
                batch_loss_ideal_test_set = 0.0

                for eval_train_batch in eval_train_dataset

                    eval_x, eval_phi, eval_psi, eval_acc, eval_len = eval_train_batch

                    eval_x_processed = fastaVocab(preprocess(x, hyperpar["max_seq_len"])) |> gpu
                    eval_y_processed = build_regressor_target_array(phi, psi, acc, hyperpar["max_seq_len"]) |> gpu
                    eval_mask_processed = prepare_regressor_mask(x, phi, psi, acc, hyperpar["max_seq_len"]) |> gpu

                    network_outputs = model(eval_x_processed)
                    batch_loss_ideal_train_set += eval_loss(network_outputs, eval_y_processed, eval_mask_processed)

                end #end IDEAL_TRAIN_EVAL loop

                for eval_test_batch in eval_test_dataset

                    eval_x, eval_phi, eval_psi, eval_acc, eval_len = eval_test_batch

                    eval_x_processed = fastaVocab(preprocess(x, hyperpar["max_seq_len"])) |> gpu
                    eval_y_processed = build_regressor_target_array(phi, psi, acc, hyperpar["max_seq_len"]) |> gpu
                    eval_mask_processed = prepare_regressor_mask(x, phi, psi, acc, hyperpar["max_seq_len"]) |> gpu

                    network_outputs = model(eval_x_processed)
                    batch_loss_ideal_test_set += eval_loss(network_outputs, eval_y_processed, eval_mask_processed)

                end #end IDEAL_TEST_EVAL loop

            ideal_train_loss = batch_loss_ideal_train_set
            ideal_test_loss = batch_loss_ideal_test_set

            Flux.trainmode!(model)

            checkpoint_model = cpu(model)
            @save "TintiNet-IGBT-Regressor-Checkpoint-$(Dates.format(Dates.now(), dformat)).bson" checkpoint_model opt
            @printf("Epoch: %d\t Batch: %d\t IdTrLos: %2.6f\t IdTsLoss: %2.6f\t Time: %s\n",
                    epoch, i, ideal_train_loss, ideal_test_loss, Dates.format(now(), "yyyy-mm-dd HH:MM:SS"))

            end # end let

            end # endif

        end # end for batch

    end # end for epoch

end # end function train!()

#####
# Uncomment next line to allow the script to undergo training.
# train!(hyperpar, training_dataset, model, ps, opt)


text:
module Propagation

export Trajectory, ODEProblem, propagate
export Manifold, EnsembleProblem
export initialstate, initialepoch, solution
export isperiodic, halo
export manifold, perturb
export stable_eigenvector, unstable_eigenvector

using ..States
using ..Calculations
using ..CoordinateFrames

import Dates: now
import LinearAlgebra: I
import Roots: find_zero

using Unitful
using Requires
using AstroTime
using Distributed
using StaticArrays
using LinearAlgebra
using DocStringExtensions
using AstrodynamicalModels
using DifferentialEquations

@template (FUNCTIONS, METHODS, MACROS) =
    """
    $(SIGNATURES)
    $(DOCSTRING)
    $(METHODLIST)
    """

@template (TYPES, CONSTANTS) =
    """
    $(TYPEDEF)
    $(DOCSTRING)
    """

include("Trajectories.jl")
include("Propagators.jl")
include("Halos.jl")
include("Manifolds.jl")

function __init__()
    @require DataFrames="a93c6f00-e57d-5684-b7b6-d8193f3e46c0" begin
        DataFrames.DataFrame(traj::Trajectory) = DataFrames.DataFrame(solution(traj))
    end
end

end # module


text:
tests = ["day" * lpad(day, 2, "0") for day = 1:25 if isfile("day" * lpad(day, 2, "0") * ".jl")]

for test in tests
    include(test * ".jl")
end


text:
using Documenter, Package

makedocs(modules = [Package], sitename = "Package.jl")

deploydocs(repo = "github.com/navidcy/Package.jl.git")


text:
@testset "Simulations" begin
    T = Float64
    @testset "Flat Interface" begin
        sys = Swalbe.SysConst(Lx=25, Ly=25, Tmax=200, tdump=100)
        h = Swalbe.run_flat(sys, "CPU", verbos=false)
        @test all(h .== 1.0)
        @test sum(h) == 25*25
        h = Swalbe.run_flat(sys, "CPU", verbos=true)
        # One dim version
        sys = Swalbe.SysConst_1D(L=25, Tmax=200, tdump=100)
        h = Swalbe.run_flat(sys, verbos=false)
        @test all(h .== 1.0)
        @test sum(h) == 25
        h = Swalbe.run_flat(sys, verbos=true)
    end
    @testset "Random Interface" begin
        sys = Swalbe.SysConst(Lx=25, Ly=25, Tmax=10000, tdump=5000)
        h = Swalbe.run_random(sys, "CPU", ϵ=0.1, verbos=false)
        difference = maximum(h) - minimum(h)
        @test difference < 0.02
        sys = Swalbe.SysConst(Lx=25, Ly=25, Tmax=100, tdump=500)
        h = Swalbe.run_random(sys, "CPU", ϵ=0.1, verbos=true)
        # One dim version
        sys = Swalbe.SysConst_1D(L=25, Tmax=10000, tdump=5000)
        h = Swalbe.run_random(sys, ϵ=0.1, verbos=false)
        difference = maximum(h) - minimum(h)
        @test difference < 0.02
        sys = Swalbe.SysConst_1D(L=25, Tmax=100, tdump=500)
        h = Swalbe.run_random(sys, ϵ=0.1, verbos=true)
    end
    @testset "Rayleigh Taylor" begin
        sys = Swalbe.SysConst(Lx=100, Ly=100, Tmax=1000, tdump=500, g=-0.002)
        _, diff = Swalbe.run_rayleightaylor(sys, "CPU", kx=4, ky=5, ϵ=0.01, verbos=true)
        @test diff[1] < diff[end]
        # TODO: Write a test to check the growthrate
        # One dim case
        sys = Swalbe.SysConst_1D(L = 512, Tmax=1000, tdump=500, g=-0.002)
        _, diff = Swalbe.run_rayleightaylor(sys, k=4, ϵ=0.01, verbos=true)
        @test diff[1] < diff[end]
    end
    @testset "Relaxing droplet" begin
        rads = 35
        sys = Swalbe.SysConst(Lx=150, Ly=150, Tmax=10000, δ=3.0)
        h, A = Swalbe.run_dropletrelax(sys, "CPU", radius=rads, verbos=false)
        # Initial droplet volume
        vol = π/3 * rads^3 * (2+cospi(1/6)) * (1-cospi(1/6))^2
        R1 = cbrt((rads^3*(2+cospi(1/6))*(1-cospi(1/6))^2)/((2+cospi(1/9))*(1-cospi(1/9))^2))
        r1 = sinpi(1/9)*R1
        # Identify the droplet, below this height is precursor layer
        drop1d = findall(h[75,:] .> 0.055)
        # Get the radius and height
        droprad = length(drop1d)/2
        droph = maximum(h)
        # Compute the droplets shape after 10k time steps
        vnum = 1/6*π*droph*(3*droprad^2 + droph^2)
        # Test that the droplet volume has not changed too much
        @test vol ≈ vnum atol = vol/100*10
        @test r1 ≈ droprad atol = r1/100*10
        @test A[1] < A[end]
        sys = Swalbe.SysConst(Lx=150, Ly=150, Tmax=100, δ=3.0, tdump=50)
        h, A = Swalbe.run_dropletrelax(sys, "CPU", radius=rads, verbos=true)
        # One dim case
        thresh = 0.05026
        sys = Swalbe.SysConst_1D(L = 150, Tmax=10000, δ=3.0)
        h, A = Swalbe.run_dropletrelax(sys, radius=rads, verbos=false)
        # Initial droplet volume
        vol = π/3 * rads^3 * (2+cospi(1/6)) * (1-cospi(1/6))^2
        R1 = cbrt((rads^3*(2+cospi(1/6))*(1-cospi(1/6))^2)/((2+cospi(1/9))*(1-cospi(1/9))^2))
        r1 = sinpi(1/9)*R1
        # Identify the droplet, below this height is precursor layer
        drop1d = findall(h .> 0.055)
        # Drop area (summation of height)
        drop_ar = zeros(sys.L)
        for i in 1:sys.L
            if h[i] > thresh
                drop_ar[i] = h[i]
            end
        end
        # Get the radius and height
        droprad = length(drop1d)/2
        droph = maximum(h)
        # Compute the droplets shape after 10k time steps
        vnum = 1/6*π*droph*(3*droprad^2 + droph^2)
        sph_r = droph/2 + length(findall(h .> thresh))^2/(8*droph)
        ar = sph_r^2/2*(π/9 - sinpi(1/9))
        # Test that the droplet volume has not changed too much
        # @test vol ≈ vnum atol = vol/100*10 # This breaks in 1d I have to investigate at some point in time
        @test ar ≈ sum(drop_ar) atol = ar/100*10
        @test r1 ≈ droprad atol = r1/100*10
        @test A[1] < A[end]
        sys = Swalbe.SysConst_1D(L = 200, Tmax=100, δ=3.0, tdump=50)
        h, A = Swalbe.run_dropletrelax(sys, radius=35, verbos=true)
    end
    @testset "Patterned droplet" begin
        sys = Swalbe.SysConst(Lx=150, Ly=150, Tmax=10000, δ=3.0)
        h = Swalbe.run_dropletpatterned(sys, "CPU", radius=35, θₛ = fill(1/9,sys.Lx, sys.Ly), verbos=false)
        # Initial droplet volume
        vol = π/3 * 35^3 * (2+cospi(1/6)) * (1-cospi(1/6))^2
        R1 = cbrt((35^3*(2+cospi(1/6))*(1-cospi(1/6))^2)/((2+cospi(1/9))*(1-cospi(1/9))^2))
        r1 = sinpi(1/9)*R1
        # Identify the droplet, below this height is precursor layer
        drop1d = findall(h[75,:] .> 0.055)
        # Get the radius and height
        droprad = length(drop1d)/2
        droph = maximum(h)
        # Compute the droplets shape after 10k time steps
        vnum = 1/6*π*droph*(3*droprad^2 + droph^2)
        # Test that the droplet volume has not changed too much
        @test vol ≈ vnum atol = vol/100*10
        @test r1 ≈ droprad atol = r1/100*10
        sys = Swalbe.SysConst(Lx=150, Ly=150, Tmax=100, tdump=50, δ=3.0)
        h = Swalbe.run_dropletpatterned(sys, "CPU", radius=35, θₛ = fill(1/9,sys.Lx, sys.Ly))
        # One dim case
        thresh = 0.05026
        sys = Swalbe.SysConst_1D(L=150, Tmax=10000, δ=3.0)
        h = Swalbe.run_dropletpatterned(sys, radius=35, θₛ = fill(1/9,sys.L), verbos=false)
        # Initial droplet volume
        vol = π/3 * 35^3 * (2+cospi(1/6)) * (1-cospi(1/6))^2
        R1 = cbrt((35^3*(2+cospi(1/6))*(1-cospi(1/6))^2)/((2+cospi(1/9))*(1-cospi(1/9))^2))
        r1 = sinpi(1/9)*R1
        # Identify the droplet, below this height is precursor layer
        drop1d = findall(h .> 0.055)
        # Drop area (summation of height)
        drop_ar = zeros(sys.L)
        for i in 1:sys.L
            if h[i] > thresh
                drop_ar[i] = h[i]
            end
        end
        # Get the radius and height
        droprad = length(drop1d)/2
        droph = maximum(h)
        # Compute the droplets shape after 10k time steps
        vnum = 1/6*π*droph*(3*droprad^2 + droph^2)
        sph_r = droph/2 + length(findall(h .> thresh))^2/(8*droph)
        ar = sph_r^2/2*(π/9 - sinpi(1/9))
        @test ar ≈ sum(drop_ar) atol = ar/100*10
        @test r1 ≈ droprad atol = r1/100*10
        sys = Swalbe.SysConst_1D(L=150, Tmax=100, tdump=50, δ=3.0)
        h = Swalbe.run_dropletpatterned(sys, radius=35, θₛ = fill(1/9,sys.L))
    end

    @testset "Sliding droplet" begin
        sys = Swalbe.SysConst(Lx=150, Ly=150, Tmax=5000, δ=2.0)
        # This will return the hydrodynamic moments: h, u, v
        mom = Swalbe.run_dropletforced(sys, "CPU", radius=35, fx=1e-4, verbos=false)
        # Now check that the center has moved
        @test findmax(mom[1])[2][1] ≠ 75
        @test findmax(mom[1])[2][2] == 75
        # Further test that the velocity is not too large
        @test all(mom[2] .< 0.1)
        @test all(mom[3] .< 0.1)
        sys = Swalbe.SysConst(Lx=150, Ly=150, Tmax=100, tdump=50, δ=3.0)
        h = Swalbe.run_dropletforced(sys, "CPU", radius=35, fx=1e-4)
        # One dim case
        sys = Swalbe.SysConst_1D(L=150, Tmax=5000, δ=2.0)
        # This will return the hydrodynamic moments: h, u, v
        mom = Swalbe.run_dropletforced(sys, radius=35, f=1e-4, verbos=false)
        # Now check that the center has moved
        @test findmax(mom[1])[2] ≠ Int(sys.L/2)
        # Further test that the velocity is not too large
        @test all(mom[2] .< 0.1)
        sys = Swalbe.SysConst_1D(L=150, Tmax=100, tdump=50, δ=3.0)
        h = Swalbe.run_dropletforced(sys, radius=35, f=1e-4)
    end
end


text:
using MyJuliaGEMM
using LinearAlgebra, Random, Test

Random.seed!(44)

for n in (48, 288, 576)
  A, B, C = rand(n, n), rand(n, n), rand(n, n)
  Cref = C + A * B
  mygemm!(C, A, B)
  @test Cref ≈ C
end

max_stars_repo_name: jkozdon/MyJuliaGEMM

max_issues_repo_name: jkozdon/MyJuliaGEMM

max_forks_repo_name: jkozdon/MyJuliaGEMM


text:
function init_observable_mappings!(m::RealBond)

    observables = OrderedDict{Symbol,Observable}()
    population_mnemonic = get(get_setting(m, :population_mnemonic))

    ############################################################################
    ## 1. GDP
    ############################################################################
    gdp_fwd_transform = function (levels)
        # FROM: Level of nominal GDP (FRED :GDP series)
        # TO:   The cyclical component of log real GDP

        levels[:rGDP] = nominal_to_real(:GDP, levels)
        levels[:lGDP] = log.(levels[:rGDP])
        _, levels[:lGDP_cyc_comp] = hpfilter(levels[:lGDP], 1600)

        levels[:lGDP_cyc_comp]
    end

    gdp_rev_transform = logleveltopct_annualized

    observables[:obs_gdp] = Observable(:obs_gdp, [:GDP__FRED, population_mnemonic, :GDPDEF__FRED],
                                       gdp_fwd_transform, gdp_rev_transform,
                                       "Log Real GDP", "")

    m.observable_mappings = observables
end


text:
using Oceananigans.Operators: ΔzC

#####
##### The basic idea: boundary condition types and data.
#####

"""
    BCType

Abstract supertype for boundary condition types.
"""
abstract type BCType end

"""
    Periodic

A type specifying a periodic boundary condition.

A condition may not be specified with a `Periodic` boundary condition.
"""
struct Periodic <: BCType end

"""
    Flux

A type specifying a boundary condition on the flux of a field.
"""
struct Flux <: BCType end

"""
    Gradient

A type specifying a boundary condition on the derivative or gradient of a field. Also
called a Neumann boundary condition.
"""
struct Gradient <: BCType end

"""
    Value

A type specifying a boundary condition on the value of a field. Also called a Dirchlet
boundary condition.
"""
struct Value <: BCType end

"""
    NoPenetration

A type specifying a no-penetration boundary condition for a velocity component that is normal to a wall.

Thus `NoPenetration` can only be applied to `u` along x, `v` along y, or `w` along z.
For all other cases --- fields located at (Cell, Cell, Cell), or `u`, `v`,
and `w` in (y, z), (x, z), and (x, y), respectively, either `Value`,
`Gradient`, or `Flux` conditions must be used.

A condition may not be specified with a `NoPenetration` boundary condition.

Note that this differs from a zero `Value` boundary condition as `Value` imposes values at the cell centers
(and could apply to tracers) while a no-penetration boundary condition only applies
to normal velocity components at a wall, where the velocity at the cell face collocated
at the wall is known and set to zero.
"""
struct NoPenetration <: BCType end

# Famous people

"""
    Dirchlet

An alias for the `Value` boundary condition type.
"""
const Dirchlet = Value

"""
    Neumann

An alias for the `Gradient` boundary condition type.
"""
const Neumann = Gradient

"""
    BoundaryCondition{C<:BCType}(condition)

Construct a boundary condition of type `C` with a `condition` that may be given by a
number, an array, or a function with signature:

    condition(i, j, grid, time, iteration, U, Φ, parameters) = # function definition

that returns a number and where `i` and `j` are indices along the boundary.

Boundary condition types include `Periodic`, `Flux`, `Value`, `Gradient`, and `NoPenetration`.
"""
struct BoundaryCondition{C<:BCType, T}
    condition :: T
end

# Some abbreviations to make life easier.
const BC   = BoundaryCondition
const FBC  = BoundaryCondition{<:Flux}
const PBC  = BoundaryCondition{<:Periodic}
const NPBC = BoundaryCondition{<:NoPenetration}
const VBC  = BoundaryCondition{<:Value}
const GBC  = BoundaryCondition{<:Gradient}
const NFBC = BoundaryCondition{Flux, Nothing}

BoundaryCondition(Tbc, c) = BoundaryCondition{Tbc, typeof(c)}(c)
bctype(bc::BoundaryCondition{C}) where C = C

# Adapt boundary condition struct to be GPU friendly and passable to GPU kernels.
Adapt.adapt_structure(to, b::BC{C, A}) where {C<:BCType, A<:AbstractArray} =
    BoundaryCondition(C, Adapt.adapt(to, parent(b.condition)))

PeriodicBC() = BoundaryCondition(Periodic, nothing)
NoPenetrationBC() = BoundaryCondition(NoPenetration, nothing)
NoFluxBC() = BoundaryCondition(Flux, nothing)

# Multiple dispatch on the type of boundary condition
getbc(bc::BC{C, <:Number}, args...)              where C = bc.condition
getbc(bc::BC{C, <:AbstractArray}, i, j, args...) where C = bc.condition[i, j]
getbc(bc::BC{C, <:Function}, args...)            where C = bc.condition(args...)

Base.getindex(bc::BC{C, <:AbstractArray}, inds...) where C = getindex(bc.condition, inds...)

Base.show(io::IO, bc::BC{C, T}) where {C, T} =
    println(io, "BoundaryCondition: type=$C, condition=$(bc.condition)")

#####
##### Wrapper for user-defined boundary condition functions
#####

"""
    BoundaryFunction{B, X1, X2, F}

A wrapper for user-defined boundary condition functions.
"""
struct BoundaryFunction{B, X1, X2, F} <: Function
    func :: F

    """
        BoundaryFunction{B, X1, X2}(func)

    A wrapper for user-defined boundary condition functions on the
    boundary specified by symbol `B` and at location `(X1, X2)`.

    Example
    =======
    julia> using Oceananigans: BoundaryCondition, BoundaryFunction, Flux, Cell

    julia> top_tracer_flux = BoundaryFunction{:z, Cell, Cell}((x, y, t) -> cos(2π*x) * cos(t))
    (::BoundaryFunction{:z,Cell,Cell,getfield(Main, Symbol("##7#8"))}) (generic function with 1 method)

    julia> top_tracer_bc = BoundaryCondition(Flux, top_tracer_flux);
    """
    function BoundaryFunction{B, X1, X2}(func) where {B, X1, X2}
        B ∈ (:x, :y, :z) || throw(ArgumentError("The boundary B at which the BoundaryFunction is
                                                to be applied must be either :x, :y, or :z."))
        new{B, X1, X2, typeof(func)}(func)
    end
end

@inline (bc::BoundaryFunction{:x, Y, Z})(j, k, grid, time, args...) where {Y, Z} =
    bc.func(ynode(Y, j, grid), znode(Z, k, grid), time)

@inline (bc::BoundaryFunction{:y, X, Z})(i, k, grid, time, args...) where {X, Z} =
    bc.func(xnode(X, i, grid), znode(Z, k, grid), time)

@inline (bc::BoundaryFunction{:z, X, Y})(i, j, grid, time, args...) where {X, Y} =
    bc.func(xnode(X, i, grid), ynode(Y, j, grid), time)

#####
##### Boundary conditions along particular coordinates
#####

"""
    CoordinateBoundaryConditions(left, right)

A set of two `BoundaryCondition`s to be applied along a coordinate x, y, or z.

The `left` boundary condition is applied on the negative or lower side of the coordinate
while the `right` boundary condition is applied on the positive or higher side.
"""
mutable struct CoordinateBoundaryConditions{L, R}
     left :: L
    right :: R
end

const CBC = CoordinateBoundaryConditions
PeriodicBCs() = CBC(PeriodicBC(), PeriodicBC())

# Here we overload setproperty! and getproperty to permit users to call
# the 'left' and 'right' bcs in the z-direction 'bottom' and 'top'
# and the 'left' and 'right' bcs in the y-direction 'south' and 'north'.
# Note that 'right' technically corresponds to face point N+1.
Base.setproperty!(cbc::CBC, side::Symbol, bc) = setbc!(cbc, Val(side), bc)
setbc!(cbc::CBC, ::Val{S}, bc) where S = setfield!(cbc, S, bc)
setbc!(cbc::CBC, ::Val{:bottom}, bc) = setfield!(cbc, :left, bc)
setbc!(cbc::CBC, ::Val{:top}, bc) = setfield!(cbc, :right, bc)
setbc!(cbc::CBC, ::Val{:south}, bc) = setfield!(cbc, :left, bc)
setbc!(cbc::CBC, ::Val{:north}, bc) = setfield!(cbc, :right, bc)

Base.getproperty(cbc::CBC, side::Symbol) = getbc(cbc, Val(side))
getbc(cbc::CBC, ::Val{S}) where S = getfield(cbc, S)
getbc(cbc::CBC, ::Val{:bottom}) = getfield(cbc, :left)
getbc(cbc::CBC, ::Val{:top}) = getfield(cbc, :right)
getbc(cbc::CBC, ::Val{:south}) = getfield(cbc, :left)
getbc(cbc::CBC, ::Val{:north}) = getfield(cbc, :right)

#####
##### Boundary conditions for Fields
#####

"""
    FieldBoundaryConditions

An alias for `NamedTuple{(:x, :y, :z)}` that represents a set of three `CoordinateBoundaryCondition`s
applied to a field along x, y, and z.
"""
const FieldBoundaryConditions = NamedTuple{(:x, :y, :z)}

show_field_boundary_conditions(bcs::FieldBoundaryConditions, padding="") =
    string("Oceananigans.FieldBoundaryConditions (NamedTuple{(:x, :y, :z)}), with boundary conditions", '\n',
           padding, "├── x: ", typeof(bcs.x), '\n',
           padding, "├── y: ", typeof(bcs.y), '\n',
           padding, "└── z: ", typeof(bcs.z))

Base.show(io::IO, fieldbcs::FieldBoundaryConditions) = print(io, show_field_boundary_conditions(fieldbcs))

"""
    FieldBoundaryConditions(x, y, z)

Construct a `FieldBoundaryConditions` using a `CoordinateBoundaryCondition` for each of the
`x`, `y`, and `z` coordinates.
"""
FieldBoundaryConditions(x, y, z) = FieldBoundaryConditions((x, y, z))

"""
    HorizontallyPeriodicBCs(;   top = BoundaryCondition(Flux, nothing),
                             bottom = BoundaryCondition(Flux, nothing))

Construct `FieldBoundaryConditions` with `Periodic` boundary conditions in the x and y
directions and specified `top` (+z) and `bottom` (-z) boundary conditions for u, v,
and tracer fields.

`HorizontallyPeriodicBCs` cannot be applied to the the vertical velocity w.
"""
function HorizontallyPeriodicBCs(;    top = BoundaryCondition(Flux, nothing),
                                   bottom = BoundaryCondition(Flux, nothing))

    x = PeriodicBCs()
    y = PeriodicBCs()
    z = CoordinateBoundaryConditions(bottom, top)

    return FieldBoundaryConditions(x, y, z)
end

"""
    ChannelBCs(; north = BoundaryCondition(Flux, nothing),
                 south = BoundaryCondition(Flux, nothing),
                   top = BoundaryCondition(Flux, nothing),
                bottom = BoundaryCondition(Flux, nothing))

Construct `FieldBoundaryConditions` with `Periodic` boundary conditions in the x
direction and specified `north` (+y), `south` (-y), `top` (+z) and `bottom` (-z)
boundary conditions for u, v, and tracer fields.

`ChannelBCs` cannot be applied to the the vertical velocity w.
"""
function ChannelBCs(;  north = BoundaryCondition(Flux, nothing),
                       south = BoundaryCondition(Flux, nothing),
                         top = BoundaryCondition(Flux, nothing),
                      bottom = BoundaryCondition(Flux, nothing)
                    )

    x = PeriodicBCs()
    y = CoordinateBoundaryConditions(south, north)
    z = CoordinateBoundaryConditions(bottom, top)

    return FieldBoundaryConditions(x, y, z)
end

#####
##### Boundary conditions for model solutions
#####

default_tracer_bcs(tracers, solution_bcs) = DefaultTracerBoundaryConditions(solution_bcs[1])

"""
    SolutionBoundaryConditions(tracers, proposal_bcs)

Construct a `NamedTuple` of `FieldBoundaryConditions` for a model with
fields `u`, `v`, `w`, and `tracers` from the proposal boundary conditions
`proposal_bcs`, which must contain the boundary conditions on `u`, `v`, and `w`
and may contain some or all of the boundary conditions on `tracers`.
"""
SolutionBoundaryConditions(tracers, proposal_bcs) =
    with_tracers(tracers, proposal_bcs, default_tracer_bcs, with_velocities=true)

"""
    HorizontallyPeriodicSolutionBCs(u=HorizontallyPeriodicBCs(), ...)

Construct `SolutionBoundaryConditions` for a horizontally-periodic model
configuration with solution fields `u`, `v`, `w`, `T`, and `S` specified by keyword arguments.

By default `HorizontallyPeriodicBCs` are applied to `u`, `v`, `T`, and `S`
and `HorizontallyPeriodicBCs(top=NoPenetrationBC(), bottom=NoPenetrationBC())` is applied to `w`.

Use `HorizontallyPeriodicBCs` when constructing non-default boundary conditions for `u`, `v`, `w`, `T`, `S`.
"""
function HorizontallyPeriodicSolutionBCs(;
    u = HorizontallyPeriodicBCs(),
    v = HorizontallyPeriodicBCs(),
    w = HorizontallyPeriodicBCs(top=NoPenetrationBC(), bottom=NoPenetrationBC()),
    tracerbcs...)

    return merge((u=u, v=v, w=w), tracerbcs)
end

"""
    ChannelSolutionBCs(u=ChannelBCs(), ...)

Construct `SolutionBoundaryConditions` for a reentrant channel model
configuration with solution fields `u`, `v`, `w`, `T`, and `S` specified by keyword arguments.

By default `ChannelBCs` are applied to `u`, `v`, `T`, and `S`
and `ChannelBCs(top=NoPenetrationBC(), bottom=NoPenetrationBC())` is applied to `w`.

Use `ChannelBCs` when constructing non-default boundary conditions for `u`, `v`, `w`, `T`, `S`.
"""
function ChannelSolutionBCs(;
    u = ChannelBCs(),
    v = ChannelBCs(north=NoPenetrationBC(), south=NoPenetrationBC()),
    w = ChannelBCs(top=NoPenetrationBC(), bottom=NoPenetrationBC()),
    tracerbcs...)

    return merge((u=u, v=v, w=w), tracerbcs)
end

# Default semantics
const BoundaryConditions = HorizontallyPeriodicSolutionBCs

#####
##### Tracer, tendency and pressure boundary condition "translators":
#####
#####   * Default boundary conditions on tracers are periodic or no flux and
#####     can be derived from boundary conditions on any field
#####
#####   * Boundary conditions on tendency terms are
#####     derived from the boundary conditions on their repsective fields.
#####
#####   * Boundary conditions on pressure are derived from boundary conditions
#####     on the north-south horizontal velocity, v.
#####

# Tracers
DefaultTracerBC(::BC) = BoundaryCondition(Flux, nothing)
DefaultTracerBC(::PBC) = PeriodicBC()

DefaultTracerCoordinateBCs(bcs) =
    CoordinateBoundaryConditions(DefaultTracerBC(bcs.left), DefaultTracerBC(bcs.right))

DefaultTracerBoundaryConditions(field_bcs) =
    FieldBoundaryConditions(Tuple(DefaultTracerCoordinateBCs(bcs) for bcs in field_bcs))

# Tendencies
TendencyBC(::BC) = BoundaryCondition(Flux, nothing)
TendencyBC(::PBC) = PeriodicBC()
TendencyBC(::NPBC) = NoPenetrationBC()

TendencyCoordinateBCs(bcs) =
    CoordinateBoundaryConditions(TendencyBC(bcs.left), TendencyBC(bcs.right))

TendencyFieldBCs(field_bcs) =
    FieldBoundaryConditions(Tuple(TendencyCoordinateBCs(bcs) for bcs in field_bcs))

TendenciesBoundaryConditions(solution_bcs) =
    NamedTuple{propertynames(solution_bcs)}(Tuple(TendencyFieldBCs(bcs) for bcs in solution_bcs))

# Pressure boundary conditions are either zero flux (Neumann) or Periodic.
# Note that a zero flux boundary condition is simpler than a zero gradient boundary condition.
PressureBC(::BC) = BoundaryCondition(Flux, nothing)
PressureBC(::PBC) = PeriodicBC()

function PressureBoundaryConditions(vbcs)
    x = CoordinateBoundaryConditions(PressureBC(vbcs.x.left), PressureBC(vbcs.x.right))
    y = CoordinateBoundaryConditions(PressureBC(vbcs.y.left), PressureBC(vbcs.y.right))
    z = CoordinateBoundaryConditions(PressureBC(vbcs.z.left), PressureBC(vbcs.z.right))
    return (x=x, y=y, z=z)
end

#####
##### Model boundary conditions for pressure, and solution, and tendency solutions
#####

const ModelBoundaryConditions = NamedTuple{(:solution, :tendency, :pressure)}

function ModelBoundaryConditions(tracers, proposal_bcs::NamedTuple)
    solution_boundary_conditions = SolutionBoundaryConditions(tracers, proposal_bcs)

    model_boundary_conditions = (solution = solution_boundary_conditions,
                                 tendency = TendenciesBoundaryConditions(solution_boundary_conditions),
                                 pressure = PressureBoundaryConditions(solution_boundary_conditions.v))
    return model_boundary_conditions
end

ModelBoundaryConditions(tracers, model_boundary_conditions::ModelBoundaryConditions) =
    model_boundary_conditions

function show_solution_boundary_conditions(bcs, padding)
    stringtuple = Tuple(string(
                  padding, "├── ", field, ": ",
                  show_field_boundary_conditions(getproperty(bcs, field), padding * "│   "), '\n')
                  for field in propertynames(bcs)[1:end-1])
    return string("Oceananigans.SolutionBoundaryConditions ",
                  "(NamedTuple{(:u, :v, :w, ...)}) with boundary conditions ", '\n', stringtuple...,
                  padding, "└── ", propertynames(bcs)[end], ": ",
                  show_field_boundary_conditions(bcs[end], padding * "    "))
end

Base.show(io::IO, bcs::ModelBoundaryConditions) =
    print(io,
          "Oceananigans.ModelBoundaryConditions (NamedTuple{(:solution, :tendency, :pressure)}) with ", '\n',
          "├── solution: ",
          show_solution_boundary_conditions(bcs.solution, "│   "), '\n',
          "├── tendency: ",
          show_solution_boundary_conditions(bcs.tendency, "│   "), '\n',
          "└── pressure: ",
          show_field_boundary_conditions(bcs.pressure, "    "))

#####
##### Algorithm for adding fluxes associated with non-trivial flux boundary conditions.
##### Inhomogeneous Value and Gradient boundary conditions are handled by filling halos.
#####

# Avoid some computation / memory accesses for Value, Gradient, Periodic, NoPenetration,
# and no-flux boundary conditions --- every boundary condition that does *not* prescribe
# a non-trivial flux.
const NotFluxBC = Union{VBC, GBC, PBC, NPBC, NFBC}
apply_z_bcs!(Gc, arch, grid, ::NotFluxBC, ::NotFluxBC, args...) = nothing
apply_y_bcs!(Gc, arch, grid, ::NotFluxBC, ::NotFluxBC, args...) = nothing

"""
    apply_z_bcs!(Gc, arch, grid, top_bc, bottom_bc, args...)

Apply flux boundary conditions to a field `c` by adding the associated flux divergence to
the source term `Gc` at the top and bottom.
"""
function apply_z_bcs!(Gc, arch, grid, bottom_bc, top_bc, args...)
    @launch device(arch) config=launch_config(grid, :xy) _apply_z_bcs!(Gc, grid, bottom_bc, top_bc, args...)
    return nothing
end

function apply_y_bcs!(Gc, arch, grid, left_bc, right_bc, args...)
    @launch device(arch) config=launch_config(grid, :xz) _apply_y_bcs!(Gc, grid, left_bc, right_bc, args...)
    return nothing
end

"""
    _apply_z_bcs!(Gc, grid, top_bc, bottom_bc, args...)

Apply a top and/or bottom boundary condition to variable `c`.
"""
function _apply_z_bcs!(Gc, grid, bottom_bc, top_bc, args...)
    @loop_xy i j grid begin
        apply_z_bottom_bc!(Gc, bottom_bc, i, j, grid, args...)
           apply_z_top_bc!(Gc, top_bc,    i, j, grid, args...)
    end
end

function _apply_y_bcs!(Gc, grid, left_bc, right_bc, args...)
    @loop_xz i k grid begin
         apply_y_left_bc!(Gc, left_bc,  i, k, grid, args...)
        apply_y_right_bc!(Gc, right_bc, i, k, grid, args...)
    end
end

# Fall back functions for boundary conditions that are not of type Flux.
@inline apply_z_top_bc!(args...) = nothing
@inline apply_z_bottom_bc!(args...) = nothing

@inline apply_y_right_bc!(args...) = nothing
@inline apply_y_left_bc!(args...) = nothing

# Shortcuts for 'no-flux' boundary conditions.
@inline apply_z_top_bc!(Gc, ::NFBC, args...) = nothing
@inline apply_z_bottom_bc!(Gc, ::NFBC, args...) = nothing

@inline apply_y_right_bc!(Gc, ::NFBC, args...) = nothing
@inline apply_y_left_bc!(Gc, ::NFBC, args...) = nothing

"""
    apply_z_top_bc!(Gc, top_flux::BC{<:Flux}, i, j, grid, args...)

Add the part of flux divergence associated with a top boundary condition on `c`.
Note that because

    `tendency = ∂c/∂t = Gc = - ∇ ⋅ flux`

a positive top flux is associated with a *decrease* in `Gc` near the top boundary.
If `top_bc.condition` is a function, the function must have the signature

    `top_bc.condition(i, j, grid, boundary_condition_args...)`
"""
@inline apply_z_top_bc!(Gc, top_flux::BC{<:Flux}, i, j, grid, args...) =
    @inbounds Gc[i, j, grid.Nz] -= getbc(top_flux, i, j, grid, args...) / ΔzF(i, j, grid.Nz, grid)

@inline apply_y_right_bc!(Gc, right_flux::BC{<:Flux}, i, k, grid, args...) =
    @inbounds Gc[i, grid.Ny, k] -= getbc(right_flux, i, k, grid, args...) / Δy(i, grid.Ny, k, grid)

"""
    apply_z_bottom_bc!(Gc, bottom_flux::BC{<:Flux}, i, j, grid, args...)

Add the flux divergence associated with a bottom flux boundary condition on `c`.
Note that because

    `tendency = ∂c/∂t = Gc = - ∇ ⋅ flux`

a positive bottom flux is associated with an *increase* in `Gc` near the bottom boundary.
If `bottom_bc.condition` is a function, the function must have the signature

    `bottom_bc.condition(i, j, grid, boundary_condition_args...)`
"""
@inline apply_z_bottom_bc!(Gc, bottom_flux::BC{<:Flux}, i, j, grid, args...) =
    @inbounds Gc[i, j, 1] += getbc(bottom_flux, i, j, grid, args...) / ΔzF(i, j, 1, grid)

@inline apply_z_left_bc!(Gc, bottom_flux::BC{<:Flux}, i, k, grid, args...) =
    @inbounds Gc[i, 1, k] += getbc(left_flux, i, k, grid, args...) / Δy(i, 1, k, grid)


text:
"""
    smd_combined, smd, smd_preclustered = frameconnect(smd::SMD, params::ParamStruct)

Connect repeated localizations of the same emitter in `smd`.

# Description
Repeated localizations of the same emitter present in `smd` are connected and
combined into higher precision localizations of that emitter.  This is done by
1) forming pre-clusters of localizations, 2) estimating rate parameters from
the pre-clusters, 3) solving a linear assignment problem for connecting
localizations in each pre-cluster, and 4) combining the connected localizations
using their MLE position estimate assuming Gaussian noise.
"""
function frameconnect(smd::SMD, params::ParamStruct = ParamStruct())
    # Generate pre-clusters of localizations in `smd`.
    smd_preclustered = precluster(smd, params)
    clusterdata = organizeclusters(smd_preclustered)

    # Estimate rate parameters.
    params.k_on, params.k_off, params.k_bleach, params.p_miss =
        estimateparams(smd_preclustered, clusterdata)

    # Estimate the underlying density of emitters.
    params.initialdensity =
        estimatedensities(smd_preclustered, clusterdata, params)

    # Connect localizations in `smd` by solving the LAP.
    nframes = isempty(smd.nframes) ? maximum(smd.framenum) : smd.nframes
    smd.connectID = connectlocalizations(
        smd_preclustered.connectID, clusterdata, params, nframes)

    # Combine the connected localizations into higher precision localizations.
    smd_combined = combinelocalizations(smd)

    return smd_combined, smd, smd_preclustered, params
end


text:
const lines = readlines()

parse_dot(l) = CartesianIndex((1 .+ parse.(Int, split(l, ",")))...)
function parse_fold(l)
    axis, v = split(l, "=")
    axis = Symbol(axis[end])
    v = parse(Int, v) + 1
    axis, v
end

const dots = parse_dot.(lines[begin:793])
const folds = parse_fold.(lines[795:end])

const paper = begin
    p = falses(maximum(dots).I)
    p[dots] .= true
    p
end

function fold(m, ::Val{:x}, pos)
    xsize = size(m, 1)
    n = copy(m[begin:pos-1, :])
    for I in CartesianIndices(n)
        xx, yy = I.I
        if (2pos - xx) ≤ xsize
            n[I] |= m[2pos-xx, yy]
        end
    end
    n
end

function fold(m, ::Val{:y}, pos)
    ysize = size(m, 2)
    n = copy(m[:, begin:pos-1])
    for I in CartesianIndices(n)
        xx, yy = I.I
        if (2pos - yy) ≤ ysize
            n[I] |= m[xx, 2pos-yy]
        end
    end
    n
end

fold(m, s::Symbol, pos) = fold(m, Val{s}(), pos)

# Part One
fold(paper, folds[1]...) |> count |> println

# Part Two
let pp = paper
    for f in folds
        pp = fold(pp, f...)
    end
    for r in eachcol(pp)
        println(join([rr ? "#" : " " for rr in r]))
    end
end

max_stars_repo_name: pLOPeGG/JuliaOfCode

max_issues_repo_name: pLOPeGG/JuliaOfCode

max_forks_repo_name: pLOPeGG/JuliaOfCode


text:
# * Basis construction

struct FEDVR{T,R<:Real,O<:AbstractVector} <: Basis{T}
    t::Vector{R}
    order::O
    i₀::Int
    t₀::R
    eiϕ::T
    x::Vector{T}
    wⁱ::Vector{Vector{R}}
    n::Vector{T}
    elems::Vector{UnitRange{Int}}
    function FEDVR(t::AbstractVector{R}, order::O; t₀::R=zero(R), ϕ::R=zero(R)) where {R<:Real,O<:AbstractVector}
        @assert length(order) == length(t)-1
        @assert all(order .> 1)

        i₀,eiϕ,T = if ϕ ≠ zero(R)
            findfirst(tt -> tt ≥ t₀, t),exp(im*ϕ),complex(R)
        else
            1,one(R),R
        end
        i₀ === nothing && throw(ArgumentError("Complex scaling starting point outside grid $(t)"))
        t₀ = t[i₀]

        x = Vector{Vector{T}}()
        w = Vector{Vector{R}}()
        for i in eachindex(order)
            xw = element_grid(order[i], t[i], t[i+1], t₀, i ≥ i₀ ? eiϕ : one(T))
            push!(x, xw[1])
            push!(w, xw[2])
        end

        rot = (i,v) -> i ≥ i₀ ? eiϕ*v : v
        n = [one(T) ./ .√(rot(i,wⁱ)) for (i,wⁱ) in enumerate(w)]
        for i in 1:length(order)-1
            n[i][end] = n[i+1][1] = one(T) ./ √(rot(i,w[i][end]) + rot(i+1,w[i+1][1]))
        end

        X = vcat(x[1], [x[i][2:end] for i in 2:length(order)]...)
        N = vcat(n[1], [n[i][2:end] for i in 2:length(order)]...)

        elems = [1:order[1]]

        l = order[1]
        for i = 2:length(order)
            l′ = l+order[i]-1
            push!(elems, l:l′)
            l = l′
        end

        new{T,R,O}(t,order,i₀,t₀,eiϕ,X,[w...],N,elems)
    end
end

FEDVR(t::AbstractVector{T},order::Integer; kwargs...) where T =
    FEDVR(t, Fill(order,length(t)-1); kwargs...)

const RestrictedFEDVR{T} = RestrictedQuasiArray{T,2,<:FEDVR{T}}
const AdjointRestrictedFEDVR{T} = AdjointRestrictedQuasiArray{T,2,<:FEDVR{T}}

const FEDVROrRestricted{T} = BasisOrRestricted{<:FEDVR{T}}
const AdjointFEDVROrRestricted{T} = AdjointBasisOrRestricted{<:FEDVR{T}}

# * Properties

axes(B::FEDVR) = (Inclusion(first(B.t)..last(B.t)), Base.OneTo(length(B.x)))
size(B::FEDVR) = (ℵ₁, length(B.x))
==(A::FEDVR,B::FEDVR) = A.t == B.t && A.order == B.order
Base.hash(B::FEDVR, h::UInt) = hash(B.t, hash(B.order, h))

assert_compatible_bases(A::FEDVROrRestricted, B::FEDVROrRestricted) =
    parent(A) == parent(B) ||
        throw(ArgumentError("Can only multiply FEDVRs of the same order and sharing the same nodes"))

distribution(B::FEDVR) = distribution(B.t)

order(B::FEDVR) = B.order
order(B::RestrictedFEDVR) = order(parent(B))[elements(B)]
order(B::FEDVROrRestricted, i) = order(parent(B))[i]

nel(B::FEDVR) = length(order(B))
element_boundaries(B::FEDVR) = vcat(1,last.(B.elems))

elements(B::FEDVR) = 1:nel(B)
function elements(B::RestrictedFEDVR)
    is=indices(B,2)
    a,b = first(is),last(is)
    elems = parent(B).elems
    i = findfirst(e -> a ∈ e, elems)
    j = findlast(e -> b ∈ e, elems)
    (isnothing(i) || isnothing(j)) && throw(ArgumentError("Invalid restriction $(a)..$(b)"))
    i:j
end

complex_rotate(x,B::FEDVR{T}) where {T<:Real} = x
complex_rotate(x,B::FEDVR{T}) where {T<:Complex} = x < B.t₀ ? x : B.t₀ + (x-B.t₀)*B.eiϕ

macro elem(B,v,i)
    :(@view($(esc(B)).$v[$(esc(B)).elems[$(esc(i))]]))
end

function show(io::IO, B::FEDVR{T}) where T
    write(io, "FEDVR{$(T)} basis with $(nel(B)) elements on $(axes(B,1).domain)")
    if T <: Complex
        ics = B.t₀ <= first(B.t)
        rot = printfmt(io, " with {1:s} @ {2:0.2f}°", ics ? "ICS" : "ECS", rad2deg(angle(B.eiϕ)))
        !ics && printfmt(io, " starting at {1:0.2f}", B.t₀)
    end
end

function show(io::IO, B::RestrictedFEDVR{T}) where T
    B′ = parent(B)
    is = indices(B,2)
    a,b = first(is),last(is)
    els = elements(B)
    N = length(B′.x)
    show(io, B′)
    write(io, ", restricted to elements $(els), basis functions $(a)..$(b) $(a>1 || b>N ? "⊂" : "⊆") 1..$(N)")
end

function block_orders(B::FEDVROrRestricted)
    is = indices(B,2)
    elb = element_boundaries(parent(B))
    els = elements(B)
    a,b = first(is),last(is)

    # If the restriction only includes the bridge function of the
    # first/last element covered, compute the block structure for one
    # element further in, since it simplifies.
    first(els) ≠ nel(parent(B)) && a == elb[first(els)+1] && (els = els[2:end])
    last(els) > first(els) && b == elb[last(els)] && (els = els[1:end-1])
    # Find out how much we need to "shave off" from the first and last
    # elements (inside the restriction), respectively.
    s = a-elb[first(els)]
    e = elb[last(els)+1]-b

    o = order(parent(B))[els]
    o,s,e
end

function block_structure(B::FEDVROrRestricted)
    o,s,e = block_orders(B)
    if length(o) > 1
        # First compute interior block structure; for each element,
        # the interior block and the block for the bridge function
        # connecting to /next/ element is added.
        b = o -> o > 2 ? [o-2,1] : [1]
        bs = reduce(vcat, map(b, o[2:end-1]), init=Int[])
        # Add first element.
        prepend!(bs, [o[1]-1-s,1])
        # Add last element.
        push!(bs, o[end]-1-e)
        bs
    else
        [first(o)-(s+e)]
    end
end

function block_bandwidths(B::FEDVROrRestricted, rows::Vector{<:Integer}=block_structure(B))
    nrows = length(rows)
    if nrows > 1
        bw = o -> o > 2 ? [2,1] : [1]
        # This actually generates one element "too much", but it does
        # not matter, since for the second last block, a lower
        # bandwidth of two extends outside the matrix, and that is
        # discarded.
        bws = bw.(first(block_orders(B)))
        l = vcat(1,bws[2:end]...)
        u = vcat(reverse.(bws[1:end-1])...,1)
        length(l) < nrows && (l = vcat(l,0))
        length(u) < nrows && (u = vcat(0,u))
        l,u
    else
        [0],[0]
    end
end

locs(B::FEDVR) = B.x

real_locs(B::FEDVR{<:Real}) = locs(B)

function real_locs(B::FEDVR{T}) where T
    R = real(T)
    nx = length(B.x)
    x = Vector{R}(undef, nx)
    ii = 1
    for i in eachindex(B.order)
        xw = element_grid(B.order[i], B.t[i], B.t[i+1], B.t[1], one(R))
        copyto!(view(x, ii:nx), xw[1])
        ii += length(xw[1])-1
    end
    x
end

IntervalSets.leftendpoint(B::FEDVR) = B.x[1]
IntervalSets.rightendpoint(B::FEDVR) = B.x[end]

weight(B::FEDVR, j) = inv(B.n[j])
weights(B::FEDVR) = inv.(B.n)

inverse_weight(B::FEDVR, j) = B.n[j]
inverse_weights(B::FEDVR) = B.n

vandermonde(B::FEDVR) = BandedMatrix(0 => B.n)
metric(B::FEDVROrRestricted{T}) where T = Diagonal(Ones{T}(size(B,2)))
metric_shape(B::FEDVROrRestricted) = Diagonal

# * Basis functions

# Evaluate basis function m of element i on a single coordinate x
function getindex(B::FEDVR{T}, x::Real, i::Integer, m::Integer) where T
    (x < B.t[i] || x > B.t[i+1]) && return zero(T)
    xⁱ = @elem(B,x,i)
    χ = @elem(B,n,i)[m]
    x′ = complex_rotate(x, B)
    # Replace with Fornberg interpolation?
    for j = 1:B.order[i]
        j == m && continue
        χ *= (x′ - xⁱ[j])/(xⁱ[m]-xⁱ[j])
    end
    χ
end

checkbounds(B::FEDVR{T}, x::Real, k::Integer) where T =
    x ∈ axes(B,1) && 1 ≤ k ≤ size(B,2) || throw(BoundsError())

"""
   findelement(B::FEDVR, k[, i=1, m=k])

Find the finite-element of `B` that contains the `k`th basis function,
optionally starting the search from element `i` and basis function `m`
of that element.
"""
function findelement(B::FEDVR, k, i=1, m=k)
    k ∈ axes(B,2) || throw(BoundsError(B, k))
    while i < length(B.order) && m > B.order[i]
        m -= B.order[i] - 1
        i += 1
    end
    i,m
end

# Evaluate a single basis function on a single coordinate
@inline function getindex(B::FEDVR{T}, x::Real, k::Integer) where T
    # @boundscheck checkbounds(B, x, k) # Slow
    i,m = findelement(B, k)
    x < B.t[i] && return zero(T)
    if x > B.t[i+1]
        # The last basis function of an element is also the first of
        # the next element (the bridge function); if we are past the
        # element boundary, move to next element.
        if i < length(B.t)-1 && m == B.order[i]
            i += 1
            m = 1
        else
            return zero(T)
        end
    end
    B[x,i,m]
end

# Evaluate a specific basis function in a specific element on all
# coordinates in x.
function basis_function!(χ, B::FEDVR{T}, x::AbstractVector, i, m) where T
    for j ∈ findall(in(B.t[i]..B.t[i+1]), x)
        χ[j] = B[x[j],i,m]
    end
    χ
end

# Evaluate a selection of basis functions on all coordinates in x
function getindex(B::FEDVR{T}, x::AbstractVector, sel::AbstractVector) where T
    nj = length(sel)
    χ = spzeros(T, length(x), nj)

    j = 1
    ni = length(B.t)-1
    i,m = findelement(B, sel[1])
    imax,mmax = findelement(B, sel[end])

    while j ≤ nj
        basis_function!(view(χ, :, j), B, x, i, m)

        m += 1
        if m > B.order[i] && i < ni
            i += 1
            m = 1
        else
            j += 1
        end
    end

    χ
end

# Evaluate a single basis function on all the coordinates in x
function getindex(B::FEDVROrRestricted{T}, x::AbstractVector, j::Integer) where T
    χ = spzeros(T, length(x))
    i,m = findelement(parent(B), j+first(indices(B,2))-1)
    basis_function!(χ, parent(B), x, i, m)

    # If basis function j is a bridge function between elements i and
    # i+1, also evaluate its contribution in the next element
    pB = parent(B)
    i < length(pB.t)-1 && m == pB.order[i] &&
        basis_function!(χ, parent(B), x, i+1, 1)

    χ
end

getindex(B::FEDVROrRestricted, x::AbstractVector, ::Colon) =
    getindex(parent(B), x, indices(B,2))

getindex(B::RestrictedFEDVR, x::AbstractVector, sel::AbstractVector) =
    getindex(parent(B), x, indices(B,2)[sel])

# * Types

const FEDVRArray{T,N,B<:FEDVROrRestricted} = FuncArray{T,N,B}
const FEDVRVector{T,B<:FEDVROrRestricted} = FuncVector{T,B}
const FEDVRMatrix{T,B<:FEDVROrRestricted} = FuncMatrix{T,B}
const FEDVRVecOrMat{T,B<:FEDVROrRestricted} = FuncVecOrMat{T,B}

const AdjointFEDVRArray{T,N,B<:FEDVROrRestricted} = AdjointFuncArray{T,N,B}
const AdjointFEDVRVector{T,B<:FEDVROrRestricted} = AdjointFuncVector{T,B}
const AdjointFEDVRMatrix{T,B<:FEDVROrRestricted} = AdjointFuncMatrix{T,B}
const AdjointFEDVRVecOrMat{T,B<:FEDVROrRestricted} = AdjointFuncVecOrMat{T,B}

# * Mass matrix

@simplify function *(Ac::QuasiAdjoint{<:Any,<:FEDVROrRestricted}, B::FEDVROrRestricted)
    one(eltype(B))*combined_restriction(parent(Ac),B)
end

# * Basis inverses

@simplify function *(A⁻¹::PInvQuasiMatrix{<:Any,<:Tuple{<:FEDVR}},
                     B::FEDVR)
    A = parent(A⁻¹)
    A == B || throw(ArgumentError("Cannot multiply basis with inverse of other basis"))
    I
end

@simplify function *(A::FEDVR,
                     B⁻¹::PInvQuasiMatrix{<:Any,<:Tuple{<:FEDVR}})
    B = parent(B⁻¹)
    A == B || throw(ArgumentError("Cannot multiply basis with inverse of other basis"))
    I
end

@simplify function *(A⁻¹::PInvQuasiMatrix{<:Any,<:Tuple{<:FEDVR}},
                     v::FEDVRArray)
    A = parent(A⁻¹)
    B,c = v.args
    A == B || throw(ArgumentError("Cannot multiply basis with inverse of other basis"))
    c
end

@simplify function *(v::AdjointFEDVRArray,
                     B⁻¹::PInvQuasiMatrix{<:Any,<:Tuple{<:FEDVR}})
    c,Ac = v.args
    B = parent(B⁻¹)
    parent(Ac) == B || throw(ArgumentError("Cannot multiply basis with inverse of other basis"))
    c
end

# * Function interpolation

function Base.:(\ )(B::FEDVR, f::BroadcastQuasiArray)
    axes(f,1) == axes(B,1) ||
        throw(DimensionMismatch("Function on $(axes(f,1).domain) cannot be interpolated over basis on $(axes(B,1).domain)"))

    T = promote_type(eltype(B), eltype(f))
    v = zeros(T, size(B,2))
    for i ∈ 1:nel(B)
        @. v[B.elems[i]] += B.wⁱ[i]*f[@elem(B,x,i)]
    end
    v .*= B.n
    v
end

function Base.:(\ )(B::RestrictedFEDVR, f::BroadcastQuasiArray)
    axes(f,1) == axes(B,1) ||
        throw(DimensionMismatch("Function on $(axes(f,1).domain) cannot be interpolated over basis on $(axes(B,1).domain)"))

    B′ = parent(B)
    is = indices(B,2)
    a,b = first(is),last(is)

    n = size(B′,2)
    T = promote_type(eltype(B), eltype(f))
    v = zeros(T, size(B,2))
    for i ∈ elements(B)
        sel = B′.elems[i]
        # Find which basis functions of finite-element `i` should be
        # evaluated.
        subsel = if a<sel[1] && b > sel[end]
            # Element is completely within the restricted basis.
            Colon()
        else
            # Element straddles restriction; find subset of functions
            # that are within the restriction.
            s = min(max(a,sel[1]),sel[end])
            e = max(min(b,sel[end]),sel[1])
            findfirst(isequal(s),sel):findfirst(isequal(e),sel)
        end

        @. v[sel[subsel] .- (a-1)] += @view((B′.wⁱ[i]*f[@elem(B′,x,i)])[subsel])
    end
    v .*= @view(B′.n[a:b])
    v
end


text:
# This file is part of Kpax3. License is MIT.

abstract type StateList end

mutable struct AminoAcidStateList <: StateList
  state::Vector{AminoAcidState}
  logpp::Vector{Float64}
  rank::Vector{Int}
end

function AminoAcidStateList(data::Matrix{UInt8},
                            partition::Vector{Int},
                            priorR::PriorRowPartition,
                            priorC::AminoAcidPriorCol,
                            settings::KSettings)
  state = Array{AminoAcidState}(undef, settings.popsize)
  logpp = zeros(Float64, settings.popsize)

  state[1] = AminoAcidState(data, partition, priorR, priorC, settings)
  logpp[1] = state[1].logpp

  if settings.verbose
    Printf.@printf("Creating initial population... ")
  end

  R = zeros(Int, size(data, 2))
  for i in 2:settings.popsize
    copyto!(R, partition)
    modifypartition!(R, state[1].k)
    state[i] = AminoAcidState(data, R, priorR, priorC, settings)
    logpp[i] = state[i].logpp
  end

  rank = sortperm(logpp, rev=true)

  if settings.verbose
    Printf.@printf("done\n")
  end

  AminoAcidStateList(state, logpp, rank)
end

function AminoAcidStateList(popsize::Int,
                            s::AminoAcidState)
  state = Array{AminoAcidState}(undef, popsize)
  logpp = zeros(Float64, popsize)
  rank = Int[i for i in 1:popsize]

  for i in 1:popsize
    state[i] = copystate(s)
    logpp[i] = state[i].logpp
  end

  AminoAcidStateList(state, logpp, rank)
end

function copystatelist!(dest::AminoAcidStateList,
                        src::AminoAcidStateList,
                        popsize::Int)
  for i in 1:popsize
    copystate!(dest.state[i], src.state[src.rank[i]])
    dest.logpp[i] = src.logpp[src.rank[i]]
    dest.rank[i] = i
  end

  nothing
end


text:
export PolygonalPrism

const PolygonalPrism{K,K2} = Prism{Polygon{K,K2}}

# Below, if we called PolygonalPrism(c, ...) in the function body, it would call the inner
# constructor Prism{Polygon{K,K2,Nothing}}(c, ...) because PolygonalPrism = Prism{Polygon{K,K2,Nothing}},
# which is not what we want.
# To call the outer constructor of Prism, we should call Prism(c, ...) instead of PolygonalPrism(c, ...).
PolygonalPrism(c::SVector{3,<:Real},
               v::SMatrix{2,K,<:Real},  # 2D coordinates of base vertices in projected prism coordinates
               h::Real=Inf,
               a::SVector{3,<:Real}=SVector(0.0,0.0,1.0)
               ) where {K} =
    (â = normalize(a); Prism(c, Polygon(v), h, [orthoaxes(â)... â]))

PolygonalPrism(c::AbstractVector{<:Real},  # center of prism
               v::AbstractMatrix{<:Real},  # vertices of base polygon
               h::Real=Inf,  # height of prism
               a::AbstractVector{<:Real}=[0.0,0.0,1.0]) =  # axis direction of prism
    (K = size(v,1); PolygonalPrism(SVector{3}(c), SMatrix{2,K}(v), h, SVector{3}(a)))

# Return the bounds of the center cut with respect to the prism center.
function bounds_ctrcut(s::PolygonalPrism{K}) where {K}
    p = inv(s.p)  # projection matrix to prism coordinates: rows are not only unit vectors, but also orthogonal
    v = [s.b.v; @SMatrix(zeros(1,K))]  # SMatrix{3,K}: 3D vectices in prism axis coordinates
    w = p * v  # SMatrix{3,K}: vertices in external coordinates

    return minimum(w, dims=Val(2))[:,1], maximum(w, dims=Val(2))[:,1]
end


text:
"""
    natural_from_canonical(y::AbstractVector{<:Real}, σ²::AbstractVector{<:Real})

Compute natural parameters `η1` and `η2` from canonical paramters `y` and `σ²`.
Inverse of [`canonical_from_natural`](@ref).
"""
function natural_from_canonical(y::AbstractVector{<:Real}, σ²::AbstractVector{<:Real})
    η1, η2 = natural_from_canonical(y, Diagonal(σ²))
    return η1, η2.diag
end

function natural_from_canonical(y::AbstractVector{<:Real}, Σ::Diagonal{<:Real})
    return Σ \ y, -inv(Σ) / 2
end

"""
    canonical_from_natural(η1::AbstractVector{<:Real}, η2::AbstractVector{<:Real})

Compute canonical parameters `y` and `σ²` from natural parameters `η1` and `η2`.
Inverse of [`natural_from_canonical`](@ref).
"""
function canonical_from_natural(η1::AbstractVector{<:Real}, η2::AbstractVector{<:Real})
    y, Σ = canonical_from_natural(η1, Diagonal(η2))
    return y, diag(Σ)
end

function canonical_from_natural(η1::AbstractVector{<:Real}, η2::Diagonal{<:Real})
    Σ = -inv(η2) / 2
    y = Σ * η1
    return y, Σ
end

"""
    expectation_from_canonical(m::AbstractVector{<:Real}, σ²::AbstractVector{<:Real})

Compute expectation parameters `m1` and `m2` from canonical parameters `m` and `σ²`.
Inverse of [`expectation_from_canonical`](@ref).
"""
function expectation_from_canonical(m::AbstractVector{<:Real}, σ²::AbstractVector{<:Real})
    m, S = expectation_from_canonical(m, Diagonal(σ²))
    return m, diag(S)
end

function expectation_from_canonical(m::AbstractVector{<:Real}, Σ::Diagonal{<:Real})
    return m, Diagonal(m.^2 + diag(Σ))
end

"""
    canonical_from_expectation(m1::AbstractVector{<:Real}, m2::AbstractVector{<:Real})

Compute canonical parameters `m` and `σ²` from expectation parameters `m1` and `m2`.
Inverse of [`expectation_from_canonical`](@ref).
"""
function canonical_from_expectation(m1::AbstractVector{<:Real}, m2::AbstractVector{<:Real})
    m1, M2 = canonical_from_expectation(m1, Diagonal(m2))
    return m1, diag(M2)
end

function canonical_from_expectation(m1::AbstractVector{<:Real}, M2::Diagonal{<:Real})
    return m1, Diagonal(diag(M2) .- m1.^2)
end


text:
using MTKHelpers
import MTKHelpers as CP
using Test

using DifferentialEquations, ModelingToolkit
using StaticArrays, LabelledArrays, NamedArrays
using ForwardDiff
using NamedArrays
# #push!(LOAD_PATH, expanduser("~/julia/turingtools/")) # access local package repo
# import AxisArrays: AxisArray
using ComponentArrays
import ComponentArrays as CA

#include("test/samplesystem.jl")
include("samplesystem.jl")

@testset "util_componentarrays" begin
    #include("test/test_util_componentarrays.jl")
    include("test_util_componentarrays.jl")
end;

@testset "ProblemParSetter_dummy" begin
    #include("test/test_problemparsetter_dummy.jl")
    include("test_problemparsetter_dummy.jl")
end;

@testset "ProblemParSetter_sym" begin
    #include("test/test_problemparsetter_sym.jl")
    include("test_problemparsetter_sym.jl")
end;

@testset "ProblemParSetter" begin
    #include("test/test_problemparsetter.jl")
    include("test_problemparsetter.jl")
end;

@testset "utilities" begin
    include("test_util.jl")
end;

@testset "smoothstep" begin
    include("test_smoothstep.jl")
end;

@testset "cairomakie" begin
    include("test_cairomakie.jl")
end;

@testset "solution" begin
    include("test_solution.jl")
end;


text:
@testset "PlotsData_tests.jl" begin
    when = [1,2,3]
    times = [1,3]
    marks = [2, 3]

    ts = SEMPP.TimeSeries(times)
    mts = SEMPP.MarkedTimeSeries(times, marks)

    μ = 2.0
    ϕ = 1
    γ = 1
    δ = 1.0

    GPD = Distributions.GeneralizedPareto
    EGPD1 = EGPD.EGPpower

    ξ = 1.5
    β = 2
    α = 1.1
    κ = 1.5

    sepp = SEMPP.SEPPExpKern(μ, ϕ, γ)
    sempp_egpd = SEMPPExpKern(μ, ϕ, γ, δ, EGPD1, ξ, α, β, κ)
    sempp_gpd = SEMPPExpKern(μ, ϕ, γ, δ, GPD, ξ, α, β)

    @testset "pp_analysis" begin
        s, p = SEMPP.pp_analysis(sepp, ts)
        @test length(s) == 2
        @test length(p) == 2

        s, p = SEMPP.pp_analysis(sempp_gpd, mts)
        @test length(s) == 2
        @test length(p) == 2
    end

    @testset "transformed_marks_ecdf" begin
        s, p = SEMPP.transformed_marks_ecdf(sempp_egpd, mts)
        @test length(s) == 2
        @test length(p) == 2

        s, p = SEMPP.transformed_marks_ecdf(sempp_gpd, mts)
        @test length(s) == 2
        @test length(p) == 2
    end

    @testset "marks_unit_exponential_qq" begin
        s, p = SEMPP.marks_unit_exponential_qq(sempp_egpd, mts)
        @test length(s) == 2
        @test length(p) == 2

        s, p = SEMPP.marks_unit_exponential_qq(sempp_gpd, mts)
        @test length(s) == 2
        @test length(p) == 2
    end
end


text:
module NARemovers

using Random
using DataFrames

using AutoMLPipeline.AbsTypes
using AutoMLPipeline.BaseFilters
using AutoMLPipeline.Utils

import AutoMLPipeline.AbsTypes: fit!, transform!
export fit!, transform!
export NARemover

"""
    NARemover(
      Dict(
        :name => "nadetect",
        :acceptance => 0.10 # tolerable NAs percentage
      )
    )

Removes columns with NAs greater than acceptance rate.
This assumes that it processes columns of features.
The output column should not be part of input to avoid
it being excluded if it fails the acceptance critera.

Implements `fit!` and `transform!`.
"""
mutable struct NARemover <: Transformer
  name::String
  model::Dict
  args::Dict

  function NARemover(args::Dict = Dict())
	 default_args = Dict(
								:name => "nadetect",
								:acceptance => 0.10
								)
	 cargs=nested_dict_merge(default_args,args)
	 cargs[:name] = cargs[:name]*"_"*randstring(3)
	 new(cargs[:name],Dict(),cargs)
  end
end

"""
NARemover(acceptance::Float64)

Helper function for NARemover.
"""
NARemover(acceptance::Float64) = NARemover(Dict(:acceptance => acceptance))

"""
    fit!(nad::NARemover,features::DataFrame,labels::Vector=[])

Checks and exit of df is empty

# Arguments
- `nad::NARemover`: custom type
- `features::DataFrame`: input
- `labels::Vector=[]`:
"""
function fit!(nad::NARemover, features::DataFrame, labels::Vector=[])
  if features == DataFrame()
	 error("empty dataframe")
  end
  nad.model = nad.args
end

"""
    transform!(nad::NARemover,nfeatures::DataFrame)

Removes columns with NAs greater than acceptance rate.

# Arguments
- `nad::NARemover`: custom type
- `nfeatures::DataFrame`: input
"""
function transform!(nad::NARemover, nfeatures::DataFrame)
  features = deepcopy(nfeatures)
  if features == DataFrame()
	 error("empty dataframe")
  end
  sz = nrow(features)
  tol = nad.model[:acceptance]
  colnames = []
  for (colname,dat) in collect(pairs(eachcol(features)))
	 if sum(ismissing.(dat)) < tol*sz
		push!(colnames,colname)
	 end
  end
  xtr =  features[:,colnames]
  return xtr
end

end


text:
export CosSumGaussianSpin

"""
    SpinCosSumGaussian{D,V}( n_cos, n_gaussians, k, α, σ, μ, δ )

Data type for parameters of initial distribution

```math
(1+ \\cos( \\sum^{n_{cos}}_{i=1} k_i x))
\\cdot
\\sum_{j=1}^{n_{gaussians}}
\\delta_j
\\exp \\big( -\\frac{1}{2}
\\frac{(v-\\mu_j)^2}{\\sigma_j^2} \\big)
```

## Parameters

- `k` : values of the wave numbers (one array for each cosines)
- `α` : strength of perturbations
- `σ` : variance of the Gaussian (one velocity vector for each gaussian).
- `μ` : mean value of the Gaussian (one velocity vector for each gaussian).
- `δ` : portion of each Gaussian

# Example

```math
f(x,v_1,v_2)=\\frac{1}{2\\pi\\sigma_1\\sigma_2} \\exp \\Big( - \\frac{1}{2}
\\big( \\frac{v_1^2}{\\sigma_1^2} + \\frac{v_2^2}{\\sigma_2^2} \\big)
\\Big) ( 1 + \\alpha \\cos(kx)),
```

```julia
df = CosSumGaussian{1,1,3}([[k]],[α], [[σ₁,σ₂]], [[μ₁,μ₂]])
```

"""
struct CosSumGaussianSpin <: AbstractCosGaussian
    params::CosGaussianParams

    function CosSumGaussianSpin(
        k::Array{Vector{Float64},1},
        α::Vector{Float64},
        σ::Array{Vector{Float64},1},
        μ::Array{Vector{Float64},1},
        δ::Vector{Float64}=[1.0],
    )
        dims = (1, 1)
        params = CosGaussianParams(dims, k, α, σ, μ, δ)
        return new(params)
    end
end

function eval_x_density(f::CosSumGaussianSpin, x::Float64)
    fval = 1.0
    for j in 1:(f.params.n_cos)
        fval += f.params.α[j] * cos(sum(f.params.k[j] * x))
    end
    return fval
end

function (f::CosSumGaussianSpin)(x::Float64, v::Float64)
    return eval_x_density(f, x) * eval_v_density(f, v)
end


text:
# This file was generated, do not modify it. # hide
Statistics.mean(v...) = mean(v)
Statistics.mean(v::AbstractVector{<:AbstractNode}) = node(mean, v...)

yhat = mean([predict(m, Xs) for  m in machines]);


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct ApplicationGatewayBackendAddressPoolPropertiesFormat <: SwaggerModel
    backendIPConfigurations::Any # spec type: Union{ Nothing, Vector{NetworkInterfaceIPConfiguration2} } # spec name: backendIPConfigurations
    backendAddresses::Any # spec type: Union{ Nothing, Vector{ApplicationGatewayBackendAddress} } # spec name: backendAddresses
    provisioningState::Any # spec type: Union{ Nothing, ProvisioningState } # spec name: provisioningState

    function ApplicationGatewayBackendAddressPoolPropertiesFormat(;backendIPConfigurations=nothing, backendAddresses=nothing, provisioningState=nothing)
        o = new()
        validate_property(ApplicationGatewayBackendAddressPoolPropertiesFormat, Symbol("backendIPConfigurations"), backendIPConfigurations)
        setfield!(o, Symbol("backendIPConfigurations"), backendIPConfigurations)
        validate_property(ApplicationGatewayBackendAddressPoolPropertiesFormat, Symbol("backendAddresses"), backendAddresses)
        setfield!(o, Symbol("backendAddresses"), backendAddresses)
        validate_property(ApplicationGatewayBackendAddressPoolPropertiesFormat, Symbol("provisioningState"), provisioningState)
        setfield!(o, Symbol("provisioningState"), provisioningState)
        o
    end
end # type ApplicationGatewayBackendAddressPoolPropertiesFormat

const _property_map_ApplicationGatewayBackendAddressPoolPropertiesFormat = Dict{Symbol,Symbol}(Symbol("backendIPConfigurations")=>Symbol("backendIPConfigurations"), Symbol("backendAddresses")=>Symbol("backendAddresses"), Symbol("provisioningState")=>Symbol("provisioningState"))
const _property_types_ApplicationGatewayBackendAddressPoolPropertiesFormat = Dict{Symbol,String}(Symbol("backendIPConfigurations")=>"Vector{NetworkInterfaceIPConfiguration2}", Symbol("backendAddresses")=>"Vector{ApplicationGatewayBackendAddress}", Symbol("provisioningState")=>"ProvisioningState")
Base.propertynames(::Type{ ApplicationGatewayBackendAddressPoolPropertiesFormat }) = collect(keys(_property_map_ApplicationGatewayBackendAddressPoolPropertiesFormat))
Swagger.property_type(::Type{ ApplicationGatewayBackendAddressPoolPropertiesFormat }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_ApplicationGatewayBackendAddressPoolPropertiesFormat[name]))}
Swagger.field_name(::Type{ ApplicationGatewayBackendAddressPoolPropertiesFormat }, property_name::Symbol) =  _property_map_ApplicationGatewayBackendAddressPoolPropertiesFormat[property_name]

function check_required(o::ApplicationGatewayBackendAddressPoolPropertiesFormat)
    true
end

function validate_property(::Type{ ApplicationGatewayBackendAddressPoolPropertiesFormat }, name::Symbol, val)
end


text:
# takes in a folder with all sample results, create a dataset in

using Distributed

##
@everywhere using Pkg
@everywhere Pkg.activate(pwd())
Pkg.instantiate()
@everywhere using JLD2, FileIO
@everywhere using StatsBase
@everywhere using MyMCMCSampler

using PyCall
using Distributions
using DataFrames
using ArgParse
using ProgressMeter
using Base.Iterators: enumerate

py"""
import deepdish as dd
def savefile(filepath,data):
    dd.io.save(filepath,data)
"""
##
sa = ArgParseSettings()
@add_arg_table sa begin
    "--indir"
        help = "input folder containing all MCMC samples"
        arg_type = String
        default = "data/third/continuous/10small/samples"
    "--outdir"
        help = "folder to store dataset file"
        arg_type = String
        default = "data/tmp/"
    "--statistics"
        help = "which statistis to calculate"
        arg_type = String
        default = "standardized_moments"
    "--order"
        help = "how many moments to use"
        arg_type = Int
        default = 4
end

args = parse_args(ARGS, sa)
for i in procs()
    @spawnat i args |> fetch
end
indir = args["indir"]
outdir = args["outdir"]
@everywhere order = args["order"]

##
# get path for all graphs
filenames = readdir(indir)
graph_paths = [ indir*"/"*filename for filename in filenames]
@show "$(length(graph_paths)) graphs to be loaded"

# load samples and calculate statistics
if args["statistics"] == "standardized_moments"
    if order <= 2
        @everywhere stats = vcat(
            [mean],
            [x->abs(moment(x,i)).^(1/i).*sign(moment(x,i)) for i in 2:order]
        )
    else
        @everywhere stats = vcat(
            [mean, x->moment(x,2)],
            [x->abs(moment(x,i)).^(1/i) ./ (moment(x,2).^(0.5)) .* sign(moment(x,i)) for i in 3:order]
        )
    end
elseif args["statistics"] == "central_moments"
    @everywhere stats = vcat([mean], [x->abs(moment(x,i)).^(1/i) .* sign(moment(x,i)) for i in 2:order])
elseif args["statistics"] == "moments"
    @everywhere stats = [x-> abs(mean(x.^i)).^(1/i).*sign(mean(x.^i)) for i in 1:order];
end

@everywhere function load_and_compute(path)
    data = load(path)
    g = data["graph"]
    samples = data["samples"]  # (nsample, nvar, nchain)
    samples = vcat([samples[:,:,i] for i in 1:size(samples,3)]...) # (nsamples, nvar)
    @assert ndims(samples) == 2
    dict = dist2dict(g)
    dict["beta"] = data["beta"]
    dict["alpha"] = data["alpha"]
    dict["p"] = data["p"]
    dict["k2"] = data["k2"]
    dict["k3"] = data["k3"]
    nmoments = length(stats)
    nvar = size(samples,2)
    marginals = zeros(nvar, nmoments)
    for i in 1:nvar, j in 1:nmoments
        marginals[i,j]  = stats[j](samples[:,i])
    end
    return (dict, marginals)
end

results = @showprogress pmap(graph_paths) do path
    load_and_compute(path)
end

# show statistics
M = vcat([item[2] for item in results]...)
for i in 1:order
    println("------------ Statistics $i ------------")
    println(summarystats(M[:,i]))
end

## save to file
if  !isdir(outdir)
    println("$(outdir) doesn't exist, creating")
    mkpath(outdir)
end
outpath = joinpath(outdir, "data.dd")
@show "Saving dataset to $outpath"
py"savefile"(outpath, results)
@show "File saved"