text:
# GaussDCA
# ========

"""
Wrapper function to `GaussDCA.gDCA`.
You need to install GaussDCA:

```julia
using Pkg
Pkg.add(PackageSpec(url="https://github.com/carlobaldassi/GaussDCA.jl"))
```
Look into [GaussDCA.jl README](https://github.com/carlobaldassi/GaussDCA.jl) for further information.
If you use this wrapper, **please cite the GaussDCA publication and the package's doi**.

It's possible to indicate the path to the julia binary where GaussDCA is installed.
However, it's recommended to use the same version where MIToS is installed. That is
because this function use `serialize`/`deserialize` to transfer data between the processes.

**GaussDCA Publication:**
Baldassi, Carlo, Marco Zamparo, Christoph Feinauer, Andrea Procaccini, Riccardo Zecchina, Martin Weigt, and Andrea Pagnani.
"Fast and accurate multivariate Gaussian modeling of protein families: predicting residue contacts and protein-interaction partners."
PloS one 9, no. 3 (2014): e92721.
"""
function gaussdca(msa; juliapath::String=joinpath(Sys.BINDIR,Base.julia_exename()), kargs...)
    base_name = tempname()
    if Sys.iswindows()
        base_name = escape_string(base_name)
        juliapath = escape_string(juliapath)
        if !endswith(juliapath, ".exe")
            juliapath = juliapath * ".exe"
        end
    end
    script_file = base_name * ".jl"
    msa_file = base_name * ".fasta"
    jdl_file = base_name * ".jls"
    plm = fill!(columnpairsmatrix(msa), NaN)
    write(msa_file, msa, FASTA)
    try
        _create_script(script_file, msa_file, jdl_file; kargs...)
        run(`$juliapath $script_file`)
        pairedvalues = open(deserialize, jdl_file, "r")
        for (i,j,value) in pairedvalues
           plm[i,j] = value
        end
    finally
        isfile(script_file) && rm(script_file)
        isfile(msa_file) && rm(msa_file)
        isfile(jdl_file) && rm(jdl_file)
    end
    plm
end

function _create_script(script_file::String, msa_file::String, jdl_file::String; kargs...)
    if length(kargs) > 0
        for (k,v) in kargs
            @assert isa(v,Number) || isa(v,Symbol) "Argument values must be numbers or symbols"
        end
        str_kargs = "," * join([ isa(v,Symbol) ? "$k=:$v" : "$k=$v" for (k,v) in kargs],',')
    else
        str_kargs = ""
    end
    open(script_file, "w") do fh
        println(fh, "using Serialization;")
        println(fh, "using GaussDCA;")
        println(fh, "values = GaussDCA.gDCA(\"$msa_file\"$str_kargs);")
        println(fh, "open(fh -> serialize(fh, values), \"$(jdl_file)\", \"w\")")
    end
end


text:
import Base: values, find

# when ############################

when{T,N}(ta::TimeArray{T,N}, period::Function, t::Int)         = ta[find(period(ta.timestamp) .== t)]
when{T,N}(ta::TimeArray{T,N}, period::Function, t::String) = ta[find(period(ta.timestamp) .== t)]

# from, to ######################

from{T,N,D}(ta::TimeArray{T,N,D}, d::D) =
    length(ta) == 0 ? ta :
         d < ta.timestamp[1] ? ta :
         d > ta.timestamp[end] ? ta[1:0] :
         ta[searchsortedfirst(ta.timestamp, d):end]

to{T,N,D}(ta::TimeArray{T,N,D}, d::D) =
    length(ta) == 0 ? ta :
        d < ta.timestamp[1] ? ta[1:0] :
        d > ta.timestamp[end] ? ta :
        ta[1:searchsortedlast(ta.timestamp, d)]

###### find ##################

find(ta::TimeArray{Bool,1}) = find(ta.values)

###### findwhen #################

findwhen(ta::TimeArray{Bool,1}) = ta.timestamp[find(ta.values)]

###### head, tail ###########

function head{T,N,D}(ta::TimeArray{T,N,D}, n::Int=1)
    ncol          = length(ta.colnames)
    new_timestamp = ta.timestamp[1:n]
    new_values    = ta.values[1:n, 1:ncol]
    TimeArray(new_timestamp, new_values, ta.colnames, ta.meta)
end

function tail{T,N,D}(ta::TimeArray{T,N,D}, n::Int=1)
    ncol          = length(ta.colnames)
    tail_start = length(ta)-n+1
    new_timestamp = ta.timestamp[tail_start:end]
    new_values    = ta.values[tail_start:end, 1:ncol]
    TimeArray(new_timestamp, new_values, ta.colnames, ta.meta)
end

###### element wrapers ###########

timestamp{T,N}(ta::TimeArray{T,N}) = ta.timestamp
values{T,N}(ta::TimeArray{T,N})    = ta.values
colnames{T,N}(ta::TimeArray{T,N})  = ta.colnames
meta{T,N}(ta::TimeArray{T,N})      = ta.meta

ext: jl


text:
ENV["GKSwstype"] = "100"  # set 'GR environment' to 'no output' (for Travis CI)
using Documenter, LazySets
import Polyhedra, Optim, Expokit, TaylorModels, Distributions

include("init.jl")

makedocs(
    sitename = "LazySets.jl",
    modules = [LazySets, Approximations, Arrays, LazySets.Parallel],
    format = Documenter.HTML(
        prettyurls = get(ENV, "CI", nothing) == "true",
        assets = ["assets/juliareach.css"]),
    pages = [
        "Home" => "index.md",
        "Manual" => Any[
        "Getting Started" => "man/getting_started.md",
        "Polyhedral Approximations" => "man/polyhedral_approximations.md",
        "Decomposing an Affine Map" => "man/decompose_example.md",
        "Fast 2D LPs" => "man/fast_2d_LPs.md",
        "Iterative Refinement" => "man/iterative_refinement.md",
        "Interval Hulls" => "man/interval_hulls.md",
        "Convex Hulls" => "man/convex_hulls.md",
        "Operations on Sets" => "man/set_operations.md",
        "A Reachability Algorithm" => "man/reach_zonotopes.md",
        "A Hybrid Reachability Algorithm" => "man/reach_zonotopes_hybrid.md",
        "Concrete Polyhedra" => "man/concrete_polyhedra.md",
        "Parallel Approximations" => "man/parallel_approximations.md",
        "Lazy Intersections" => "man/lazy_intersections.md"
        ],
        "Library" => Any[
        "Set Interfaces" => "lib/interfaces.md",
        "Common Set Representations" => "lib/representations.md",
        "Common Set Operations" => "lib/operations.md",
        "Comparisons" => "lib/comparisons.md",
        "Conversions between set representations" => "lib/conversion.md",
        "Binary Functions on Sets" => "lib/binary_functions.md",
        "Approximations" => "lib/approximations.md",
        "Utility Functions" => "lib/utils.md",
        "Parallel" => "lib/parallel.md"
        ],
        "About" => "about.md"
    ],
    doctest = false,
    strict = true
)

deploydocs(
    repo = "github.com/JuliaReach/LazySets.jl.git"
)


text:
include("environment.jl")
include("learnedheuristic.jl")
include("rewards/rewards.jl")
include("utils.jl")

LearnedHeuristic(agent::RL.Agent) = LearnedHeuristic{DefaultStateRepresentation{DefaultFeaturization, DefaultTrajectoryState}, DefaultReward, FixedOutput}(agent)


text:
using Test

m = page_model()
include("../src/components/SulphateForcing.jl")

add_comp!(m, SulphateForcing)

p = load_parameters(m)
p["y_year_0"] = 2008.
p["y_year"] = Mimi.dim_keys(m.md, :time)
set_leftover_params!(m, p)

run(m)

forcing=m[:SulphateForcing,:fs_sulphateforcing]
forcing_compare=readpagedata(m,"test/validationdata/fs_sulphateforcing.csv")

@test forcing ≈ forcing_compare rtol=1e-3


text:
using MAT
using LinearAlgebra
using SparseArrays
using MatrixNetworks
using Clustering

## Load Data
mat = matread("../data/Florida_Bay_Dataset.mat")
A125 = mat["A125"]                    # This is the processed dataset--a directed who-eats-who graph
n = size(A125,1)
C = round.(Int64,mat["C"])            # These are four diffrent ways to classify the species
label = C[:,3]                        # These is the labels of Li et al. NeurIPs 2017, ignoring "predator" and "no predator" tags
Li_c = vec(round.(Int64,mat["Li_c"])) # This is the clustering obtained by Li et al.

## Form different graphs based on different reductions
include("../src/Motif_Adj_Mat.jl")
AFan = ExtractFanMat(A125)      # Motif-adjacency matrix
ASim = turntounweighted(A125)   # Symmetrized version of directed graph

# Save, to run other algorithms requiring MATLAB interface
matwrite("../data/Fl_bay_Afan.mat",Dict("AFan"=>AFan))

# Scale the graphs so that their total volume is the same
mSim = sum(ASim)
mFan = sum(AFan)
AFan = mSim/mFan*AFan
m = mSim

# Get weighted degree vectors for each graph
dSim = vec(sum(ASim,dims = 2))
dFan = vec(sum(AFan,dims = 2))

## Run Louvain-based algorithms
include("../src/Graph_Louvain.jl")
numpoints = 20  # number of lambda values to try

# For each graph, store the ari scores and number of clusters
SimLouv = zeros(numpoints,2)
FanLouv= zeros(numpoints,2)

Lams = range(0.05/m,5/m,length=numpoints)
SimC  = zeros(n,numpoints)
FanC = zeros(n,numpoints)

numtimes = 1          # number of times to run Louvain.
totaltimes = zeros(length(Lams),2)  # store runtimes

for i = 1:length(Lams)
    lam = Lams[i]
    println(" Lambda = $lam")
    # Run the standard Lambda-Louvain algorithm (i.e., the motif is an edge)
    start = time()
    cSim, BestObj = Many_Louvain(ASim,dSim,lam,numtimes)
    totaltimes[i,1] = time()-start
    SimC[:,i] = cSim
    sim_ari = ari(cSim,label)
    SimLouv[i,:] = [sim_ari length(unique(cSim))]

    # Run the Lambda-Louvain algorithm for the bi-fan motif
    start = time()
    cFan, BestObj = Many_Louvain(AFan,dFan,lam,numtimes)
    totaltimes[i,2] = time()-start
    FanC[:,i] = cFan
    fan_ari = ari(cFan,label)
    FanLouv[i,:] = [fan_ari length(unique(cFan))]
end

p1 = sortperm(SimLouv[:,2])
p2 = sortperm(FanLouv[:,2])
SimLouv = SimLouv[p1,:]
FanLouv = FanLouv[p2,:]
matwrite("Output/Florida_Bay_125_output.mat", Dict("FanC" => FanC,
    "FanLouv" => FanLouv,"Lams" => collect(Lams),"totaltimes"=>totaltimes))

## Extract the graclus clusterings
G = matread("Output/Graclus_Flbay.mat")
GracC = round.(Int64,G["GracFan"])
num = size(GracC,2)
Grac = zeros(num,2)
GracRun = G["times"][:,2]
for i = 1:num
    cG = vec(GracC[:,i])
    Grac[i,1] = ari(cG,label)
    Grac[i,2] = length(unique(cG))
end

## Extract the Metis clusterings
G = matread("Output/Metis_Flbay.mat")
MetC = round.(Int64,G["MetFan"])
num = size(MetC,2)
Met = zeros(num,2)
MetRun = G["times"][:,2]
for i = 1:num
    cM = vec(MetC[:,i])
    Met[i,1] = ari(cM,label)
    Met[i,2] = length(unique(cM))
end

## Extract the hMetis clusterings
G = matread("Output/hmetis_2to50_10.mat")
hMetC = round.(Int64,G["hmetC"])
num = size(hMetC,2)
hMet = zeros(num,2)
hMetRun = G["times"]
for i = 1:num
    cM = vec(hMetC[:,i])
    hMet[i,1] = ari(cM,label)
    hMet[i,2] = length(unique(cM))
end
# order by cluster size--this isn't always automatic
p = sortperm(hMet[:,2])
hMet = hMet[p,:]
## Extract the Recursive Spectral Clusterings
G = matread("Output/RecursiveSpectral_Flbay.mat")
SpecC = round.(Int64,G["SpecFan"])
num = size(MetC,2)
Spec = zeros(num,2)
SpecRun = G["times"][:,2]
for i = 1:num
    cM = vec(SpecC[:,i])
    Spec[i,1] = ari(cM,label)
    Spec[i,2] = length(unique(cM))
end

p = sortperm(Spec[:,2])
Spec = Spec[p,:]
## Plot results
using Plots

lw = 3
msw = 0
s1 = 300
s2 = 250
ms = 5
lw = 2

rightmax = 51

p = plot(FanLouv[:,2], FanLouv[:,1],color=:green,markerstrokecolor=:green,
markerstrokewidth = msw, markersize = ms,markershape = :circle,linewidth= lw, label = "HypLam")

plot!(SimLouv[:,2], SimLouv[:,1],color=:black,markerstrokecolor=:black,
markerstrokewidth = msw, markersize = ms,markershape = :circle,linewidth= lw, xlim = [0, rightmax],grid = false, label = "LamLouv",xaxis = "Number of Clusters")

plot!(Grac[:,2], Grac[:,1],color=:blue,markerstrokecolor=:blue,
markerstrokewidth = msw, markersize = ms,markershape = :circle,linewidth= lw,ylabel = "ARI Scores",size = (s1,s2), label = "Graclus")

plot!(Met[:,2], Met[:,1],color=:orange,markerstrokecolor=:orange,
markerstrokewidth = msw, markersize = ms,markershape = :circle,linewidth= lw,ylabel = "ARI Scores",size = (s1,s2), label = "Metis")

plot!(hMet[:,2], hMet[:,1],color=:purple,markerstrokewidth = msw, markerstrokecolor=:purple,
markersize = ms,markershape = :circle,linewidth= lw,ylabel = "ARI Scores",size = (s1,s2), label = "hMetis")

plot!(Spec[:,2], Spec[:,1],color=:brown,markerstrokewidth = msw, markersize = ms,markershape = :circle,linewidth= lw,ylabel = "ARI Scores",size = (s1,s2),label = "Spectral",
ylim = [0,.5],foreground_color_legend = nothing,background_color_legend=nothing,markerstrokecolor=:brown)

savefig(p,"plots/Florida_Bay_ari.pdf")

## Runtime plots
p = plot(FanLouv[:,2], totaltimes[:,2],color=:green,linewidth= lw, label = "")
scatter!(FanLouv[:,2],totaltimes[:,2],color=:green,markerstrokecolor=:green,markerstrokewidth = msw, markersize = ms, label = "HyperLam")

plot!(SimLouv[:,2], totaltimes[:,1],color=:black,linewidth= lw, xlim = [0, rightmax],grid = false, label = "",xaxis = "Number of Clusters")
scatter!(SimLouv[:,2],totaltimes[:,1],color=:black,markerstrokecolor=:black,markerstrokewidth = msw, markersize = ms, label = "LamLouvain",markershape= :circle)

plot!(Grac[:,2], GracRun,color=:blue,linewidth= lw, label = "",ylabel = "Runtime (seconds)",size = (s1,s2))
scatter!(Grac[:,2], GracRun,color=:blue,markerstrokecolor=:blue,markerstrokewidth = msw, markersize = ms, label = "Graclus")

plot!(Met[:,2], MetRun,color=:orange,linewidth= lw, label = "",size = (s1,s2),legend = false)
scatter!(Met[:,2], MetRun,color=:orange,markerstrokecolor=:orange,markerstrokewidth = msw, markersize = ms, label = "Metis")

plot!(hMet[:,2], hMetRun,color=:purple,linewidth= lw, label = "",size = (s1,s2),yaxis = :log10)
scatter!(hMet[:,2], hMetRun,color=:purple,markerstrokewidth = msw,
markersize = ms, label = "hMetis",markerstrokecolor=:purple)

plot!(Spec[:,2], SpecRun,color=:brown,linewidth= lw, label = "",size = (s1,s2),yaxis = :log10)
scatter!(Spec[:,2], SpecRun,color=:brown,markerstrokewidth = msw, markersize = ms, label = "Spectral",
foreground_color_legend = nothing,background_color_legend=nothing,markerstrokecolor=:brown)

savefig(p,"plots/Florida_Bay_runtime.pdf")

## ARI/Runtime Tradeoff Plots
ms = 5
p = scatter(totaltimes[:,2],FanLouv[:,1],color=:green,linewidth= lw,markerstrokecolor=:green,
markerstrokewidth = msw, markersize = ms, label = "HyperLam", legend = false)

scatter!(totaltimes[:,1],SimLouv[:,1],color=:black,linewidth= lw,
grid = false,xaxis = "Runtime", yaxis = "ARI Scores",markerstrokecolor=:black,
markerstrokewidth = msw, markersize = ms, label = "LamLouvain",markershape= :circle)

scatter!(GracRun,Grac[:,1],color=:blue,linewidth= lw,markerstrokecolor=:blue,
size = (s1,s2),markerstrokewidth = msw, markersize = ms, label = "Graclus")

scatter!(MetRun,Met[:,1],color=:orange,linewidth= lw,size = (s1,s2),markerstrokecolor=:orange,
legend = false,markerstrokewidth = msw, markersize = ms, label = "Metis")

scatter!(hMetRun,hMet[:,1],color=:purple,linewidth= lw,markerstrokecolor=:purple,
markerstrokewidth = msw, markersize = ms, label = "hMetis")

scatter!(SpecRun, Spec[:,1], color=:brown,linewidth= lw,size = (s1,s2),xaxis = :log10,
markerstrokewidth = msw, markersize = ms, label = "Spectral",markerstrokecolor=:brown,
foreground_color_legend = nothing,background_color_legend=nothing)

## Increase size of largest ARI dot
ms = 10
best = findmax(FanLouv[:,1])[2]
scatter!([totaltimes[best,2]],[FanLouv[best,1]],color=:green,linewidth= lw,markerstrokecolor=:green,
markerstrokewidth = msw, markersize = ms, label = "HyperLam", legend = false)

best = findmax(SimLouv[:,1])[2]
scatter!([totaltimes[best,1]],[SimLouv[best,1]],color=:black,linewidth= lw,
grid = false,xaxis = "Runtime", yaxis = "ARI Scores",markerstrokecolor=:black,
markerstrokewidth = msw, markersize = ms, label = "LamLouvain",markershape= :circle)

best = findmax(Grac[:,1])[2]
scatter!([GracRun[best]],[Grac[best,1]],color=:blue,linewidth= lw,markerstrokecolor=:blue,
size = (s1,s2),markerstrokewidth = msw, markersize = ms, label = "Graclus")

best = findmax(Met[:,1])[2]
scatter!([MetRun[best]],[Met[best,1]],color=:orange,linewidth= lw,size = (s1,s2),markerstrokecolor=:orange,
legend = false,markerstrokewidth = msw, markersize = ms, label = "Metis")

best = findmax(hMet[:,1])[2]
scatter!([hMetRun[best]],[hMet[best,1]],color=:purple,linewidth= lw,markerstrokecolor=:purple,
markerstrokewidth = msw, markersize = ms, label = "hMetis")

best = findmax(Spec[:,1])[2]
scatter!([SpecRun[best]], [Spec[best,1]], color=:brown,linewidth= lw,size = (s1,s2),xaxis = :log10,
markerstrokewidth = msw, markersize = ms, label = "Spectral",markerstrokecolor=:brown,
foreground_color_legend = nothing,background_color_legend=nothing)

savefig(p,"plots/Florida_Bay_tradeoff.pdf")

max_stars_repo_name: nveldt/ParamCC

max_issues_repo_name: nveldt/ParamCC

max_forks_repo_name: nveldt/ParamCC


text:
# ---
# title: 330. Patching Array
# id: problem330
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Greedy
# link: <https://leetcode.com/problems/patching-array/description/>
# hidden: true
# ---
#
# Given a sorted positive integer array _nums_ and an integer _n_ , add/patch
# elements to the array such that any number in range `[1, n]` inclusive can be
# formed by the sum of some elements in the array. Return the minimum number of
# patches required.
#
# **Example 1:**
#
#
#
#     Input: _nums_ = [1,3], _n_ = 6
#     Output: 1
#     Explanation:
#     Combinations of _nums_ are [1], [3], [1,3], which form possible sums of: 1, 3, 4.
#     Now if we add/patch 2 to _nums_ , the combinations are: [1], [2], [3], [1,3], [2,3], [1,2,3].
#     Possible sums are 1, 2, 3, 4, 5, 6, which now covers the range [1, 6].
#     So we only need 1 patch.
#
# **Example 2:**
#
#
#
#     Input: _nums_ = [1,5,10], _n_ = 20
#     Output: 2
#     Explanation: The two patches can be [2, 4].
#
#
# **Example 3:**
#
#
#
#     Input: _nums_ = [1,2,2], _n_ = 5
#     Output: 0
#
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
module SampleListTest

using Test
using ForneyLab
using ForneyLab: outboundType, isApplicable, prod!, unsafeMean, unsafeCov, unsafeVar, dims, bootstrap
using ForneyLab: SPSampleListOutNPP, VBSampleListOut

@testset "dims" begin
    @test dims(ProbabilityDistribution(Univariate, SampleList, s=[0.0, 1.0], w=[0.5, 0.5])) == ()
    @test dims(ProbabilityDistribution(Multivariate, SampleList, s=[[0.0], [1.0]], w=[0.5, 0.5])) == (1,)
    @test dims(ProbabilityDistribution(MatrixVariate, SampleList, s=[mat(0.0), mat(1.0)], w=[0.5, 0.5])) == (1,1)
end

@testset "unsafeMean" begin
    @test unsafeMean(ProbabilityDistribution(Univariate, SampleList, s=[1.0, 1.0], w=[0.5, 0.5])) == 1.0
    @test unsafeMean(ProbabilityDistribution(Multivariate, SampleList, s=[[1.0], [1.0]], w=[0.5, 0.5])) == [1.0]
    @test unsafeMean(ProbabilityDistribution(MatrixVariate, SampleList, s=[mat(1.0), mat(1.0)], w=[0.5, 0.5])) == mat(1.0)
end

@testset "unsafeVar" begin
    @test unsafeVar(ProbabilityDistribution(Univariate, SampleList, s=[1.0, 1.0], w=[0.5, 0.5])) == 0.0
    @test unsafeVar(ProbabilityDistribution(Multivariate, SampleList, s=[[1.0], [1.0]], w=[0.5, 0.5])) == [0.0]
end

@testset "unsafeCov" begin
    @test unsafeCov(ProbabilityDistribution(Univariate, SampleList, s=[1.0, 1.0], w=[0.5, 0.5])) == 0.0
    @test unsafeCov(ProbabilityDistribution(Multivariate, SampleList, s=[[1.0], [1.0]], w=[0.5, 0.5])) == mat(0.0)
    @test unsafeCov(ProbabilityDistribution(MatrixVariate, SampleList, s=[eye(2), eye(2)], w=[0.5, 0.5])) == zeros(4,4)
end

@testset "Univariate SampleList prod!" begin
    dist_gaussian = ProbabilityDistribution(Univariate, GaussianMeanVariance, m=0.0, v=1.0)
    dist_prod = dist_gaussian * ProbabilityDistribution(Univariate, SampleList, s=[0.0, 1.0], w=[0.5, 0.5])
    dist_true_prod = ProbabilityDistribution(Univariate, SampleList, s=[0.0, 1.0], w=[0.6224593312018546, 0.37754066879814546])

    @test dist_prod.params[:s] == dist_true_prod.params[:s]
    @test dist_prod.params[:w] == dist_true_prod.params[:w]
    @test dist_prod.params[:logintegrand]([-0.7, 1.5]) == logPdf.([dist_gaussian],[-0.7, 1.5])
end

@testset "Multivariate SampleList prod!" begin
    dist_gaussian = ProbabilityDistribution(Multivariate, GaussianMeanVariance, m=[0.0], v=mat(1.0))
    dist_prod = dist_gaussian * ProbabilityDistribution(Multivariate, SampleList, s=[[0.0], [1.0]], w=[0.5, 0.5])
    dist_true_prod = ProbabilityDistribution(Multivariate, SampleList, s=[[0.0], [1.0]], w=[0.6224593312018546, 0.37754066879814546])

    @test dist_prod.params[:s] == dist_true_prod.params[:s]
    @test dist_prod.params[:w] == dist_true_prod.params[:w]
    @test dist_prod.params[:logintegrand]([[-0.7], [1.5]]) == logPdf.([dist_gaussian],[[-0.7], [1.5]])
end

@testset "Initialization of logproposal, logintegrand, unnormalizedweights" begin
    dist_beta = ProbabilityDistribution(Univariate, Beta, a=2.0, b=1.0)
    dist_gamma = ProbabilityDistribution(Univariate, Gamma, a=2.0, b=1.0)
    dist_sample = ProbabilityDistribution(Univariate, SampleList, s=[0.0, 1.0], w=[0.5, 0.5])
    dist_betagamma = dist_beta*dist_gamma
    dist_gammasample = dist_gamma*dist_sample
    dist_betagg = dist_betagamma*dist_gamma

    @test haskey(dist_betagamma.params, :logproposal)
    @test haskey(dist_betagamma.params, :logintegrand)
    @test haskey(dist_betagamma.params, :unnormalizedweights)

    @test !haskey(dist_gammasample.params, :logproposal)
    @test haskey(dist_gammasample.params, :logintegrand)
    @test !haskey(dist_gammasample.params, :unnormalizedweights)

    @test haskey(dist_betagg.params, :logproposal)
    @test haskey(dist_betagg.params, :logintegrand)
    @test haskey(dist_betagg.params, :unnormalizedweights)
end

@testset "bootstrap" begin
    p1 = ProbabilityDistribution(Univariate, SampleList, s=[2.0], w=[1.0])
    p2 = ProbabilityDistribution(Univariate, PointMass, m=0.0)
    @test bootstrap(p1, p2) == [2.0]

    p1 = ProbabilityDistribution(Multivariate, SampleList, s=[[2.0]], w=[1.0])
    p2 = ProbabilityDistribution(MatrixVariate, PointMass, m=mat(tiny))
    @test isapprox(bootstrap(p1, p2)[1][1], 2.0, atol=1e-4)

    p1 = ProbabilityDistribution(Univariate, GaussianMeanVariance, m=2.0, v=0.0)
    p2 = ProbabilityDistribution(Univariate, SampleList, s=[0.0], w=[1.0])
    @test bootstrap(p1, p2) == [2.0]

    p1 = ProbabilityDistribution(Multivariate, GaussianMeanVariance, m=[2.0], v=mat(0.0))
    p2 = ProbabilityDistribution(MatrixVariate, SampleList, s=[mat(tiny)], w=[1.0])
    @test isapprox(bootstrap(p1, p2)[1][1], 2.0, atol=1e-4)
end

#-------------
# Update rules
#-------------

@testset "SPSampleListOutNPP" begin
    @test SPSampleListOutNPP <: SumProductRule{SampleList}
    @test outboundType(SPSampleListOutNPP) == Message{SampleList}
    @test isApplicable(SPSampleListOutNPP, [Nothing, Message{PointMass}, Message{PointMass}])

    @test ruleSPSampleListOutNPP(nothing, Message(Multivariate, PointMass, m=[0.0,2.0,5.2]), Message(Multivariate, PointMass, m=[0.4,0.5,0.1])) == Message(Univariate, SampleList, s=[0.0,2.0,5.2], w=[0.4,0.5,0.1])
    @test ruleSPSampleListOutNPP(nothing, Message(Multivariate, PointMass, m=[[0.0,2.0],[5.2,-0.6]]), Message(Multivariate, PointMass, m=[0.4,0.6])) == Message(Multivariate, SampleList, s=[[0.0,2.0],[5.2,-0.6]], w=[0.4,0.6])
    @test ruleSPSampleListOutNPP(nothing, Message(Multivariate, PointMass, m=[eye(2),eye(2)]), Message(Multivariate, PointMass, m=[0.4,0.6])) == Message(MatrixVariate, SampleList, s=[eye(2),eye(2)], w=[0.4,0.6])
end

@testset "VBSampleListOut" begin
    @test VBSampleListOut <: NaiveVariationalRule{SampleList}
    @test outboundType(VBSampleListOut) == Message{SampleList}
    @test isApplicable(VBSampleListOut, [Nothing, ProbabilityDistribution, ProbabilityDistribution])

    @test ruleVBSampleListOut(nothing, ProbabilityDistribution(Multivariate, PointMass, m=[0.0,2.0,5.2]), ProbabilityDistribution(Multivariate, PointMass, m=[0.4,0.5,0.1])) == Message(Univariate, SampleList, s=[0.0,2.0,5.2], w=[0.4,0.5,0.1])
    @test ruleVBSampleListOut(nothing, ProbabilityDistribution(Multivariate, PointMass, m=[[0.0,2.0],[5.2,-0.6]]), ProbabilityDistribution(Multivariate, PointMass, m=[0.4,0.6])) == Message(Multivariate, SampleList, s=[[0.0,2.0],[5.2,-0.6]], w=[0.4,0.6])
    @test ruleVBSampleListOut(nothing, ProbabilityDistribution(Multivariate, PointMass, m=[eye(2),eye(2)]), ProbabilityDistribution(Multivariate, PointMass, m=[0.4,0.6])) == Message(MatrixVariate, SampleList, s=[eye(2),eye(2)], w=[0.4,0.6])
end

#------------
# Integration
#------------

g(x) = x

@testset "Sampling code generation" begin
    # Define a model
    fg = FactorGraph()
    N = 4
    @RV s ~ SampleList(collect(1.0:N), ones(N)/N)
    @RV x ~ GaussianMeanVariance(s, 0.0)
    @RV y ~ Nonlinear{Sampling}(x, g=g)

    # Define an algorithm
    pfz = PosteriorFactorization(fg)
    algo = messagePassingAlgorithm(y)
    code = algorithmSourceCode(algo)

    @test occursin("ruleSPSampleListOutNPP", code)
    @test occursin("ruleSPGaussianMeanVarianceOutNSP", code)
    @test occursin("ruleSPNonlinearSOutNM", code)
end

# Generated algorithm code
function step!(data::Dict, marginals::Dict=Dict(), messages::Vector{Message}=Array{Message}(undef, 3))
    messages[1] = ruleSPSampleListOutNPP(nothing, Message(Multivariate, PointMass, m=[1.0, 2.0, 3.0, 4.0]), Message(Multivariate, PointMass, m=[0.25, 0.25, 0.25, 0.25]))
    messages[2] = ruleSPGaussianMeanVarianceOutNSP(nothing, messages[1], Message(Univariate, PointMass, m=0.0))
    messages[3] = ruleSPNonlinearSOutNM(g, nothing, messages[2])

    marginals[:y] = messages[3].dist

    return marginals
end

@testset "Sampling code execution" begin
    # Execute the algorithm code
    marginals = step!(Dict())

    @test marginals[:y] == ProbabilityDistribution(Univariate, SampleList, s=collect(1.0:4.0), w=ones(4)/4)
end

@testset "Differential Entropy estimate" begin
    fg = FactorGraph()

    @RV x ~ Gamma(3.4, 1.2)
    @RV s1 ~ Nonlinear{Sampling}(x, g=g)
    @RV y1 ~ Poisson(s1)
    @RV y2 ~ Poisson(x)

    @RV z ~ Gamma(3.4, 1.2)
    @RV s2 ~ Nonlinear{Sampling}(z, g=g)
    @RV y3 ~ Poisson(z)
    @RV y4 ~ Poisson(s2)

    @RV w ~ Gamma(3.4, 1.2)
    @RV y5 ~ Poisson(w)
    @RV y6 ~ Poisson(w)

    placeholder(y1,:y1)
    placeholder(y2,:y2)
    placeholder(y3,:y3)
    placeholder(y4,:y4)
    placeholder(y5,:y5)
    placeholder(y6,:y6)

    pfz = PosteriorFactorization(fg)
    algo = messagePassingAlgorithm([x,z,w])
    source_code = algorithmSourceCode(algo)
    eval(Meta.parse(source_code));

    v_1, v_2 = 7, 5
    data = Dict(:y1 => v_1, :y2 => v_2, :y3 => v_1, :y4 => v_2, :y5 => v_1, :y6 => v_2)
    marginals = step!(data)

    @test isa(marginals[:x].params[:entropy], Float64)
    @test isa(marginals[:z].params[:entropy], Float64)
end

end


text:
"""
Two ways of implementing options:

<something, ?> OR ?[<something> and <something else>]

The latter is more powerful, but a bit trickier to parse.
For now we'll go with the former, since you can express a more complex optional
section as a separate rule.
"""

"""
Defines the RegEx pattern used to embed grammar rules within one another.
"""
const rule_regex = r"<[\w_]+(,\s*[\?|\+|\*])?>"

"""
The base type for all grammar rules.
"""
abstract Rule

type Grammar
  rules::Dict{AbstractString, Rule}
  # caching the root rule might reduce the number of look-ups.

  Grammar() = new(Dict{AbstractString, Rule}())
end

"""
Parses a single option for a grammar rule and returns a Rule object encoding
that rule.
"""
rule(g::Grammar, r::AbstractString) =
  search(r, rule_regex) != 0:-1 ? NonTerminalRule(g, r) : TerminalRule(parse(r))

"""
Creates a modified version of a given grammar rule, according to some provided
modifier string. This function is used to implement EBNF constructs, such as
ONE-OR-MORE, ZERO-OR-MORE and OPTIONAL.
"""
function rule(r::Rule, modifier::AbstractString)
  if      modifier == "+"
    OneOrMoreRule(r)
  elseif  modifier == "*"
    ZeroOrMoreRule(r)
  elseif  modifier == "?"
    OptionalRule(r)
  else
    error("Unrecognised rule modifier: $(modifier).")
  end
end

"""
Creates a single grammar rule from a provided rule definition string.
"""
rule(g::Grammar, name::AbstractString, def::AbstractString) =
  g.rules[name] = rule(g, def)

"""
Creates a named OR rule for a given grammar, using a list of possible options,
each given by a rule definition string.
"""
rule(g::Grammar, name::AbstractString, defs...) =
  g.rules[name] = OrRule(Rule[rule(g, def) for def in defs])

"""
A rule reference is used to indirectly point to another rule, via its name.
"""
immutable RuleReference <: Rule
  name::AbstractString
end

"""
An OR rule allows a given grammar rule to be interpreted in a number of
different ways.
"""
immutable OrRule <: Rule
  num_options::Int
  options::Vector{Rule}

  OrRule(options::Vector{Rule}) = new(length(options), options)
end

"""
The ONE-OR-MORE rule (<x, +>) allows a grammar rule to be applied more than
once and at least once.
"""
immutable OneOrMoreRule <: Rule
  rule::Rule
end

"""
The ZERO-OR-MORE rule (<x, *>) allows a grammar rule to be applied
an arbitrary number of times, or not applied at all.
"""
immutable ZeroOrMoreRule <: Rule
  rule::Rule
end

"""
The OPTIONAL rule (<x, ?>) allows a grammar rule to be optionally applied at
most once.
"""
immutable OptionalRule <: Rule
  rule::Rule
end

"""
Due to Julia's current implementation of anonymous functions, this type isn't
quite as efficient as it could be.
"""
immutable NonTerminalRule <: Rule
  non_terminals::Vector{Rule}
  builder::Function # (Grammar, NonTerminalRule, Task) -> Expr

  NonTerminalRule(nts::Vector{Rule}, builder::Function) =
    new(nts, builder)

  function NonTerminalRule(g::Grammar, r::AbstractString)
    # Create an array to hold the non-terminal symbols within the given rule,
    # and another to hold their associated (possibly modified) rules.
    tags = AbstractString[]
    nts = Rule[]

    # Replace each grammar symbol in the given string with an associated
    # placeholder tag.
    function replacer(tag::AbstractString)
      push!(tags, string(tag))
      "__WALLACE_GRAMMAR_TAG_$(length(tags))__"
    end
    r = replace(r, rule_regex, replacer)

    # Now parse the rule definition to a Julia expression and then back into
    # Julia code (in the form a string).
    r = expr_to_def_s(Base.parse(r))

    # Replace each placeholder with a call to the appropriate function, and
    # build the list of rules used by the symbols within this non-terminal.
    for (i, tag) in enumerate(tags)
      parts = split(tag[2 : (length(tag) - 1)], r",\s+")
      tag_rule_name, modifier = parts[1], Base.get(parts, 2, "")

      # Create a reference to the rule for the symbol.
      tag_rule = RuleReference(tag_rule_name)

      # Apply any modifiers to the rule for this symbol.
      tag_rule = modifier == "" ? tag_rule : rule(tag_rule, modifier)

      push!(nts, tag_rule)

      # Replace the placeholder for this symbol within the Expr object with a
      # call to derive, along with the appropriate rule.
      r = replace(r, ":__WALLACE_GRAMMAR_TAG_$(i)__", "derive(g, r.non_terminals[$(i)], nxt)")
    end

    # Construct the lambda function for this non-terminal.
    r = "(g::Grammar, r::NonTerminalRule, nxt::Task) -> $(r)"
    builder = eval(Base.parse(r))

    # Build the rule object.
    new(nts, builder)
  end
end

"""
A terminal rule simply returns some predetermined constant symbol, expression,
string or integer.
"""
immutable TerminalRule <: Rule
  value

  TerminalRule(val) = new(val)
end

"""
Selects an interpretation of a given OR rule according to an index provided by
a given codon sequence.
"""
derive(g::Grammar, r::OrRule, nxt::Task) =
  derive(g, g.options[(consume(nxt) % num_options) + 1], nxt)

derive(g::Grammar, r::RuleReference, nxt::Task) =
  derive(g, g.rules[r.name], nxt)

derive(g::Grammar, r::TerminalRule, nxt::Task) =
  r.value

derive(g::Grammar, r::NonTerminalRule, nxt::Task) =
  r.builder(g, nxt)

derive(g::Grammar, nxt::Task) =
  derive(g, g.rules["root"], nxt)

"""
Converts a given Expr into a string containing Julia code capable of
reproducing that Expr.
"""
expr_to_def_s(ex::Expr) =
  "Expr($(join(vcat([":$(ex.head)"], map(expr_to_def_s, ex.args)), ", ")))"
expr_to_def_s(sym::Symbol) = ":$(sym)"
expr_to_def_s(a::Any) = string(a)

"""
TESTING!
"""
g = Grammar()
rule(g, "root", "<exp>")
# rule(g, "digit", 1:9) <-- Need to add support for Range next.
rule(g, "num", "blah!") # this one is the trickiest to deal with.
rule(g, "val", "x", "y", "<num>")
rule(g, "op", "<exp> * <exp>", "<exp> - <exp>", "<exp> + <exp>")
rule(g, "exp", "(<exp>)", "<val>", "<op>")


text:
#=
Stuff to do at the start of a develepment session:
using Profile
using Revise
using BenchmarkTools
] activate .
=#
using DelimitedFiles
using DataFrames
using CSV
using CPIMC
using CPIMC.Estimators
using CPIMC.PlaneWaves
using CPIMC.UniformElectronGas
using CPIMC.DefaultUpdates
using StaticArrays
using OnlineStats

#=
import CPIMC.DefaultUpdates: isuseful
function update!(m::Model, e::Ensemble, c::Configuration, updates)
    usefull_update_inds = findall(up -> isuseful(c, up), updates)
    @assert !isempty(usefull_update_inds)
    i = rand(usefull_update_inds)
    dv, Δ = updates[i](m, e, c)
    dv *= length(usefull_update_inds)/length(findall(up -> isuseful(apply_step(c, Δ), up), updates))
    if Δ == Step()
        return (i, :trivial)
    end

    if rand() < dv
        apply_step!(c, Δ)

        return (i, :accept)
    else
        return (i, :reject)
    end
end
=#

W_off_diag(e::Ensemble, avg_K::Float64) = -avg_K/e.β
abs_E_madelung(N::Int, λ::Float64) = 2.83729747948527 * pi/2.0 * N * λ
E_int_from_Hartree(E_Ha::Float64, λ::Float64) = E_Ha /(16/((2*pi)^4 * (λ/2)^2) * 0.5)
E_int_from_Hartree(E_Ha::Float64, e::Ensemble) = E_Ha /(16/((2*pi)^4 * (e.λ/2)^2) * 0.5)
E_Ry(E_internal::Float64, λ::Float64) = E_internal * 16/((2*pi)^4 * λ^2)
E_Ry(E_internal::Float64, e::Ensemble) = E_internal * 16/((2*pi)^4 * e.λ^2)
E_Ha(E_internal::Float64, λ::Float64) = E_internal * 16/((2*pi)^4 * λ^2) * 0.5
E_Ha(E_internal::Float64, e::Ensemble) = E_internal * 16/((2*pi)^4 * e.λ^2) * 0.5
Et_Ry(E_internal::Float64, e::Ensemble) = E_Ry(E_internal-abs_E_madelung(e.N, e.λ),e.λ)
Et_Ha(E_internal::Float64, e::Ensemble) = E_Ha(E_internal-abs_E_madelung(e.N, e.λ),e.λ)

#To run on a Linux System use "julia --threads NT run_threads.jl", where NT is the
#desired number of Threads.
#Inside the Code you can use Threads.nthreads() to check how many Threads
#the Programm uses.
function main()
    # MC options
    NMC = 10^5
    cyc = 50
    N_Runs = 24
    NEquil = 2*10^4
    # system parameters
    θ = 0.125
    rs = 2.0

    # use 7 particles
    S = sphere_with_same_spin(PlaneWave((0,0,0),Up),dk=1)
    #S = sphere(PlaneWave((0,0,0),Up),dk=2)
    N = length(S)
    c = Configuration(S)
    d = dimension(c.occupations)
    ξ = fractional_spin_polarization(S)

    println("parameters:")
    println("θ: ", θ)
    println("rs: ", rs)
    println("N: ", N)
    println("ξ: ", ξ)

    e = CEnsemble(λ(N,rs,d), β(θ,N,ξ,d), N)

    updates = [move_particle!, add_type_B!, remove_type_B!, add_type_C!, remove_type_C!, add_type_D!, remove_type_D!, add_type_E!, remove_type_E!, shuffle_indices!]

    measurements = Dict(
      :sign => (Variance(), signum)
    , :Ekin => (Variance(), Ekin)
    , :W_diag => (Variance(), W_diag)
    , :K => (Variance(), K)
    , :K_fermion => (Variance(), K)
    , :T1c => (Variance(), right_type_1_count)
    , :occs => (Group([Variance() for i in 1:100]), occupations)
    )

    measurements_Mean = Dict(
      :sign => (Mean(), signum)
    , :Ekin => (Mean(), Ekin)
    , :W_diag => (Mean(), W_diag)
    , :K => (Mean(), K)
    , :K_fermion => (Mean(), K)
    , :T1c => (Mean(), right_type_1_count)
    , :occs => (Group([Mean() for i in 1:100]), occupations)
    )

    println("starting mc process using $(Threads.nthreads()) threads... ")

    measurements_of_runs = Set{Dict{Symbol,Tuple{OnlineStat,Function}}}()

    total_rates = Dict()

    Threads.@threads for t in 1:N_Runs
        me = deepcopy(measurements_Mean)
        push!(measurements_of_runs,me)
        c_run = Configuration(copy(c.occupations))
            equilibrate_diagonal!(UEG(), e, c_run)
        rates = sweep!(UEG(), e, c_run, updates, me, NMC, cyc, NEquil)
        lock(ReentrantLock()) do
            merge!(+, total_rates, rates)
        end
    end
    #Addup counters

    println(" finished.")

    println("parameters:")
    println("θ: ", θ)
    println("rs: ", rs)
    println("N: ", N)

    println("measurements:")
    println("=============")

    # average over the uncorrelated mean values of the single runs
    for m in measurements_of_runs
        for (key,(stat,obs)) in m
            if key == :occs
                fit!(first(measurements[key]), mean.(m[:occs][1].stats))
            else
                fit!(first(measurements[key]), mean(stat))
            end
        end
    end

    # print measurements
    avg_sign = mean(first(measurements[:sign]))
    for (k,(f,m)) in measurements
        if typeof(f) == Variance{Float64,Float64,EqualWeight}
            if in(k,[:sign, :K, :T1c])
                println(k, "\t", mean(f), " +/- ", std(f)/sqrt(N_Runs-1))
            else
                println(k, "\t", mean(f)/avg_sign, " +/- ", std(f)/sqrt(N_Runs-1)/avg_sign)
            end
        end
    end

    # print addidtional observables
    μW_diag = mean(first(measurements[:W_diag]))/avg_sign
    ΔW_diag = std(first(measurements[:W_diag]))/sqrt(N_Runs-1)/avg_sign
    μW_off_diag = W_off_diag(e::Ensemble, mean(first(measurements[:K_fermion]))/avg_sign)
    ΔW_off_diag = abs(W_off_diag(e::Ensemble, std(first(measurements[:K_fermion]))/sqrt(N_Runs-1)/avg_sign))
    μT = mean(first(measurements[:Ekin]))/avg_sign
    ΔT = std(first(measurements[:Ekin]))/sqrt(N_Runs-1)/avg_sign
    μW = μW_diag + μW_off_diag
    ΔW = ΔW_diag + ΔW_off_diag
    μE = μW + μT
    ΔE = ΔW + ΔT
    μWt_Ha = Et_Ha(μW, e)
    ΔWt_Ha = E_Ha(ΔW, e)
    μT_Ha = E_Ha(μT,e.λ)
    ΔT_Ha = E_Ha(ΔT,e.λ)
    μE_Ha = μT_Ha + μWt_Ha
    ΔE_Ha = ΔT_Ha + ΔWt_Ha

    println("W_off_diag", "\t", μW_off_diag, " +/- ", ΔW_off_diag)
    println("W", "\t", μW, " +/- ", ΔW)
    println("E", "\t", μE, " +/- ", ΔE)
    println("W_t_Ha", "\t", μWt_Ha, " +/- ", ΔWt_Ha)
    println("T_Ha", "\t", μT_Ha, " +/- ", ΔT_Ha)

    println("")

    println("occupations:")
    println("============")
    println(mean.(measurements[:occs][1].stats))
    println(std.(measurements[:occs][1].stats))

    print_rates(total_rates)

    # create resultsfile
    # add measurements to file
    df = DataFrame(sign = 1)
    for (k,(f,m)) in measurements
        if typeof(f) == Variance{Float64,Float64,EqualWeight}
            if in(k,[:sign, :K])
                df[!,k] .= mean(f)
                df[!,Symbol(:Δ, k)] .= std(f)/sqrt(N_Runs-1)
            else
                df[!,k] .= mean(f)/avg_sign
                df[!,Symbol(:Δ, k)] .= std(f)/sqrt(N_Runs-1)/avg_sign
            end
        end
    end

    # add additional variables to file
    df[!,:W] .= μW
    df[!,:ΔW] .= ΔW
    df[!,:E] .= μE
    df[!,:ΔE] .= ΔE
    df[!,:Wt_Ha] .= μWt_Ha
    df[!,:ΔWt_Ha] .= ΔWt_Ha
    df[!,:T_Ha] .= μT_Ha
    df[!,:ΔT_Ha] .= ΔT_Ha
    df[!,:E_Ha] .= μE_Ha
    df[!,:ΔE_Ha] .= ΔE_Ha

    # create occupation numbers file
    open("examples/out/occNums_$(N)_th$(replace(string(θ),"." => ""))_rs$(replace(string(rs),"." => ""))_steps$((NMC*N_Runs/cyc)).dat", "w") do io
        writedlm(io, zip(mean.(measurements[:occs][1].stats), std.(measurements[:occs][1].stats)/(NMC*N_Runs/cyc)))
    end

    CSV.write("examples/out/results_N$(N)_th$(θ)_rs$(rs)_steps$((NMC*N_Runs/cyc)).csv",df)
end

main()


text:
"""
bio_seq : logic associated with dealing with proteins/translation and
other sequence functions unique to their biological interpretation.
In general our order of aa's is ARNDCQEGHILKMFPSTWYV- .
This is the "alphabetical" order and comes from cross-comptability with Gremlin-deived MRFs.
"""

module bio_seq

using Unicode

islower(s) = all(c->islowercase(c) | isspace(c), s);

nucs = "ACGT"
quick_rev = Dict{String, String}("A" => "T", "C" => "G", "G" => "C", "T" => "A")
codon_table = Dict{String, Char}("..." => '-', "TGT"=>'C',"GAC"=>'D',"TTC"=>'F',"TAG"=>'*',"GTG"=>'V',"CCT"=>'P', "GCT"=>'A',"GGC"=>'G',"CGG"=>'R',
									 "ATT"=>'I',"GAT"=>'D',"CAG"=>'Q',"ATG"=>'M',"CTC"=>'L',"TCT"=>'S',"CGT"=>'R',"ACG"=>'T',"AGA"=>'R',
									 "TGG"=>'W',"TCG"=>'S',"TTA"=>'L',"AGT"=>'S',"CGA"=>'R',"TGC"=>'C',"CAA"=>'Q',"TTG"=>'L',"AAT"=>'N',
									 "AAC"=>'N',"TAA"=>'*',"TAC"=>'Y',"CCA"=>'P',"ACT"=>'T',"TAT"=>'Y',"CAC"=>'H',"CCG"=>'P',"GCA"=>'A',
									 "GAA"=>'E',"ACA"=>'T',"GTA"=>'V',"AGG"=>'R',"AGC"=>'S',"TCA"=>'S',"GGT"=>'G',"GCC"=>'A',"TGA"=>'*',
									 "GTC"=>'V',"TTT"=>'F',"CAT"=>'H',"AAG"=>'K',"AAA"=>'K',"CTT"=>'L',"CTA"=>'L',"ATC"=>'I',"ACC"=>'T',
									 "TCC"=>'S',"CCC"=>'P',"GAG"=>'E',"GCG"=>'A',"CTG"=>'L',"CGC"=>'R',"GTT"=>'V',"GGA"=>'G',"ATA"=>'I',
									 "GGG"=>'G')
aa_to_num = Dict{Char, Int64}('A'=>1, 'R'=>2, 'N'=>3, 'D'=>4, 'C'=>5, 'Q'=>6,
									'E'=>7, 'G'=>8, 'H'=>9, 'I'=>10, 'L'=>11, 'K'=>12, 'M'=>13, 'F'=>14, 'P'=>15,
									'S'=>16, 'T'=>17, 'W'=>18, 'Y'=>19, 'V'=>20, '-'=>21, '*'=>21)

function load_fasta(file_name)
	in_file = open(file_name)
	in_read = readlines(in_file)
	close(in_file)

	seqs = Dict{String, String}()
	cur_head = ""
	cur_seq = String[]
	for line in in_read
		if !(isempty(line)) && line[1] == '>'
			if cur_head != ""
				seqs[cur_head] = join(cur_seq)
			end
			cur_head = strip(line[2:end])
			empty!(cur_seq)
		else
			push!(cur_seq, strip(line))
		end
	end
	seqs[cur_head] = join(cur_seq)
	return seqs
end

function convert_protein(prot_seq)
	prot_len = length(prot_seq)
	big_vector = zeros(Float64, prot_len * 21)
	#and fill it up like a very sparse row vector. One-hot.
	#Order is ARNDCQEGHILKMFPSTWYV-
	prot_fill = Dict{Char, Int64}('A'=>1, 'R'=>2, 'N'=>3, 'D'=>4, 'C'=>5, 'Q'=>6, 'E'=>7, 'G'=>8, 'H'=>9, 'I'=>10,
									'L'=>11, 'K'=>12, 'M'=>13, 'F'=>14, 'P'=>15, 'S'=>16, 'T'=>17, 'W'=>18, 'Y'=>19, 'V'=>20, '-'=>21, 'X'=>21, '*'=>21)
	chr_count = 0
	for chr in prot_seq
		big_vector[(21 * chr_count) + prot_fill[chr]] = 1.0
		chr_count += 1
	end
	return big_vector' #more convenient for it to be a row vector later.
end

function translate(seq, skip=1)
	#adding ellipses for dashes we can't change.
	seq = replace(seq, "*"=>"") #just skip the insertions which we won't change here.
	codons = [seq[x:min(x+2, length(seq))] for x in skip:3:length(seq)]
	out_seq = Char[]
	for codon in codons
		if length(codon) < 3
			continue
		end
		if islower(codon)
			push!(out_seq, lowercase(codon_table[uppercase(codon)]))
		else
			push!(out_seq, codon_table[uppercase(codon)])
		end
	end
	return join(out_seq)
end

function rev_comp(seq)
	return(join([quick_rev[string(B)] for B in reverse(seq)], ""))
end

function get_two_base()
	two_base = [join([one, two]) for one in nucs for two in nucs]
end

function get_three_base() #64.
	three_base = [join([one, two, three]) for one in nucs for two in nucs for three in nucs]
	return three_base
end

function get_four_base() #256 combos.
	four_base = [join([one, two, three, four]) for one in nucs for two in nucs for three in nucs for four in nucs]
	return four_base
end

function get_six_base()
	six_base = [join([one, two, three, four, five, six]) for one in nucs for two in nucs for three in nucs for four in nucs for five in nucs for six in nucs]
	return six_base
end

function precompute_dicodons() #64 * 64 codons = 4096 options in dictionary.
	#Let's precompute all the possible codon pairs to translate very quickly.
	all_bases = get_six_base()
	di_codon_dict = Dict{AbstractString, Tuple{Int64, Int64}}()
	for di_codon in all_bases
		di_codon_dict[di_codon] = translate_to_tuple(di_codon)
	end
	return di_codon_dict
end

function get_base_dict()
	mut_dict = Dict{Int64, Array{AbstractString, 1}}()
	mut_dict[2] = get_two_base()
	mut_dict[3] = get_three_base()
	mut_dict[4] = get_four_base()
	return mut_dict
end

function translate_to_tuple(seq)
	#adding ellipses for dashes we can't change.
	#seq = replace(seq, "*"=>"") #just skip the insertions which we won't change here.
	codons = [seq[x:min(x+2, length(seq))] for x in 1:3:length(seq)]
	out_seq = Char[]
	for codon in codons
		codon = uppercase(codon)
		if length(codon) < 3
			continue
		end
		push!(out_seq, codon_table[uppercase(codon)])
	end
	out_tuple::Tuple{Int64, Int64} = (aa_to_num[out_seq[1]], aa_to_num[out_seq[2]])
	return out_tuple
end

function translate_constrained_maybe_map(nuc, trns; do_map = false)
	gap_poses = Dict{Int64, Int64}()
	for (i, j) in nuc.gap_pos
		if i in keys(gap_poses)
			gap_poses[i] += j
		else
			gap_poses[i] = j
		end
	end

	true_nuc = Char[]
	i_count = 1
	explicit_mapping = Dict{Int64, Int64}()
	aa_counter = 0
	adder = 0

	for unmapped_n in 1:length(nuc.nucs)
		n = unmapped_n + (trns - 1)
		aux_n = n + adder
		if !(n in nuc.skip_trn)
			#then we ignore n for our purposes.
			explicit_mapping[aux_n] = div(aa_counter, 3) + 1
			aa_counter += 1
			push!(true_nuc, nuc.nucs[unmapped_n])
		end
		if n in keys(gap_poses)
			for ii in 1:gap_poses[n]
				for iii in 1:3
					aa_counter += 1
					push!(true_nuc, '.')
				end
			end
		end
	end

	ret_seq = translate(join(true_nuc), 1)
	gap_counter = 0
	if do_map
		return ret_seq, explicit_mapping
	else
		return ret_seq
	end
end

end


text:
using Distributions

module BayesChp
export bchp

function hazard_func(r, lambda)
    (1.0 / lambda) * ones(size(r))
end

function gauss_pdf(ms, val, sigma = 1)
  res = Float64[]
  for m in ms
    d = Normal(m, sigma)
    pdfgauss = pdf(d, val)
    push!(res, pdfgauss)
  end
  res
end
#println( gauss_pdf([1,2,3], 2.9) )

function estimate_mus(x)
  n = length(x)
  mus = Float64[]
  for i = 1:(n-1)
    subvec = x[(n-i) : end]
    push!(mus, mean(subvec))
  end
  #if n == 1
  #  push!(mus, x[1])
  #end
  mus
end
r = estimate_mus([1, 1, 5, 6])
print(r)

function bchp(X, mu0, mdist)
  T = length(X)
  R = zeros(T+1, T+1)
  R[1,1] = 1
  maxes = zeros(T+1)
  for t = 1 : T-1

    mus = estimate_mus(X[1:(t+1)])
    predprobs = gauss_pdf(mus, X[t], 1.0)

    H = hazard_func(collect(1:t), mdist)

    # Growth probabilities
    R[2:t+1, t + 1] = R[1:t, t] .* predprobs .* (1 - H)

    # Probability that there was a changepoint
    R[1, t+1] = sum( R[1:t, t] .* predprobs .* H )

    # Normalize
    R[:, t+1] = R[:, t+1] ./ sum(R[:, t+1])
    #
    maxes[t] = indmax(R[:, t])
  end
  R
  maxes
end

end
sig = readdlm("D:/MyTemp/gitlocal/recurrency-code/src.jl/sigtest.txt")

R = bchp(sig, 1.1, 100)
writedlm("D:/MyTemp/gitlocal/recurrency-code/src.jl/probs100.txt", R)

# R = bchp(sig, 1.1, 50)
# writedlm("D:/MyTemp/gitlocal/recurrency-code/src.jl/probs50.txt", R)
#
# R = bchp(sig, 1.1, 200)
# writedlm("D:/MyTemp/gitlocal/recurrency-code/src.jl/probs200.txt", R)

max_stars_repo_name: av-maslov/Recurrency

max_issues_repo_name: av-maslov/Recurrency

max_forks_repo_name: av-maslov/Recurrency


text:
using Test
using H3
using H3Geometry
using LibGEOS
using GeoInterface

# verify basics
@test H3Geometry.geoToH3(42, -110, 6) == H3.API.geoToH3(H3.Lib.GeoCoord(deg2rad(42.0), deg2rad(-110.0)), 6)
@test H3Geometry.geoToH3([42.0, -110], 6) == H3.API.geoToH3(H3.Lib.GeoCoord(deg2rad(42.0), deg2rad(-110.0)), 6)
@test H3Geometry.geoToH3(LibGEOS.Point(42.0, -110), 6) == H3.API.geoToH3(H3.Lib.GeoCoord(deg2rad(42.0), deg2rad(-110.0)), 6)

# verify that H3.API.geoToH3 is overloaded
@test H3.API.geoToH3(42, -110, 6) == H3.API.geoToH3(H3.Lib.GeoCoord(deg2rad(42.0), deg2rad(-110.0)), 6)

# validate that reasonable dot usage works
@test begin
    @info "Convert combined with broadcast"
    p = coordinates(boundary(LibGEOS.Polygon([[[0.0, 0.0], [1, 1], [2, 0], [0, 0]]])))
    H3.API.geoToH3.(p, 5) == [0x085754e67fffffff, 0x0857541affffffff, 0x0857542b7fffffff, 0x085754e67fffffff]
end

# verify that a hexagon appears in roughly the right box and has the right shape
@test begin
    @info "Hexagonal sanity check"
    p0 = H3.API.geoToH3(46, -110,5)
    p1 = H3Geometry.h3ToPolygon(p0)
    p2 = LibGEOS.Polygon([[[45.85, -110.15], [46.1, -110.15], [46.1, -109.85], [45.85, -109.85], [45.85, -110.15]]])
    length(coordinates(boundary(p1))) == 7 && area(difference(p1, p2)) == 0.0 && area(p1) / area(p2) > 0.1
end

# verify that we can do unions of H3index values
@test begin
    @info "Unions of H3index"
    p0 = H3.API.geoToH3(46, -110, 5)
    p1 = H3Geometry.h3ToPolygon(p0)
    p2 = H3.API.kRing(p0, 1)
    p3 = H3Geometry.h3ToPolygon(p2)
    p4 = H3Geometry.h3ToPolygon(H3.API.geoToH3(46, -110, 4))
    isapprox(area(p3), 7 * area(p1), atol=1e-4) && isapprox(area(p3), area(p4), atol=1e-3)
end

# verify minimal overspray for polyfill with cover=0
@test begin
    @info "Polyfill overspray"
    p0 = LibGEOS.Polygon([[[45.5, -110.2], [46.1, -110.2], [46.1, -109.85], [45.5, -109.85], [45.5, -110.2]]])
    p1 = H3Geometry.polyfill(p0, 6, cover=0)
    ratio =  area(H3Geometry.h3ToPolygon(p1)) / area(p0)
    ratio < 1 && ratio >= 0.99
end

# verify minimal underspray for polyfill with cover=1.2
@test begin
    @info "Polyfill underspray"
    p0 = LibGEOS.Polygon([[[45.5, -110.2], [46.1, -110.2], [46.1, -109.85], [45.5, -109.85], [45.5, -110.2]]])
    p1 = H3Geometry.h3ToPolygon(H3Geometry.polyfill(p0, 7, cover=1.2))
    ratio =  area(p1) / area(p0)
    a1 = area(difference(p1, p0))
    a2 = area(difference(p0, p1))
    ratio > 1.1 && a1 < 0.15 && a2 < 1e-9
end


text:
qml_file  = path*"/SampleFiles/XML/ISC_2011-tohoku-oki.xml"
uw_file   = path*"/SampleFiles/UW/99011116541"
suds_file = path*"/SampleFiles/SUDS/eq_wvm1.sud"

printstyled("  read_quake wrapper\n", color=:light_green)

printstyled("    QML\n", color=:light_green)
H, R = read_qml(qml_file)
Ev1 = SeisEvent(hdr = H[1], source = R[1])
Ev2 = read_quake("qml", qml_file)
@test Ev2.hdr.src == abspath(qml_file)
@test Ev2.source.src == abspath(qml_file)
Ev2.hdr.src = Ev1.hdr.src
Ev2.source.src = Ev1.source.src
@test Ev1 == Ev2

printstyled("    SUDS\n", color=:light_green)
Ev1 = readsudsevt(suds_file)
Ev2 = read_quake("suds", suds_file)
@test Ev1 == Ev2

printstyled("    UW\n", color=:light_green)
Ev1 = readuwevt(uw_file, full=true)
Ev2 = read_quake("uw", uw_file, full=true)
@test Ev1 == Ev2


text:
using WiltonInts84
using FixedSizeArrays
using Base.Test

p1 = Vec(0.0, 0.0, 0.0)
p2 = Vec(1.0, 0.0, 0.0)
p3 = Vec(0.0, 1.0, 0.0)
h = 0.1

x = -0.5 : 0.01 : 1.5
y = -0.5 : 0.01 : 1.5

r, R = 0.0, 0.15
fails = 0
M = zeros(typeof(h), length(x), length(y))
for i in eachindex(x)
    for j in eachindex(y)
        c = Vec(x[i],y[j],h)
        try
            I, K = wiltonints(p1,p2,p3,c,r,R,Val{0})
            M[i,j] = I[2]
        catch
            @show i j c
            fails += 1
        end
    end
end

# There are still some cases where the package fails,
# but at least lets make sure that number does not increase...
@test fails == 20


text:
################## To IMatrix ######################
IMatrix{N,T}(::AbstractMatrix) where {N,T} = IMatrix{N,T}()
IMatrix{N}(A::AbstractMatrix{T}) where {N,T} = IMatrix{N,T}()
IMatrix(A::AbstractMatrix{T}) where {T} =
    IMatrix{size(A, 1) == size(A, 2) ? size(A, 2) : throw(DimensionMismatch()),T}()

################## To Diagonal ######################
Diagonal(A::PermMatrix) = Diagonal(diag(A))
Diagonal(A::IMatrix{N,T}) where {N,T} = Diagonal{T}(ones(T, N))
Diagonal{T}(::IMatrix{N}) where {T,N} = Diagonal{T}(ones(T, N))

################## To SparseMatrixCSC ######################
SparseMatrixCSC{Tv,Ti}(A::IMatrix{N}) where {Tv,Ti<:Integer,N} =
    SparseMatrixCSC{Tv,Ti}(I, N, N)
SparseMatrixCSC{Tv}(A::IMatrix) where {Tv} = SparseMatrixCSC{Tv,Int}(A)
SparseMatrixCSC(A::IMatrix{N,T}) where {N,T} = SparseMatrixCSC{T,Int}(I, N, N)
function SparseMatrixCSC(M::PermMatrix)
    n = size(M, 1)
    #SparseMatrixCSC(n, n, collect(1:n+1), M.perm, M.vals)
    order = invperm(M.perm)
    SparseMatrixCSC(n, n, collect(1:n+1), order, M.vals[order])
end

@static if VERSION < v"1.3-"

    function SparseMatrixCSC(D::Diagonal{T}) where {T}
        m = length(D.diag)
        return SparseMatrixCSC(m, m, Vector(1:(m+1)), Vector(1:m), Vector{T}(D.diag))
    end

end

SparseMatrixCSC{Tv,Ti}(M::PermMatrix{Tv,Ti}) where {Tv,Ti} = SparseMatrixCSC(M)
SparseMatrixCSC(coo::SparseMatrixCOO) = sparse(coo.is, coo.js, coo.vs, coo.m, coo.n)

################## To Dense ######################
Matrix{T}(::IMatrix{N}) where {T,N} = Matrix{T}(I, N, N)
Matrix(::IMatrix{N,T}) where {N,T} = Matrix{T}(I, N, N)

function Matrix{T}(X::PermMatrix) where {T}
    n = size(X, 1)
    Mf = zeros(T, n, n)
    @simd for i = 1:n
        @inbounds Mf[i, X.perm[i]] = X.vals[i]
    end
    return Mf
end
Matrix(X::PermMatrix{T}) where {T} = Matrix{T}(X)

function Matrix(coo::SparseMatrixCOO{T}) where {T}
    mat = zeros(T, coo.m, coo.n)
    for (i, j, v) in zip(coo.is, coo.js, coo.vs)
        mat[i, j] += v
    end
    mat
end

################## To PermMatrix ######################
PermMatrix{Tv,Ti}(::IMatrix{N}) where {Tv,Ti,N} =
    PermMatrix{Tv,Ti}(Vector{Ti}(1:N), ones(Tv, N))
PermMatrix{Tv}(X::IMatrix) where {Tv} = PermMatrix{Tv,Int}(X)
PermMatrix(X::IMatrix{N,T}) where {N,T} = PermMatrix{T,Int}(X)

PermMatrix{Tv,Ti}(A::PermMatrix) where {Tv,Ti} =
    PermMatrix(Vector{Ti}(A.perm), Vector{Tv}(A.vals))

function _findnz(A::AbstractMatrix)
    I = findall(!iszero, A)
    getindex.(I, 1), getindex.(I, 2), A[I]
end

_findnz(A::AbstractSparseArray) = findnz(A)

function PermMatrix{Tv,Ti}(A::AbstractMatrix) where {Tv,Ti}
    i, j, v = _findnz(A)
    j == collect(1:size(A, 2)) || throw(ArgumentError("This is not a PermMatrix"))
    order = invperm(i)
    PermMatrix{Tv,Ti}(Vector{Ti}(order), Vector{Tv}(v[order]))
end
PermMatrix(A::AbstractMatrix{T}) where {T} = PermMatrix{T,Int}(A)
PermMatrix(A::SparseMatrixCSC{Tv,Ti}) where {Tv,Ti} = PermMatrix{Tv,Ti}(A) # inherit indice type
PermMatrix{Tv,Ti}(A::Diagonal{Tv}) where {Tv,Ti} =
    PermMatrix(Vector{Ti}(1:size(A, 1)), A.diag)
#PermMatrix(A::Diagonal{T}) where T = PermMatrix{T, Int}(A)
# lazy implementation
function PermMatrix{Tv,Ti,Vv,Vi}(
    A::AbstractMatrix,
) where {Tv,Ti<:Integer,Vv<:AbstractVector{Tv},Vi<:AbstractVector{Ti}}
    pm = PermMatrix(PermMatrix{Tv,Ti}(A))
    PermMatrix(Vi(pm.perm), Vv(pm.vals))
end

############## To SparseMatrixCOO ##############
function SparseMatrixCOO(A::Matrix{Tv}; atol = 1e-12) where {Tv}
    m, n = size(A)
    is = Int[]
    js = Int[]
    vs = Tv[]
    for j = 1:n
        for i = 1:m
            if abs(A[i, j]) > atol
                push!(is, i)
                push!(js, j)
                push!(vs, A[i, j])
            end
        end
    end
    SparseMatrixCOO(is, js, vs, m, n)
end

import Base: convert
convert(T::Type{<:PermMatrix}, m::AbstractMatrix) = m isa T ? m : T(m)
convert(T::Type{<:IMatrix}, m::AbstractMatrix) = m isa T ? m : T(m)
convert(T::Type{<:Diagonal}, m::AbstractMatrix) = m isa T ? m : T(m)


text:
"""
    JuMP.set_objective_function(model::InfiniteModel,
                                func::JuMP.AbstractJuMPScalar)

Extend [`JuMP.set_objective_function`](@ref) to set the objective expression of
infinite model `model`. Errors if `func` contains infinite variables and/or
parameters. Also errors if `func` contains invalid variables.

**Example**
```jldoctest; setup = :(using InfiniteOpt, JuMP; model = InfiniteModel(); @hold_variable(model, x))
julia> set_objective_function(model, 2x + 1)

julia> objective_function(model)
2 x + 1
```
"""
function JuMP.set_objective_function(model::InfiniteModel,
                                     func::JuMP.AbstractJuMPScalar)
    # check the function
    if isa(func, InfiniteExpr) || isa(func, ParameterExpr)
        error("Objective function cannot contain infinite parameters/variables.")
    end
    JuMP.check_belongs_to_model(func, model)
    # update the function
    model.objective_function = func
    # delete old mappings
    for vindex in keys(model.var_in_objective)
        model.var_in_objective[vindex] = false
    end
    for mindex in keys(model.meas_in_objective)
        model.meas_in_objective[mindex] = false
    end
    # update new mappings
    vrefs = _all_function_variables(func)
    for vref in vrefs
        if isa(vref, InfOptVariableRef)
            model.var_in_objective[JuMP.index(vref)] = true
        elseif isa(vref, MeasureRef)
            model.meas_in_objective[JuMP.index(vref)] = true
        end
    end
    set_optimizer_model_ready(model, false)
    return
end

"""
    JuMP.set_objective_function(model::InfiniteModel, func::Real)

Extend [`JuMP.set_objective_function`](@ref) to set the objective expression of
`model` with a number.

**Example**
```jldoctest; setup = :(using InfiniteOpt, JuMP; model = InfiniteModel())
julia> set_objective_function(model, 3)

julia> objective_function(model)
3
```
"""
function JuMP.set_objective_function(model::InfiniteModel, func::Real)
    # update function
    model.objective_function = JuMP.GenericAffExpr{Float64, HoldVariableRef}(func)
    # delete old mappings
    for vindex in keys(model.var_in_objective)
        model.var_in_objective[vindex] = false
    end
    for mindex in keys(model.meas_in_objective)
        model.meas_in_objective[mindex] = false
    end
    set_optimizer_model_ready(model, false)
    return
end

"""
    JuMP.set_objective_sense(model::InfiniteModel, sense::MOI.OptimizationSense)

Extend [`JuMP.set_objective_sense`](@ref JuMP.set_objective_sense(::JuMP.Model, ::MOI.OptimizationSense))
to set the objective sense of infinite model `model`.

**Example**
```jldoctest; setup = :(using InfiniteOpt, JuMP, MathOptInterface; model = InfiniteModel(), const MOI = MathOptInterface)
julia> set_objective_sense(model, MOI.MIN_SENSE)

julia> objective_sense(model)
MIN_SENSE::OptimizationSense = 0
```
"""
function JuMP.set_objective_sense(model::InfiniteModel,
                                  sense::MOI.OptimizationSense)
    model.objective_sense = sense
    set_optimizer_model_ready(model, false)
    return
end

"""
    JuMP.set_objective(model::InfiniteModel, sense::MOI.OptimizationSense,
                       func::Union{JuMP.AbstractJuMPScalar, Real})

Extend `JuMP.set_objective` to set the objective of infinite model
`model`. Errors if `func` contains infinite variables and/or parameters, or if
it does not belong to the model.

**Example**
```jldoctest; setup = :(using InfiniteOpt, JuMP, MathOptInterface; model = InfiniteModel(), const MOI = MathOptInterface; @hold_variable(model, x))
julia> set_objective(model, MOI.MIN_SENSE, 2x + 1)

julia> objective_function(model)
2 x + 1
```
"""
function JuMP.set_objective(model::InfiniteModel, sense::MOI.OptimizationSense,
                            func::Union{JuMP.AbstractJuMPScalar, Real})
    JuMP.set_objective_sense(model, sense)
    JuMP.set_objective_function(model, func)
    return
end

# Fallback
function JuMP.set_objective(model::InfiniteModel, sense::MOI.OptimizationSense,
                            func)
    error("The objective function `$(func)` is not supported.")
end

"""
    JuMP.objective_sense(model::InfiniteModel)::MOI.OptimizationSense

Extend [`JuMP.objective_sense`](@ref JuMP.objective_sense(::JuMP.Model)) to
return the objective sense of the infinite model `model`.

**Example**
```jldoctest; setup = :(using InfiniteOpt, JuMP; model = InfiniteModel(), @objective(model, Min, 1))
julia> objective_sense(model)
MIN_SENSE::OptimizationSense = 0
```
"""
function JuMP.objective_sense(model::InfiniteModel)::MOI.OptimizationSense
    return model.objective_sense
end

"""
    JuMP.objective_function_type(model::InfiniteModel)::Type{<:JuMP.AbstractJuMPScalar}

Extend [`JuMP.objective_function_type`](@ref JuMP.objective_function_type(::JuMP.Model))
to return the objective function type of infinite model `model`.

**Example**
```jldoctest; setup = :(using InfiniteOpt, JuMP; model = InfiniteModel(), @objective(model, Min, 1))
julia> objective_function_type(model)
GenericAffExpr{Float64,HoldVariableRef}
```
"""
function JuMP.objective_function_type(model::InfiniteModel)::Type{<:JuMP.AbstractJuMPScalar}
    return typeof(model.objective_function)
end

"""
    JuMP.objective_function(model::InfiniteModel)::JuMP.AbstractJuMPScalar

Extend [`JuMP.objective_function`](@ref JuMP.objective_function(::JuMP.Model))
to return the objective of infinite model `model`.

**Example**
```jldoctest; setup = :(using InfiniteOpt, JuMP; model = InfiniteModel(), @objective(model, Min, 1))
julia> objective_function(model)
1
```
"""
function JuMP.objective_function(model::InfiniteModel)::JuMP.AbstractJuMPScalar
    return model.objective_function
end

# Extend objective_function for infinite models
function JuMP.objective_function(model::InfiniteModel, FT::Type)
    # InexactError should be thrown, this is needed in `objective.jl`
    if !(model.objective_function isa FT)
        throw(InexactError(:objective_function, FT,
                           typeof(model.objective_function)))
    end
    return model.objective_function::FT
end

"""
    JuMP.set_objective_coefficient(model::InfiniteModel,
                                   variable::GeneralVariableRef,
                                   coefficient::Real)

Extend [`JuMP.set_objective_coefficient`](@ref JuMP.set_objective_coefficient(::JuMP.Model, ::JuMP.VariableRef, ::Real))
Set the linear objective coefficient associated with `variable` to `coefficient`.
Errors if the function type is unsupported.

**Example**
```jldoctest; setup = :(using InfiniteOpt, JuMP; model = InfiniteModel())
julia> @hold_variable(model, x)
x

julia> @hold_variable(model, y)
y

julia> @objective(model, x + y)
x + y

julia> set_objective_coefficient(model, y, 2)

julia> objective_function(model)
x + 2 y
```
"""
function JuMP.set_objective_coefficient(model::InfiniteModel,
                                        variable::GeneralVariableRef,
                                        coeff::Real)
    new_expr = _set_variable_coefficient!(JuMP.objective_function(model),
                                          variable, coeff)
    JuMP.set_objective_function(model, new_expr)
    return
end


text:
using DICOMClient

using HTTP
using Test

@testset "version" begin
    @test DICOMClient.package_version() isa VersionNumber
    @test DICOMClient.package_version() > v"0"
end

@testset "_add_trailing_slash" begin
    @test DICOMClient._add_trailing_slash(HTTP.URI("https://juliahealth.org")) == HTTP.URI("https://juliahealth.org/")
    @test DICOMClient._add_trailing_slash(HTTP.URI("https://juliahealth.org/")) == HTTP.URI("https://juliahealth.org/")
end

@testset "parse_multipart_image" begin
    @test_throws ArgumentError DICOMClient.parse_multipart_image(Val(:jpeg), UInt8[0xff, 0xd8, 0xff, 0xd9, 0xff, 0xd8])
end


text:
x = 0
for i = 0:999
    for j = 0:999
        for k = 0:99
            x += i * j * k
        end
    end
end
print(x)

max_stars_repo_name: vagoff/shootout

max_issues_repo_name: vagoff/shootout

max_forks_repo_name: vagoff/shootout


text:
using CryptoGroups: @bin_str, @hex_str, EC2N, PB, specialize, ECPoint, <|, generator, octet
using CryptoGroups.Specs: octet2int
# Elliptic Curve Domain Parameter Setup

basis = PB(hex"80000000 00000000 00000000 00000000 00000000 00000201", 191)

curve_spec = EC2N(basis;   ### Need to think about adding a proper methods
                  a = hex"2866537B 67675263 6A68F565 54E12640 276B649E F7526267",
                  b = hex"2E45EF57 1F00786F 67B0081B 9495A3D9 5462F5DE 0AA185EC",
                  G = hex"04 36B3DAF8 A23206F9 C4F299D7 B21A9C36 9137F2C8 4AE1AA0D 765BE734 33B3F95E 332932E7 0EA245CA 2418EA0E F98018FB",
                  n = 1569275433846670190958947355803350458831205595451630533029,
                  h = 2
)

P = specialize(ECPoint, curve_spec)
G = P <| generator(curve_spec)

# Key Generation

d = 1275552191113212300012030439187146164646146646466749494799
Q = d*G

@test octet(Q) == hex"04 5DE37E75 6BD55D72 E3768CB3 96FFEB96 2614DEA4 CE28A2E7 55C0E0E0 2F5FB132 CAF416EF 85B229BB B8E13520 03125BA1"

# Signature Generation

M = "abc"
H = Hash("SHA1")
e = octet2int(H(M))

## Elliptic curve computaion

k = 1542725565216523985789236956265265265235675811949404040041
R = k*G

x = octet(gx(R)) # Instead I could have octet(gx(R))
y = octet(gy(R))

@test x == hex"438E5A11 FB55E4C6 5471DCD4 9E266142 A3BDF2BF 9D5772D5"
@test y == hex"2AD603A0 5BD1D177 649F9167 E6F475B7 E2FF590C 85AF15DA"

x̄ = octet2int(x)
ȳ = octet2int(y)

@test x̄ == 1656469817011541734314669640730254878828443186986697061077

n = order(P)
r = x̄ % n

@test r == 87194383164871543355722284926904419997237591535066528048

## Modular computation

s = modinv(k, n) * (e + d*r) % n

@test s == 308992691965804947361541664549085895292153777025772063598

println("""The signature is:
    r = $r
    s = $s
""")

# Signature verification

(r′, s′) = r, s

M′ = "abc"
H′ = Hash("sha256")
e′ = octet2int(H(M′))

@test 1 < r′ < n - 1
@test 1 < s′ < n - 1

c = modinv(s′, n)

@test c == 952933666850866331568782284754801289889992082635386177703

u₁ = e′*c % n
u₂ = r′*c % n

@test u₁ == 1248886407154707854022434516084062503301792374360994400066
@test u₂ == 527017380977534012168222466016199849611971141652753464154

W = u₁*G + u₂*Q
x₁ = octet(gx(W))
y₁ = octet(gy(W))

@test x₁ == hex"438E5A11 FB55E4C6 5471DCD4 9E266142 A3BDF2BF 9D5772D5"
@test y₁ == hex"2AD603A0 5BD1D177 649F9167 E6F475B7 E2FF590C 85AF15DA"

x̄₁ = octet2int(x₁)
ν = x̄₁ % n

@test ν == r′


text:
function GeoInterface.geotype(g::AbstractGeometry)
    gtype = getgeomtype(g)
    if gtype in (GDAL.wkbPoint, GDAL.wkbPoint25D, GDAL.wkbPointM, GDAL.wkbPointZM)
        return :Point
    elseif gtype in (GDAL.wkbMultiPoint, GDAL.wkbMultiPoint25D, GDAL.wkbMultiPointM, GDAL.wkbMultiPointZM)
        return :MultiPoint
    elseif gtype in (GDAL.wkbLineString, GDAL.wkbLineString25D, GDAL.wkbLineStringM, GDAL.wkbLineStringZM)
        return :LineString
    elseif gtype == GDAL.wkbLinearRing
        return :LinearRing
    elseif gtype in (GDAL.wkbMultiLineString, GDAL.wkbMultiLineString25D, GDAL.wkbMultiLineStringM, GDAL.wkbMultiLineStringZM)
        return :MultiLineString
    elseif gtype in (GDAL.wkbPolygon, GDAL.wkbPolygon25D, GDAL.wkbPolygonM, GDAL.wkbPolygonZM)
        return :Polygon
    elseif gtype in (GDAL.wkbMultiPolygon, GDAL.wkbMultiPolygon25D, GDAL.wkbMultiPolygonM, GDAL.wkbMultiPolygonZM)
        return :MultiPolygon
    elseif gtype in (GDAL.wkbGeometryCollection, GDAL.wkbGeometryCollection25D, GDAL.wkbGeometryCollectionM, GDAL.wkbGeometryCollectionZM)
        return :GeometryCollection
    else
        @warn "unknown geometry type" gtype
        return :Unknown
    end
end

function GeoInterface.coordinates(g::AbstractGeometry)
    gtype = getgeomtype(g)
    ndim = getcoorddim(g)
    if gtype in (GDAL.wkbPoint, GDAL.wkbPoint25D, GDAL.wkbPointM, GDAL.wkbPointZM)
        if ndim == 2
            return Float64[getx(g,0), gety(g,0)]
        elseif ndim == 3
            return Float64[getx(g,0), gety(g,0), getz(g,0)]
        else
            @assert ndim == 0
            @warn("Empty Point")
        end
    elseif gtype in (GDAL.wkbMultiPoint, GDAL.wkbMultiPoint25D, GDAL.wkbMultiPointM, GDAL.wkbMultiPointZM)
        return Vector{Float64}[
            GeoInterface.coordinates(getgeom(g,i-1)) for i in 1:ngeom(g)
        ]
    elseif gtype in (GDAL.wkbLineString, GDAL.wkbLineString25D, GDAL.wkbLineStringM, GDAL.wkbLineStringZM) || gtype == GDAL.wkbLinearRing
        return Vector{Float64}[
            collect(getpoint(g,i-1)[1:ndim]) for i in 1:ngeom(g)
        ]
    elseif gtype in (
            GDAL.wkbMultiLineString, GDAL.wkbMultiLineString25D, GDAL.wkbMultiLineStringM, GDAL.wkbMultiLineStringZM,
            GDAL.wkbPolygon, GDAL.wkbPolygon25D, GDAL.wkbPolygonM, GDAL.wkbPolygonZM
        )
        return Vector{Vector{Float64}}[
            GeoInterface.coordinates(getgeom(g,i-1)) for i in 1:ngeom(g)
        ]
    elseif gtype in (GDAL.wkbMultiPolygon, GDAL.wkbMultiPolygon25D, GDAL.wkbMultiPolygonM, GDAL.wkbMultiPolygonZM)
        return Vector{Vector{Vector{Float64}}}[
            GeoInterface.coordinates(getgeom(g,i-1)) for i in 1:ngeom(g)
        ]
    end
end


text:
function ABCDQR_simula(nt::Int)
	#
	# simulated model to test the functions
	#
	# nt: number of time instants
	# no: number of outputs, no in [1,2]
	# ni: number of inputs, ni in [1,2]

	# random seed for the simulations
	srand(99)

	# simulated mechanical system
	#
	# |    k1    ------     k2    ------
	# |---###---|  m1  |---###---|  m2  |
	# |    c1    ------     c2    ------
	# |           @  @             @  @
	# |========================================

	# vibration parameters
	# ----------------------------
	# nyquist frequency (Hz)
	fnq = 8
	# sampling frequency (Hz)
	fs = 2*fnq
	# time interval (s)
	dt = 1/fs

	# matrices for the mechanical system
	# ----------------------------------------
	# mass matrix
	m1 = 35
	m2 = 17.5
	M = [m1 0;0 m2]
	Minv = [1/m1 0;0 1/m2] # inv(M)

	# stiffness matrix
	k1 = 8750
	k2 = 3500
	K = [k1+k2 -k2;-k2 k2]
	# theoretical natural frecuencies: 1.71 Hz, 3.31 Hz

	# eigenvalues and eigenectors
	D,V = eig(K,M)

	# modal matrices
	Mm = V'*M*V
	Km = V'*K*V

	# damping ratios
	zm = [0.02 0.02]

	# modal damping matrix
	Gm = zeros(2,2)
	for i in 1:2
		Gm[i,i] = 2*sqrt(Mm[i,i]*Km[i,i])*zm[i]
	end

	# damping matrix
	G = (V')\Gm/V

	# matrices for the discrete state space model
	# --------------------------------------------------------------
	# Matrix A
	Ac = [zeros(2,2) eye(2);-Minv*K -Minv*G] # Continuous Matrix A
	A = expm(Ac*dt) # discrete Matrix A

	# Matrix C (continuous = discrete)
	C = [-Minv*K -Minv*G]

	# Matrix B
	# Js is the input location matrix
	J = eye(2)
	Bc = [zeros(2,2);Minv*J] # continuous matrix B
	B = (A-eye(4,4))*inv(Ac)*Bc # discrete matrix B

	# matrix D (continuous = discrete)
	D = Minv*J

	# input definition
	# ---------------------------------------
	# system input u_{t} = a_{t} + a_{t-1} (MA1)
	u = 100*randn(2,nt) # if 100 is not used, Q and R are too little

	# state equation, no noise
	# ---------------------------------
	x = zeros(4,nt)
	for t=1:nt-1
		x[:,t+1] = A*x[:,t] + B*u[:,t]
	end
	P = cov(x')
	Q = 0.01*P # std approx 0.1
	U = chol(Q)
	w = U'*randn(4,nt)

	# state equation with noise
	# ------------------------------------------
	x = zeros(4,nt)
	for t = 1:nt-1
		x[:,t+1] = A*x[:,t] + B*u[:,t] + w[:,t]
	end

	# observation equation
	# -------------------------------------
	y = C*x + D*u
	L = cov(y')
	R = 0.01*L # std approx 0.1
	U = chol(R)
	v = U'*randn(2,nt)
	y = y + v

	#
	S = cov(w',v')

	return Dict("y"=>y,"u"=>u,"A"=>A,"B"=>B,"C"=>C,"D"=>D,"Q"=>Q,"R"=>R,"S"=>S)

end


text:
export DoorKey

using Random

mutable struct DoorKey{W<:GridWorldBase, R} <: AbstractGridWorld
    world::W
    agent_pos::CartesianIndex{2}
    agent::Agent
    goal_reward::Float64
    reward::Float64
    rng::R
end

function DoorKey(;n = 7, agent_start_pos = CartesianIndex(2,2), agent_start_dir = RIGHT, goal_pos = CartesianIndex(n-1, n-1), rng = Random.GLOBAL_RNG)
    door = Door(:yellow)
    key = Key(:yellow)
    objects = (EMPTY, WALL, GOAL, door, key)
    world = GridWorldBase(objects, n, n)

    world[WALL, [1,n], 1:n] .= true
    world[WALL, 1:n, [1,n]] .= true
    world[GOAL, goal_pos] = true

    goal_reward = 1.0
    reward = 0.0

    env = DoorKey(world, agent_start_pos, Agent(;dir = agent_start_dir), goal_reward, reward, rng)

    reset!(env, agent_start_pos = agent_start_pos, agent_start_dir = agent_start_dir, goal_pos = goal_pos)

    return env
end

function (env::DoorKey)(::MoveForward)
    env.reward = 0.0
    door = env.world.objects[end - 1]
    key = env.world.objects[end]
    dir = get_dir(env.agent)
    dest = dir(env.agent_pos)

    if env.world[key, dest]
        if PICK_UP(env.agent, key)
            env.world[key, dest] = false
            env.world[EMPTY, dest] = true
        end
        env.agent_pos = dest
    elseif env.world[door, dest] && env.agent.inventory !== key
        nothing
    elseif env.world[door, dest] && env.agent.inventory === key
        env.agent_pos = dest
    elseif !env.world[WALL,dest]
        env.agent_pos = dest
        if env.world[GOAL, env.agent_pos]
            env.reward = env.goal_reward
        end
    end

    env
end

RLBase.get_terminal(env::DoorKey) = env.world[GOAL, env.agent_pos]

function RLBase.reset!(env::DoorKey; agent_start_pos = CartesianIndex(2, 2), agent_start_dir = RIGHT, goal_pos = CartesianIndex(size(env.world)[end] - 1, size(env.world)[end] - 1))
    n = size(env.world)[end]
    door = env.world.objects[end - 1]
    key = env.world.objects[end]
    env.reward = 0.0
    env.agent_pos = agent_start_pos
    agent = get_agent(env)
    set_dir!(agent, agent_start_dir)

    door_pos = CartesianIndex(rand(env.rng, 2:n-1), rand(env.rng, 3:n-2))
    @assert agent_start_pos[2] < door_pos[2] "Agent should start on the left side of the door"
    @assert goal_pos[2] > door_pos[2] "Goal should be placed on the right side of the door"
    env.world[WALL, :, door_pos[2]] .= true
    env.world[door, door_pos] = true
    env.world[WALL, door_pos] = false

    key_pos = CartesianIndex(rand(env.rng, 2:n-1), rand(env.rng, 2:door_pos[2]-1))
    while key_pos == agent_start_pos
        key_pos = CartesianIndex(rand(env.rng, 2:n-1), rand(env.rng, 2:door_pos[2]-1))
    end
    env.world[key, key_pos] = true

    env.world[EMPTY, :, :] .= .!(.|((env.world[x, :, :] for x in [WALL, GOAL, door, key])...))
end


text:
using Yao, YaoExtensions
using Test, YaoArrayRegister

@test_throws ErrorException YaoExtensions.PauliString([X, Y, H])

g = YaoExtensions.PauliString(X, Y, Z)
@test collect(subblocks(g)) == [X, Y, Z]
@test g[1] == X
@test g[end] == Z

@test collect(g) == [X, Y, Z]
@test eltype(g) == ConstGate.PauliGate

r = rand_state(3)
@test apply!(copy(r), g) ≈ apply!(copy(r), kron(X, Y, Z))

g = YaoExtensions.PauliString([X, Y, Z])

@test mat(g) ≈ mat(kron(3, X, Y, Z))
@test chsubblocks(g, [X, X, X]) == YaoExtensions.PauliString(X, X, X)

g[3] = I2
@test occupied_locs(g) == (1, 2)

@test ishermitian(g) == ishermitian(mat(g))
@test isreflexive(g) == isreflexive(mat(g))
@test isunitary(g) == isunitary(mat(g))
@test length(g) == 3

@test cache_key(g) == [cache_key(g[1]), cache_key(g[2]), cache_key(g[3])]


text:
@testset "invalid" begin
    c = lxc(raw"""\begin{:sectiom}""")
    @test_throws ErrorException lx_begin(c,0)
end

@testset "section" begin
    counter = F.globvar("section_counter")
    c = lxc(raw"""\begin{:section, title="aaa"}""")
    h = lx_begin(c,0)
    @test F.globvar("section_counter") == counter + 1
    @test isapproxstr(strip(h), """
        ~~~

                  aaa

        ~~~""")
end


text:
##################
# Script         #
##################

include("Create_Master_Matrices.jl")

# This code can be run if you have .coords files in your directory
grid_coordinates("testfiles", 4, 4, (2,0))  #apply 4x4 grid cells to the coordinates before calculating GD
create_master_matrices("testfiles", "gridcells")

#apply 10-degree wide latitudinal bands to the coordinates
latband_coordinates("testfiles", 10)    # apply 10-degree wide latitudinal bands before calculating GD
create_master_matrices("testfiles", "latbands")

# and so forth if you have data for biomes, equal-area-cells etc
#create_master_matrices("testfiles", "anthromes")
#create_master_matrices("testfiles", "biomes")
#create_master_matrices("testfiles", "equal_area")

max_stars_repo_path: Miraldo et al. 2016/Create_matrices_script.jl

max_stars_repo_name: mkborregaard/codes-for-published-papers

max_issues_repo_name: DalavanCloud/codes-for-published-papers

max_forks_repo_name: DalavanCloud/codes-for-published-papers


text:
# using nctools
# using Test

@testset "nc_write" begin
  cellsize = 2
  range = [70, 140, 15, 55]

  lon = range[1]+cellsize/2:cellsize:range[2]
  lat = range[3]+cellsize/2:cellsize:range[4]

  nlon = length(lon)
  nlat = length(lat)
  ntime = 10
  dat = rand(Float64, nlon, nlat, ntime)

  # time = 1:size(dat2, 3)
  dims = [
    NcDim("lon", lon, Dict("longname" => "Longitude", "units" => "degrees east"))
    NcDim("lat", lat, Dict("longname" => "Latitude", "units" => "degrees north"))
    NcDim("time", 1:ntime)
  ]
  # timatts = Dict("longname" => "Time",
  #           "units"    => "hours since 01-01-2000 00:00:00");
  # ntime = size(dat)[3]
  fn = "temp_HI.nc"
  isfile(fn) && rm(fn)

  ## test for nc_write `type`
  nc_write(dat, fn, dims, Dict("longname" => "Heatwave Index");
    varname = "HI", overwrite = true,
    type = Float32)
  @test nc_read(fn) |> eltype == Float32
  @test nc_read(fn, type = Float64) |> eltype == Float64

  ## test for overwrite
  nc_write(dat, fn, dims, Dict("longname" => "Heatwave Index"); varname = "HI", overwrite = true)
  nc_write!(dat, fn, dims; varname = "HI2") # test for multiple variables

  nc_info(fn)
  @test nc_bands(fn) == ["HI", "HI2"]

  data = nc_read(fn, "HI")
  @test data == dat

  isfile(fn) && rm(fn)
end


text:
function EuropeanSpreadOptionStVolMC(K, T, S1, S2, sig1, sig2, div1, div2, alpha1, alpha2,
        Vbar1, Vbar2, xi1, xi2, rhoz, r, N, M)

        #rhoz to be entered as a korrelation matrix

        #N = 1

        dt = T/N
        alpha1dt = alpha1 * dt
        alpha2dt = alpha2 * dt

        xi1sdt = xi1 * sqrt(dt)
        xi2sdt = xi2 * sqrt(dt)
        sdt = sqrt(dt)

        lnS1 = log(S1)
        lnS2 = log(S2)

        z = zeros(Float64, 4)

        ## variable to store the sum of the option values
        sum_CT = 0
        sum_CT2 = 0

        ## outer loop number of simulations

        for j = 1:M

                lnSt1 = lnS1
                lnSt2 = lnS2
                Vt1 = Vbar1
                Vt2 = Vbar2

                for i = 1:N

                        # require generateCorrRandom.jl
                        z = generateCorrRandom(rhoz)

                        Vt1 = Vt1 + alpha1dt * (Vbar1 - Vt1) + xi1sdt * sqrt(Vt1) * z[3]
                        Vt2 = Vt2 + alpha2dt * (Vbar2 - Vt2) + xi2sdt * sqrt(Vt2) * z[4]

                        lnSt1 = lnSt1 + (r-div1-0.5*Vt1)*dt + sqrt(Vt1) * sdt * z[1]
                        lnSt2 = lnSt2 + (r-div2-0.5*Vt2)*dt + sqrt(Vt2) * sdt * z[2]

                end

                St1 = exp(lnSt1)
                St2 = exp(lnSt2)
                CT = max(0 , St1 - St2 - K)
                sum_CT = sum_CT + CT
                sum_CT2 = sum_CT2 + CT*CT

        end

        call_value = sum_CT / M * exp(-r * T)
        SD = sqrt((sum_CT2 - sum_CT * sum_CT / M ) * exp(-2 * r * T) / (M - 1))
        SE = SD / sqrt(M)

end


text:
#!/usr/bin/env julia
# -*- coding: utf-8 -*-

# ============================================================================
# Brille - Braille alphabet encode/decode module for Julia
# Copyright (C) 2019 by Ralf Kilian
# Distributed under the MIT License (https://opensource.org/licenses/MIT)
#
# GitHub: https://github.com/urbanware-org/brille
# GitLab: https://gitlab.com/urbanware-org/brille
# ============================================================================

__version__ = "1.0.0"

BRAILLE_ALPHABET = Dict([
                    ("a", "⠁"), ("b", "⠃"), ("c", "⠉"), ("d", "⠙"),
                    ("e", "⠑"), ("f", "⠋"), ("g", "⠛"), ("h", "⠓"),
                    ("i", "⠊"), ("j", "⠚"), ("k", "⠅"), ("l", "⠇"),
                    ("m", "⠍"), ("n", "⠝"), ("o", "⠕"), ("p", "⠏"),
                    ("q", "⠟"), ("r", "⠗"), ("s", "⠎"), ("t", "⠞"),
                    ("u", "⠥"), ("v", "⠧"), ("w", "⠺"), ("x", "⠭"),
                    ("y", "⠽"), ("z", "⠵"), ("1", "⠼⠁"), ("2", "⠼⠃"),
                    ("3", "⠼⠉"), ("4", "⠼⠙"), ("5", "⠼⠑"), ("6", "⠼⠋"),
                    ("7", "⠼⠛"), ("8", "⠼⠓"), ("9", "⠼⠊"), ("0", "⠼⠚")])

function get_version()
    """
        Return the version of this module.
    """
    return __version__
end

function string2braille(string)
    """
        Convert (encode) a string into the Braille alphabet.
    """
    output = ""

    while occursin("  ", string)
        string = replace(string, "  " => " ")
    end

    for char in lowercase(strip(string))
        try
            char = BRAILLE_ALPHABET["$char"]
            output = "$output$char"
        catch
            output = "$output$char"
        end
    end

    return strip(output)
end

function braille2string(braille_code)
    """
        Convert (decode) the Braille alphabet into a string.
    """
    output = ""

    while occursin("  ", braille_code)
        braille_code = replace(braille_code, "  " => " ")
    end

    braille = Dict(value => key for (key, value) in BRAILLE_ALPHABET)
    for char in braille_code
        try
            char = braille["$char"]
            output = "$output$char"
        catch
            output = "$output$char"
        end
    end

    return strip(output)
end

# EOF

max_stars_repo_name: urbanware-org/brille

max_issues_repo_name: urbanware-org/brille

max_forks_repo_name: urbanware-org/brille