text:
# Links pointing to the latest ConceptNetNumberbatch version (v"17.06")
const CONCEPTNET_MULTI_LINK = "https://conceptnet.s3.amazonaws.com/downloads/2017/numberbatch/numberbatch-17.06.txt.gz"
const CONCEPTNET_EN_LINK = "https://conceptnet.s3.amazonaws.com/downloads/2017/numberbatch/numberbatch-en-17.06.txt.gz"
const CONCEPTNET_HDF5_LINK = "https://conceptnet.s3.amazonaws.com/precomputed-data/2016/numberbatch/17.06/mini.h5"

# Accepted languages (map from conceptnet to Languages.Language)
const LANGUAGES = Dict(:en=>Languages.English(),
                       :fr=>Languages.French(),
                       :de=>Languages.German(),
                       :it=>Languages.Italian(),
                       :fi=>Languages.Finnish(),
                       :nl=>Languages.Dutch(),
                       :af=>Languages.Dutch(),
                       :pt=>Languages.Portuguese(),
                       :es=>Languages.Spanish(),
                       :ru=>Languages.Russian(),
                       :sh=>Languages.Serbian(),# and Languages.Croatian()
                       :sw=>Languages.Swedish(),
                       :cs=>Languages.Czech(),
                       :pl=>Languages.Polish(),
                       :bg=>Languages.Bulgarian(),
                       :eo=>Languages.Esperanto(),
                       :hu=>Languages.Hungarian(),
                       :el=>Languages.Greek(),
                       :no=>Languages.Nynorsk(),
                       :sl=>Languages.Slovene(),
                       :ro=>Languages.Romanian(),
                       :vi=>Languages.Vietnamese(),
                       :lv=>Languages.Latvian(),
                       :tr=>Languages.Turkish(),
                       :da=>Languages.Danish(),
                       :ar=>Languages.Arabic(),
                       :fa=>Languages.Persian(),
                       :ko=>Languages.Korean(),
                       :th=>Languages.Thai(),
                       :ka=>Languages.Georgian(),
                       :he=>Languages.Hebrew(),
                       :te=>Languages.Telugu(),
                       :et=>Languages.Estonian(),
                       :hi=>Languages.Hindi(),
                       :lt=>Languages.Lithuanian(),
                       :uk=>Languages.Ukrainian(),
                       :be=>Languages.Belarusian(),
                       :sw=>Languages.Swahili(),
                       :ur=>Languages.Urdu(),
                       :ku=>Languages.Kurdish(),
                       :az=>Languages.Azerbaijani(),
                       :ta=>Languages.Tamil()
                       # add more mappings here if needed
                       # AND supported by Languages.jl
                      )

# Regular expression on which to split text into tokens
const DEFAULT_SPLITTER = r"(,|\n|\r|\:|\\|\/|;|\.|\[|\]|\{|\}|\'|\`|\"|\"|\?|\!|\=|\~|\&|\s+)"


text:
using Unitful, StaticArrays, RecursiveArrayTools, LightLattices, Test

###############################################################
#
#   Cell tests
#
#

tcell = TrivialCell{3, Float64}()
@testset "Trivial cell" begin
    @test_throws BoundsError tcell[2]
    @test tcell[1] == SVector{3, Float64}(0.0,0.0,0.0)
    @test relative_coordinate(tcell, 1, 1) == SVector{3, Float64}(0.0,0.0,0.0)
	@test num_of_groups(tcell) == 1
	@test group_size(tcell, 1) == 1
	@test_throws ErrorException group_size(tcell, 2)
end

dcell = HomogeneousCell([[0,0,0],[0.25,0.25,0.25]], :diamond)
dcell_tuples = HomogeneousCell([(0,0.0,0),(0.25,0.25,0.25)], :diamond)
dcell_unlabeled = HomogeneousCell([[0,0,0],[0.25,0.25,0.25]])
@testset "Homogeneous unit cell of diamond" begin
    @test_throws BoundsError dcell[3]
    @test dcell[2] == SVector{3,Float64}(0.25,0.25,0.25)
    @test relative_coordinate(dcell, 2, 1) == SVector{3,Float64}(0.25,0.25,0.25)
	@test typeof(dcell_tuples) == typeof(dcell)
	@test dcell_tuples.cell_vectors == dcell.cell_vectors
	@test dcell_unlabeled.label == nothing
	@test num_of_groups(dcell) == 1
	@test group_size(dcell, 1) == 2
	@test_throws ErrorException group_size(dcell, 2)
end

### Fluorapatite basis_cell requires introduction of primitive vectors.

const af = 9.462u"Å"
const cf = 6.849u"Å"

const x = 0.369
const y = 0.389

fpvecs = hcat(af*[0.5, 0.5*sqrt(3), 0.0],
              af*[0.5, -0.5*sqrt(3), 0.0],
              cf*[0.0, 0.0, 1.0]
             ) |> SMatrix{3,3}
cell_vectors_raw1 = [[0.0, 0.0, 0.25], [0.0, 0.0, 0.75]]
cell_vectors_raw2 = [[x, y, 0.25], [-y, x-y, 0.25], [y-x, -x, 0.25],
                  [-x, -y, 0.75], [y, y-x, 0.75], [x-y, x, 0.75]]
fcell_vectors = ArrayPartition([fpvecs*vec for vec in cell_vectors_raw1], [fpvecs*vec for vec in cell_vectors_raw2])
fcell = InhomogeneousCell([fpvecs*vec for vec in cell_vectors_raw1], [fpvecs*vec for vec in cell_vectors_raw2]; label = :fluorapatite_magnetic)
@testset "Inhomogeneous unit cell for magnetic sublattice of fluorapatite" begin
    @test_throws BoundsError fcell[9]
    @test_throws BoundsError fcell[3,1]
    @test_throws BoundsError fcell[7, 2]
    @test fcell[2] == fcell_vectors[2]
    @test fcell[2, 1] == fcell_vectors[2]
    @test fcell[3] == fcell_vectors[2,1]
    @test fcell[1,2] == fcell_vectors[2,1]
    @test relative_coordinate(fcell, 5, 2) == fcell_vectors[2,3] - fcell_vectors[1,2]
    @test relative_coordinate(fcell, (3,2), 2) == fcell_vectors[2,3] - fcell_vectors[1,2]
    @test relative_coordinate(fcell, 5, (2,1)) == fcell_vectors[2,3] - fcell_vectors[1,2]
    @test relative_coordinate(fcell, (3,2), (2,1)) == fcell_vectors[2,3] - fcell_vectors[1,2]
	@test num_of_groups(fcell) == 2
	@test group_size(fcell, 1) == 2
	@test group_size(fcell, 2) == 6
	@test_throws ErrorException group_size(fcell, 3)
end

### Tests for switch_coord_type

mock_h_cell = HomogeneousCell([[0,0,0], [1,1,1]])
mock_ih_cell = InhomogeneousCell([[0,0,0]], [[1,1,1]])
tcell_transformed = switch_coord_type(tcell, Int)
mock_h_cell_transformed = switch_coord_type(mock_h_cell, Float64)
mock_ih_cell_transformed = switch_coord_type(mock_ih_cell, Float64)
@testset "Switching of coordinate type" begin
	@test switch_coord_type(tcell, Float64) == tcell
	@test switch_coord_type(mock_h_cell, Int) === mock_h_cell
	@test switch_coord_type(mock_ih_cell, Int) === mock_ih_cell
    @test typeof(tcell_transformed) == TrivialCell{3,Int}
    @test typeof(mock_h_cell_transformed) == HomogeneousCell{3, Float64, Nothing}
    @test typeof(mock_ih_cell_transformed) <: InhomogeneousCell{3,Float64, 2, Nothing}
    @test tcell_transformed[1] == SVector(0,0,0)
    @test mock_h_cell_transformed[2] == SVector(1.0,1.0,1.0)
    @test mock_ih_cell_transformed[1,2] == SVector(1.0,1.0,1.0)
end

###############################################################
#
#   Lattice tests
#
#

### Fluorine sublattice of CaF2 is a cubic lattice.

cubic_lattice_f = RegularLattice((11,11,11), 2.725u"Å"; periodic=false)
cubic_lattice_p = RegularLattice((11,11,11), 2.725u"Å")
@testset "Fluorine sublattice of CaF2 as an example of cubic lattice" begin
    @test cubic_lattice_f.central_cell == CartesianIndex(6,6,6)
    @test_throws BoundsError cubic_lattice_f[11,11,12]
    @test_throws BoundsError cubic_lattice_f[11,11,11,2]
    @test_throws BoundsError cubic_lattice_f[CartesianIndex(11,11,11),2]
    @test cubic_lattice_p[13,24,182] ≈ 2.725u"Å"*SVector(2, 2, 6)
    @test cubic_lattice_p[13,24,182] ≈ 2.725u"Å"*SVector(2, 2, 6)
    @test cubic_lattice_p[CartesianIndex(13,24,182)] ≈ 2.725u"Å"*SVector(2,2,6)
    @test cubic_lattice_p[CartesianIndex(13,24,182), 1] ≈ 2.725u"Å"*SVector(2,2,6)
    @test relative_coordinate(cubic_lattice_f, CartesianIndex(2,2,6), CartesianIndex(1,1,1)) ≈ 2.725u"Å"*SVector(1,1,5)
    @test relative_coordinate(cubic_lattice_f, (CartesianIndex(2,2,6),1), (CartesianIndex(1,1,1),1)) ≈ 2.725u"Å"*SVector(1,1,5)
    @test relative_coordinate(cubic_lattice_f, (2,2,6,1), (1,1,1,1)) ≈ 2.725u"Å"*SVector(1,1,5)
    @test relative_coordinate(cubic_lattice_p, CartesianIndex(13,24,182), CartesianIndex(1,1,1)) ≈ 2.725u"Å"*SVector(1,1,5)
    @test relative_coordinate(cubic_lattice_p, (CartesianIndex(13,24,182),1), (CartesianIndex(1,1,1),1)) ≈ 2.725u"Å"*SVector(1,1,5)
    @test relative_coordinate(cubic_lattice_p, (13,24,182,1), (1,1,1,1)) ≈ 2.725u"Å"*SVector(1,1,5)
	@test num_of_groups(cubic_lattice_p) == 1
	@test group_size(cubic_lattice_p, 1) == 11^3
	@test_throws ErrorException group_size(cubic_lattice_p, 2)
	@test length(cubic_lattice_p) == 11^3
end

### We consider dimensionless diamond lattice here.

dpvecs = 0.5*hcat([0,1,1],[1,1,0],[1,0,1]) |> SMatrix{3,3}
diamond_lattice_f = RegularLattice((11,11,11), dpvecs, dcell; label = :fcc, periodic = false)
diamond_lattice_p = RegularLattice((11,11,11), dpvecs, dcell; label = :fcc)
I1 = (CartesianIndex(13,24,182),1)
I2 = (CartesianIndex(1,1,1),2)
I1t = (13,24,182,1)
I2t = (1,1,1,2)
@testset "Dimensionless diamond lattice as an example of lattice with homogeneous cell" begin
    @test_throws BoundsError diamond_lattice_f[13,24,182,1]
    @test_throws BoundsError diamond_lattice_p[13,24,182,3]
    @test_throws BoundsError diamond_lattice_f[CartesianIndex(13,24,182),1]
    @test_throws BoundsError diamond_lattice_p[CartesianIndex(13,24,182),3]
    @test diamond_lattice_p[CartesianIndex(13,24,182), 1] == dpvecs*SVector(2,2,6)
    @test diamond_lattice_p[CartesianIndex(13,24,182), 2] == dpvecs*SVector(2,2,6) + SVector(0.25,0.25,0.25)
    @test diamond_lattice_p[13,24,182, 1] == dpvecs*SVector(2,2,6)
    @test diamond_lattice_p[13,24,182, 2] == dpvecs*SVector(2,2,6) + SVector(0.25,0.25,0.25)
    @test relative_coordinate(diamond_lattice_p, I1, I2) == dpvecs*SVector(1,1,5) - SVector(0.25,0.25,0.25)
    @test relative_coordinate(diamond_lattice_p, I1, I2) == -relative_coordinate(diamond_lattice_p, I2, I1)
    @test relative_coordinate(diamond_lattice_p, I1t, I2t) == dpvecs*SVector(1,1,5) - SVector(0.25,0.25,0.25)
    @test relative_coordinate(diamond_lattice_p, I1t, I2t) == -relative_coordinate(diamond_lattice_p, I2, I1)
	@test num_of_groups(diamond_lattice_p) == 1
	@test group_size(diamond_lattice_p, 1) == 11^3*2
	@test_throws ErrorException group_size(diamond_lattice_p, 2)
	@test length(diamond_lattice_p) == 11^3*2
end

### Fluorapatite lattice

fluorapatite_lattice_p = RegularLattice((11,11,11), fpvecs, fcell; label = :hexagonal)
I1 = (CartesianIndex(13,24,182), 5,2)
I2 = (CartesianIndex(1,1,1), 2,1)
I1t = (13,24,182, 5,2)
I2t = (1,1,1, 2,1)
@testset "Fluorapatite magnetic sublattice as an example of lattice with inhomogeneous cell." begin
    @test_throws BoundsError fluorapatite_lattice_p[CartesianIndex(13,24,182), 3,1]
    @test_throws BoundsError fluorapatite_lattice_p[13,24,182, 3,1]
    @test fluorapatite_lattice_p[CartesianIndex(13,24,182), 5,2] == fpvecs*SVector(2,2,6) + fcell[5,2]
    @test fluorapatite_lattice_p[13,24,182, 5,2] == fpvecs*SVector(2,2,6) + fcell[5,2]
    @test relative_coordinate(fluorapatite_lattice_p, I1, I2) ≈ fpvecs*SVector(1,1,5) + fcell[5,2] - fcell[2,1]
    @test relative_coordinate(fluorapatite_lattice_p, I1, I2) ≈ - relative_coordinate(fluorapatite_lattice_p, I2, I1)
    @test relative_coordinate(fluorapatite_lattice_p, I1t, I2t) ≈ fpvecs*SVector(1,1,5) + fcell[5,2] - fcell[2,1]
    @test relative_coordinate(fluorapatite_lattice_p, I1t, I2t) ≈ - relative_coordinate(fluorapatite_lattice_p, I2t, I1t)
	@test num_of_groups(fluorapatite_lattice_p) == 2
	@test group_size(fluorapatite_lattice_p, 1) == 11^3*2
	@test group_size(fluorapatite_lattice_p, 2) == 11^3*6
	@test_throws ErrorException group_size(fluorapatite_lattice_p, 3)
	@test length(fluorapatite_lattice_p) == 11^3*8
end

### Sorting utilities

using LightLattices: takes_precedence

@testset "Comparing of the indices" begin
	@test takes_precedence(1,2) == true
	@test takes_precedence(2,1) == false

	@test takes_precedence(CartesianIndex(1), CartesianIndex(2)) == true
	@test takes_precedence(CartesianIndex(2), CartesianIndex(1)) == false

	@test takes_precedence(CartesianIndex(2,3,1), CartesianIndex(1,1,2)) == true
	@test takes_precedence(CartesianIndex(1,1,2), CartesianIndex(2,3,1)) == false

	@test takes_precedence((CartesianIndex(1,1,2), 3), (CartesianIndex(2,3,1),4)) == true
	@test takes_precedence((CartesianIndex(2,3,1),4), (CartesianIndex(1,1,2), 3)) == false

	@test takes_precedence((1,1,2, 3), (2,3,1,4)) == true
	@test takes_precedence((2,3,1,4), (1,1,2, 3)) == false

	@test takes_precedence((CartesianIndex(2,3,1),4,1), (CartesianIndex(1,1,2), 3, 2)) == true
	@test takes_precedence((CartesianIndex(1,1,2), 3, 2), (CartesianIndex(2,3,1),4,1)) == false

	@test takes_precedence((2,3,1,4,1), (1,1,2, 3, 2)) == true
	@test takes_precedence((1,1,2, 3, 2), (2,3,1,4,1)) == false
end

ext: jl


text:
# include("admm_test.jl")

# test_admm_with_autodiff_smallinput(gd_iter=3,dbg=true)
# test_admm_without_autodiff_smallinput(gd_iter=3)
#test_admm_without_autodiff_largeinput(gd_iter=3,dbg=true)

using LinearAlgebra

# @testset "$(format("ADMM Algorithm: Small Input[Gaussian + Bernoulli]"))" begin
#   let
#     truth_matrix     = rand([(FixedRankMatrix(Distributions.Gaussian(5, 10), rank = 5), 200, 100),
#                              (FixedRankMatrix(Distributions.Bernoulli(0.5), rank = 5),  200, 100)])
#     sample_model     = provide(Sampler{BernoulliModel}(), rate = 0.8)
#     input_matrix     = sample_model.draw(truth_matrix)
#     display(input_matrix)
#     @time completed_matrix, type_tracker = complete(A      = input_matrix,
#                                               maxiter      = 200,
#                                               ρ            = 0.3,
#                                               use_autodiff = false,
#                                               gd_iter      = 3,
#                                               debug_mode   = false)
#     predicted_bernoulli = predict(Bernoulli(),
#                                   forward_map(Bernoulli(), completed_matrix[type_tracker[:Bernoulli]]))
#     truth_bernoulli = truth_matrix[type_tracker[:Bernoulli]]
#     summary_bernoulli = provide(Diagnostics{Any}(),
#                                 reference = truth_bernoulli,
#                                 input_data = predicted_bernoulli)
#     display(summary_bernoulli)
#     predicted_gaussian = predict(Gaussian(),
#                                  forward_map(Gaussian(), completed_matrix[type_tracker[:Gaussian]]))
#     truth_gaussian = truth_matrix[type_tracker[:Gaussian]]
#     summary_gaussian = provide(Diagnostics{Any}(),
#                                reference = truth_gaussian,
#                                input_data = predicted_gaussian)
#     display(summary_gaussian)
#   end
# end

# @testset "$(format("ADMM Algorithm: Small Input[Gaussian + Poisson]"))" begin
#   let
#     truth_matrix     = rand([(FixedRankMatrix(Distributions.Gaussian(5, 10), rank = 5), 200, 100),
#                              (FixedRankMatrix(Distributions.Poisson(10), rank = 10), 200, 100)])
#     sample_model     = provide(Sampler{BernoulliModel}(), rate = 0.8)
#     input_matrix     = sample_model.draw(truth_matrix)
#     display(input_matrix)
#     completed_matrix, type_tracker = complete(A            = input_matrix,
#                                               maxiter      = 200,
#                                               ρ            = 0.3,
#                                               use_autodiff = false,
#                                               gd_iter      = 5,
#                                               debug_mode   = false)
#     predicted_poisson = predict(Poisson(),
#                                 forward_map(Poisson(), completed_matrix[type_tracker[:Poisson]]))
#     truth_poisson = truth_matrix[type_tracker[:Poisson]]
#     summary_poisson = provide(Diagnostics{Any}(),
#                               reference = truth_poisson,
#                               input_data = predicted_poisson)
#     display(summary_poisson)
#     predicted_gaussian = predict(Gaussian(),
#                                  forward_map(Gaussian(), completed_matrix[type_tracker[:Gaussian]]))
#     truth_gaussian = truth_matrix[type_tracker[:Gaussian]]
#     summary_gaussian = provide(Diagnostics{Any}(),
#                                reference = truth_gaussian,
#                                input_data = predicted_gaussian)
#     display(summary_gaussian)
#   end
# end

# @testset "$(format("ADMM Algorithm: Small Input[Gaussian + Poisson + Bernoulli]"))" begin
#   let
#     truth_matrix     = rand([(FixedRankMatrix(Distributions.Gaussian(5, 10), rank = 5), 300, 100),
#                              (FixedRankMatrix(Distributions.Poisson(10), rank = 5), 300, 100),
#                              (FixedRankMatrix(Distributions.Bernoulli(0.5), rank = 5), 300, 100)])
#     sample_model     = provide(Sampler{BernoulliModel}(), rate = 0.8)
#     input_matrix     = sample_model.draw(truth_matrix)
#     display(input_matrix)
#     completed_matrix, type_tracker = complete(A            = input_matrix,
#                                               maxiter      = 200,
#                                               ρ            = 0.3,
#                                               use_autodiff = false,
#                                               gd_iter      = 5,
#                                               debug_mode   = false)
#     predicted_poisson = predict(Poisson(),
#                                 forward_map(Poisson(), completed_matrix[type_tracker[:Poisson]]))
#     truth_poisson = truth_matrix[type_tracker[:Poisson]]
#     summary_poisson = provide(Diagnostics{Any}(),
#                               reference = truth_poisson,
#                               input_data = predicted_poisson)
#     display(summary_poisson)
#     predicted_gaussian = predict(Gaussian(),
#                                  forward_map(Gaussian(), completed_matrix[type_tracker[:Gaussian]]))
#     truth_gaussian = truth_matrix[type_tracker[:Gaussian]]
#     summary_gaussian = provide(Diagnostics{Any}(),
#                                reference = truth_gaussian,
#                                input_data = predicted_gaussian)
#     display(summary_gaussian)

#     predicted_bernoulli = predict(Bernoulli(),
#                                   forward_map(Bernoulli(), completed_matrix[type_tracker[:Bernoulli]]))
#     truth_bernoulli = truth_matrix[type_tracker[:Bernoulli]]
#     summary_bernoulli = provide(Diagnostics{Any}(),
#                                 reference = truth_bernoulli,
#                                 input_data = predicted_bernoulli)
#     display(summary_bernoulli)
#   end
# end

# @testset "$(format("ADMM Algorithm: Medium Input[Gaussian + Poisson + Bernoulli]"))" begin
#   let
#     truth_matrix     = rand([(FixedRankMatrix(Distributions.Gaussian(5, 10), rank = 5), 300, 100),
#                              (FixedRankMatrix(Distributions.Poisson(10), rank = 5), 300, 100),
#                              (FixedRankMatrix(Distributions.Bernoulli(0.5), rank = 5), 300, 100)])
#     sample_model     = provide(Sampler{BernoulliModel}(), rate = 0.8)
#     input_matrix     = sample_model.draw(truth_matrix)
#     display(input_matrix)
#     completed_matrix, type_tracker = complete(A            = input_matrix,
#                                               maxiter      = 200,
#                                               ρ            = 0.3,
#                                               use_autodiff = false,
#                                               gd_iter      = 3,
#                                               debug_mode   = false)
#     predicted_poisson = predict(Poisson(),
#                                 forward_map(Poisson(), completed_matrix[type_tracker[:Poisson]]))
#     truth_poisson = truth_matrix[type_tracker[:Poisson]]
#     summary_poisson = provide(Diagnostics{Any}(),
#                               reference = truth_poisson,
#                               input_data = predicted_poisson)

#     predicted_gaussian = predict(Gaussian(),
#                                  forward_map(Gaussian(), completed_matrix[type_tracker[:Gaussian]]))
#     truth_gaussian = truth_matrix[type_tracker[:Gaussian]]
#     summary_gaussian = provide(Diagnostics{Any}(),
#                                reference = truth_gaussian,
#                                input_data = predicted_gaussian)

#     predicted_bernoulli = predict(Bernoulli(),
#                                   forward_map(Bernoulli(), completed_matrix[type_tracker[:Bernoulli]]))
#     truth_bernoulli = truth_matrix[type_tracker[:Bernoulli]]
#     summary_bernoulli = provide(Diagnostics{Any}(),
#                                 reference = truth_bernoulli,
#                                 input_data = predicted_bernoulli)
#     @info("Gaussian")
#     display(summary_gaussian)
#     @info("Poisson")
#     display(summary_poisson)
#     @info("Bernoulli")
#     display(summary_bernoulli)
#   end
# end

function relative_l2_error(x, y)
  return LinearAlgebra.norm(x -  y, 2)^2 / LinearAlgebra.norm(y)^2
end

# @testset "$(format("ADMM Algorithm: Small Input[Gaussian + Poisson + Bernoulli + Gamma]"))" begin
#   let
#     truth_matrix     = rand([(FixedRankMatrix(Distributions.Gaussian(5, 10), rank = 3), 400, 100),
#                              (FixedRankMatrix(Distributions.Poisson(10), rank = 3), 400, 100),
#                              (FixedRankMatrix(Distributions.Bernoulli(0.5), rank = 3), 400, 100),
#                              (FixedRankMatrix(Distributions.Gamma(10, 2), rank = 3), 400, 100)])
#     sample_model     = provide(Sampler{BernoulliModel}(), rate = 0.8)
#     input_matrix     = sample_model.draw(truth_matrix)
#     display(input_matrix)
#     completed_matrix, type_tracker = complete(A            = input_matrix,
#                                               maxiter      = 200,
#                                               ρ            = 0.3,
#                                               use_autodiff = false,
#                                               gd_iter      = 3,
#                                               debug_mode   = false)
#     predicted_poisson = predict(Poisson(),
#                                 forward_map(Poisson(), completed_matrix[type_tracker[:Poisson]]))
#     truth_poisson = truth_matrix[type_tracker[:Poisson]]
#     summary_poisson = provide(Diagnostics{Any}(),
#                               reference = truth_poisson,
#                               input_data = predicted_poisson)
#     predicted_gaussian = predict(Gaussian(),
#                                  forward_map(Gaussian(), completed_matrix[type_tracker[:Gaussian]]))
#     truth_gaussian = truth_matrix[type_tracker[:Gaussian]]
#     summary_gaussian = provide(Diagnostics{Any}(),
#                                reference = truth_gaussian,
#                                input_data = predicted_gaussian)
#     predicted_gamma = predict(Gamma(),
#                               forward_map(Gamma(), completed_matrix[type_tracker[:Gamma]]))
#     truth_gamma = truth_matrix[type_tracker[:Gamma]]
#     summary_gamma = provide(Diagnostics{Any}(),
#                             reference = truth_gamma,
#                             input_data = predicted_gamma)
#     predicted_bernoulli = predict(Bernoulli(),
#                                   forward_map(Bernoulli(), completed_matrix[type_tracker[:Bernoulli]]))
#     truth_bernoulli = truth_matrix[type_tracker[:Bernoulli]]
#     summary_bernoulli = provide(Diagnostics{Any}(),
#                                 reference = truth_bernoulli,
#                                 input_data = predicted_bernoulli)
#     @info("Gaussian")
#     display(summary_gaussian)
#     @info("Poisson")
#     display(summary_poisson)
#     @info("Bernoulli")
#     display(summary_bernoulli)
#     @info("Gamma")
#     display(summary_gamma)
#   end
# end

# # 1. rectangular matrices sufficient condition
# # 2. big scale simulation
# # 3. fix the small dimension
# # 4. time performance
# # 5. beroulli / uniform
# # 6. weighted missing pattern

# @testset "$(format("ADMM Algorithm: Small Input[Gaussian + Poisson + Bernoulli + Gamma]"))" begin
#   let
#     truth_matrix     = rand([(FixedRankMatrix(Distributions.Gaussian(5, 10), rank = 4), 400, 15),
#                              (FixedRankMatrix(Distributions.Poisson(10), rank = 4), 400, 15),
#                              (FixedRankMatrix(Distributions.Bernoulli(0.5), rank = 4), 400, 15),
#                              (FixedRankMatrix(Distributions.Gamma(10, 2), rank = 4), 400, 15)])
#     sample_model     = provide(Sampler{BernoulliModel}(), rate = 0.8)
#     input_matrix     = sample_model.draw(truth_matrix)
#     display(input_matrix)
#     completed_matrix, type_tracker = complete(A            = input_matrix,
#                                               maxiter      = 200,
#                                               ρ            = 0.3,
#                                               use_autodiff = false,
#                                               gd_iter      = 10,
#                                               debug_mode   = false)
#     predicted_poisson = predict(Poisson(),
#                                 forward_map(Poisson(), completed_matrix[type_tracker[:Poisson]]))
#     truth_poisson = truth_matrix[type_tracker[:Poisson]]
#     summary_poisson = provide(Diagnostics{Any}(),
#                               reference = truth_poisson,
#                               input_data = predicted_poisson)
#     predicted_gaussian = predict(Gaussian(),
#                                  forward_map(Gaussian(), completed_matrix[type_tracker[:Gaussian]]))
#     truth_gaussian = truth_matrix[type_tracker[:Gaussian]]
#     summary_gaussian = provide(Diagnostics{Any}(),
#                                reference = truth_gaussian,
#                                input_data = predicted_gaussian)
#     predicted_gamma = predict(Gamma(),
#                               forward_map(Gamma(), completed_matrix[type_tracker[:Gamma]]))
#     truth_gamma = truth_matrix[type_tracker[:Gamma]]
#     summary_gamma = provide(Diagnostics{Any}(),
#                             reference = truth_gamma,
#                             input_data = predicted_gamma)
#     predicted_bernoulli = predict(Bernoulli(),
#                                   forward_map(Bernoulli(), completed_matrix[type_tracker[:Bernoulli]]))
#     truth_bernoulli = truth_matrix[type_tracker[:Bernoulli]]
#     summary_bernoulli = provide(Diagnostics{Any}(),
#                                 reference = truth_bernoulli,
#                                 input_data = predicted_bernoulli)
#     @info("Gaussian")
#     display(summary_gaussian)
#     @info("Poisson")
#     display(summary_poisson)
#     @info("Bernoulli")
#     display(summary_bernoulli)
#     @info("Gamma")
#     display(summary_gamma)
#   end
# end

# @testset "$(format("ADMM Algorithm: Medium Input[Gaussian + Poisson + Bernoulli + Gamma]"))" begin
#   let
#     truth_matrix     = rand([(FixedRankMatrix(Distributions.Gaussian(5, 10), rank = 3), 1600, 400),
#                              (FixedRankMatrix(Distributions.Poisson(10), rank = 3), 1600, 400),
#                              (FixedRankMatrix(Distributions.Bernoulli(0.5), rank = 3), 1600, 400),
#                              (FixedRankMatrix(Distributions.Gamma(10, 2), rank = 3), 1600, 400)])
#     sample_model     = provide(Sampler{BernoulliModel}(), rate = 0.8)
#     input_matrix     = sample_model.draw(truth_matrix)
#     display(input_matrix)
#     manual_type_input = Array{Symbol}(undef, 1600, 1600)
#     manual_type_input[:, 1:400]     .= :Gaussian
#     manual_type_input[:, 401:800]   .= :Poisson
#     manual_type_input[:, 801:1200]  .= :Bernoulli
#     manual_type_input[:, 1201:1600] .= :Gamma
#     completed_matrix, type_tracker = complete(A            = input_matrix,
#                                               maxiter      = 200,
#                                               ρ            = 0.3,
#                                               use_autodiff = false,
#                                               gd_iter      = 3,
#                                               debug_mode   = false,
#                                               type_assignment = manual_type_input)

#     predicted_poisson = predict(Poisson(),
#                                 forward_map(Poisson(), completed_matrix[type_tracker[:Poisson]]))
#     truth_poisson = truth_matrix[type_tracker[:Poisson]]
#     summary_poisson = provide(Diagnostics{Any}(),
#                               reference = truth_poisson,
#                               input_data = predicted_poisson)
#     predicted_gaussian = predict(Gaussian(),
#                                  forward_map(Gaussian(), completed_matrix[type_tracker[:Gaussian]]))
#     truth_gaussian = truth_matrix[type_tracker[:Gaussian]]
#     summary_gaussian = provide(Diagnostics{Any}(),
#                                reference = truth_gaussian,
#                                input_data = predicted_gaussian)
#     predicted_gamma = predict(Gamma(),
#                               forward_map(Gamma(), completed_matrix[type_tracker[:Gamma]]))
#     truth_gamma = truth_matrix[type_tracker[:Gamma]]
#     summary_gamma = provide(Diagnostics{Any}(),
#                             reference = truth_gamma,
#                             input_data = predicted_gamma)
#     predicted_bernoulli = predict(Bernoulli(),
#                                   forward_map(Bernoulli(), completed_matrix[type_tracker[:Bernoulli]]))
#     truth_bernoulli = truth_matrix[type_tracker[:Bernoulli]]
#     summary_bernoulli = provide(Diagnostics{Any}(),
#                                 reference = truth_bernoulli,
#                                 input_data = predicted_bernoulli)
#     @info("Gaussian")
#     display(summary_gaussian)
#     @info("Poisson")
#     display(summary_poisson)
#     @info("Bernoulli")
#     display(summary_bernoulli)
#     @info("Gamma")
#     display(summary_gamma)
#   end
# end

# @testset "$(format("ADMM Algorithm: Small Input[Gamma]"))" begin
#   using MatrixCompletion
#   import Distributions
#   import LinearAlgebra
#   truth_matrix     = rand([(FixedRankMatrix(Distributions.Gamma(10, 2), rank = 5), 1600, 400)])
#   sample_model     = provide(Sampler{BernoulliModel}(), rate = 0.8)
#   input_matrix     = sample_model.draw(truth_matrix)
#   manual_type_input = Array{Symbol}(undef, 1600, 400)
#   manual_type_input .= :Gamma
#   completed_matrix, type_tracker = complete(A            = input_matrix,
#                                             maxiter      = 200,
#                                             ρ            = 0.3,
#                                             use_autodiff = false,
#                                             gd_iter      = 10,
#                                             debug_mode   = false,
#                                             type_assignment = manual_type_input)

#   predicted_matrix = predict(Gamma(), forward_map(Gamma(), completed_matrix))
#   error_matrix = abs.(predicted_matrix - truth_matrix)
#   total_l2_error = LinearAlgebra.norm(error_matrix, 2)^2
#   relative_error = relative_l2_error(predicted_matrix, truth_matrix)
#   @info("completed")
#   display(completed_matrix[1:10, 1:10])
#   @info("predicted")
#   display(predicted_matrix[1:10, 1:10])
#   @info("truth")
#   display(truth_matrix[1:10, 1:10])
#   @info("error")
#   display(error_matrix[1:10, 1:10])
#   @show(relative_error)
#   @show(total_l2_error)

# end

# @testset "$(format("ADMM Algorithm: Small Input[Gaussian + Gamma]"))" begin
#   let
#     truth_matrix     = rand([(FixedRankMatrix(Distributions.Gaussian(5, 10), rank = 5), 200, 100),
#                              (FixedRankMatrix(Distributions.Gamma(5, 0.5), rank = 5), 200, 100)])
#     sample_model     = provide(Sampler{BernoulliModel}(), rate = 0.8)
#     input_matrix     = sample_model.draw(truth_matrix)
#     manual_type_input = Array{Symbol}(undef, 200, 200)
#     manual_type_input[:, 1:100] .= :Gaussian
#     manual_type_input[:, 101:200] .= :Gamma
#     display(input_matrix)
#     completed_matrix, type_tracker = complete(A            = input_matrix,
#                                               maxiter      = 200,
#                                               σ            = 0.3,
#                                               use_autodiff = false,
#                                               gd_iter      = 3,
#                                               debug_mode   = false,
#                                               type_assignment = manual_type_input)
#     predicted_gamma = predict(Gamma(),
#                                 forward_map(Gamma(), completed_matrix[type_tracker[:Gamma]]))
#     truth_gamma = truth_matrix[type_tracker[:Gamma]]
#     display(truth_gamma - predicted_gamma)
#     error_matrix = abs.(truth_gamma - predicted_gamma)
#     relative_error = LinearAlgebra.norm(error_matrix,2)^2 / LinearAlgebra.norm(truth_gamma) ^ 2

#     @show(relative_l2_error(predicted_gamma, truth_gamma))
#     # summary_gamma = provide(Diagnostics{Any}(),
#     #                           reference = truth_gamma,
#     #                           input_data = predicted_gamma)
#     # display(summary_gamma)
#     predicted_gaussian = predict(Gaussian(),
#                                  forward_map(Gaussian(), completed_matrix[type_tracker[:Gaussian]]))
#     truth_gaussian = truth_matrix[type_tracker[:Gaussian]]
#     summary_gaussian = provide(Diagnostics{Any}(),
#                                reference = truth_gaussian,
#                                input_data = predicted_gaussian)
#     display(summary_gaussian)
#   end
# end

include("./sub_test_runner_admm_negative_binomial.jl")
# include("./sub_test_runner_admm_bernoulli.jl")
# include("./sub_test_runner_admm_gaussian.jl")
# include("./sub_test_runner_admm_poisson.jl")
# include("./sub_test_runner_admm_gamma.jl")

# @testset "$(format("ADMM Algorithm: Small Input[Gaussian + Bernoulli, AutoDiff]"))" begin
#   let
#     truth_matrix     = rand([(FixedRankMatrix(Distributions.Gaussian(5, 10), rank = 2), 200, 100),
#                              (FixedRankMatrix(Distributions.Bernoulli(0.5), rank = 2), 200, 100)])
#     sample_model     = provide(Sampler{BernoulliModel}(), rate = 0.8)
#     input_matrix     = sample_model.draw(truth_matrix)
#     display(input_matrix)
#     completed_matrix, type_tracker = complete(A            = input_matrix,
#                                               maxiter      = 200,
#                                               σ            = 0.3,
#                                               use_autodiff = true,
#                                               gd_iter      = 3,
#                                               debug_mode   = false)
#     predicted_bernoulli = predict(Bernoulli(),
#                                   forward_map(Bernoulli(), completed_matrix[type_tracker[:Bernoulli]]))
#     truth_bernoulli = truth_matrix[type_tracker[:Bernoulli]]
#     summary_bernoulli = provide(Diagnostics{Any}(),
#                                 reference = truth_bernoulli,
#                                 input_data = predicted_bernoulli)
#     display(summary_bernoulli)
#     predicted_gaussian = predict(Gaussian(),
#                                  forward_map(Gaussian(), completed_matrix[type_tracker[:Gaussian]]))
#     truth_gaussian = truth_matrix[type_tracker[:Gaussian]]
#     summary_gaussian = provide(Diagnostics{Any}(),
#                                reference = truth_gaussian,
#                                input_data = predicted_gaussian)
#     display(summary_gaussian)
#   end
# end


text:
function run_precompile_projects(argv=ARGS)

    argset = ArgParse.ArgParseSettings()
    ArgParse.@add_arg_table! argset begin
        "rootdir"
            help="Root directory"
            arg_type = String
            default = pwd()
    end

    parsed_args = ArgParse.parse_args(argv, argset)
    rootdir = parsed_args["rootdir"]

    _walk_pkgs(;rootdir) do path, proj
        try;
                Pkg.activate(path)
                Pkg.precompile()
            catch err; @error("Precompilation Fails", err)
        end
    end
end


text:
using Distributions
using Distances
using Random
using LinearAlgebra
using Zygote
using Zygote: @adjoint

include("./kernel.jl")
include("./transform.jl")
include("./covariate.jl")

abstract type AbstractParameter{T} end

struct Parameter{T}
    param::Vector{T}
    probs::Vector{Float64}
end

mutable struct ScalarParameter{T<:Real} <: AbstractParameter{T}
    param::T
    dist::Distribution{Univariate,Continuous}
end

struct VectorParameter{T<:Real} <: AbstractParameter{T}
    param::Vector{T}
    dist::Vector{Distribution{Univariate,Continuous}}
end

logjacobian(g, λ, y) = sum(x -> log∂(g, λ, x), y)

Zygote.@adjoint function mul!(
    C::AbstractMatrix,
    A::AbstractMatrix,
    B::AbstractMatrix,
    α::Real,
    β::Real,
)
    tmp, back = pullback(C, A, B, α, β) do C, A, B, α, β
        (A * B) .* α .+ C .* β
    end
    C .= tmp
    return C, c::AbstractMatrix -> back(c)
end

Zygote.@adjoint function ldiv!(A, b)
    tmp, back = pullback(A, b) do A, b
        A \ b
    end
    b .= tmp
    return b, c::AbstractArray -> back(c)
end
Zygote.@adjoint function ldiv!(c, A, b)
    tmp, back = pullback(A, b) do A, b
        A \ b
    end
    c .= tmp
    return b, c::AbstractArray -> (nothing, back(c)...)
end

Zygote.@adjoint function cholesky!(A)
    tmp, back = pullback(A) do A
        cholesky(A)
    end
    return tmp, c::AbstractArray -> back(c)
end

function update_Kxx′!(
    Kxx′::AbstractMatrix,
    B::AbstractMatrix,
    Kux::AbstractMatrix,
    W::AbstractMatrix,
    W′::AbstractMatrix
)
    mul!(Kxx′, W', B, -1, 1)
    mul!(Kxx′, Kux', W′, -1, 1)
    return Kxx′
end

function update_Kx′!(
    Kx′::AbstractMatrix,
    B::AbstractMatrix,
    Kux′::AbstractMatrix,
    W′::AbstractMatrix
) # B = Kux′ .- Ku * W′
    mul!(Kx′, (W′)', B, -1, 1)
    mul!(Kx′, (Kux′)', W′, -1, 1)
    return Kx′
end

struct Border
    U::Matrix{Float64} # p observation locations
    Fu::LU{Float64, Matrix{Float64}} # p x p covariates at U
    X::Matrix{Float64} # n - p observation locations
    W::Matrix{Float64} # Fu \ Fx
    y::Vector{Float64} # length n vector
end
function border(f, U, X, y)
    Fu = lu!(f(U))
    W = ldiv!(Fu, f(X))
    return Border(U, Fu, X, W, y)
end

struct Block
    Ku::Matrix{Float64}
    Kux::Matrix{Float64}
    B::Matrix{Float64}
    rKx::Matrix{Float64}
end
function block(p, m)
    Ku = Matrix{Float64}(undef, p, p)
    Kux = Matrix{Float64}(undef, p, m)
    B = Matrix{Float64}(undef, p, m)
    rKx = Matrix{Float64}(undef, m, m)
    return Block(Ku, Kux, B, rKx)
end

struct RightHandSide
    z::Vector{Float64}
    rz::Vector{Float64}
    rc::Vector{Float64}
end
function righthandside(p, m)
    z = Vector{Float64}(undef, p + m)
    rz = Vector{Float64}(undef, m)
    rc = Vector{Float64}(undef, m)
    return RightHandSide(z, rz, rc)
end

function logabsdet(Fu, rKx)
    s = (-1.0) ^ (size(Fu, 1) % 2)
    return 2 * logabsdet(Fu)[1] + logdet(rKx), s
end

function reduced_system(k, θ, ϵ, br, bl)
    p, m = size(br.W)
    Ku = jitter!(covariance!(bl.Ku, k, θ, br.U, br.U), ϵ)
    Kux = covariance!(bl.Kux, k, θ, br.U, br.X)
    Kx = jitter!(covariance!(bl.rKx, k, θ, br.X, br.X), ϵ)
    B = mul!(copyto!(bl.B, Kux), bl.Ku, br.W, -1, 1)
    rKx = cholesky!(Symmetric(update_Kx′!(Kx, B, Kux, br.W)))
    return rKx
end

function solve_reduced(g, λ, rKx, br, rhs)
    p = size(br.U, 2)
    rhs.z .= (x -> transform(g, λ, x)).(br.y)
    @views rz = mul!(copyto!(rhs.rz, rhs.z[p+1:end]), br.W', rhs.z[1:p], -1, 1)
    return ldiv!(rKx, copyto!(rhs.rc, rz)), rz
end

function logprob(k, θ, ϵ, g, λ, br, bl, rhs)
    p, m = size(br.W)
    lnJ = logjacobian(g, λ, br.y)
    rKx = reduced_system(k, θ, ϵ, br, bl)
    rc, rz = solve_reduced(g, λ, rKx, br, rhs)
    lnQ = log(dot(rc, rz))
    lnDet = logabsdet(br.Fu, rKx)[1]

    return m / (m + p) * lnJ - m / 2 * lnQ - 1 / 2 * lnDet
end

max_stars_repo_name: dbindel/btg

max_issues_repo_name: dbindel/btg

max_forks_repo_name: dbindel/btg


text:
struct Telegraphers <: Equation end
@declare_dofs Telegraphers [:I, :V, :L, :C] 2

function evaluate_flux(eq::Telegraphers, celldofs, cellflux)
    s = TelegraphersShortcuts()
    cellflux .= 0.0
    L = celldofs[s.L]
    C = celldofs[s.C]
    for dimension=1:2
        if dimension == 1
            cellflux[1,s.I] = 1/L * celldofs[s.V]
            cellflux[1,s.V] = 1/C * celldofs[s.I]
        else
            # Only defined in x-dir!
            cellflux[2,1] = 0.0
            cellflux[2,2] = 0.0
        end
    end
end

function max_eigenval(eq::Telegraphers, celldofs, normalidx)
    s = TelegraphersShortcuts()
    L = celldofs[s.L]
    C = celldofs[s.C]
    return sqrt(1/(C*L))
end

function get_initial_values(eq::Telegraphers, scenario::PlanarWaves, global_position; t=0.0)
    pxg, pyg = global_position
    f1(x) = sin(2*pi * x)
    f2(x) = 0.1 * cos(2*pi * x)

    C = 2.0
    L = 1.5
    c = sqrt(1/(L*C)) # wavespeed
    Z_0 = sqrt(L/C) # impedance

    V = f1(pxg - c*t) + f2(pxg + c*t)
    I = f1(pxg - c*t)/Z_0 - f2(pxg + c*t)/Z_0

    [I, V, L, C]
end

function is_periodic_boundary(eq::Telegraphers, scenario::PlanarWaves)
    true
end

function is_analytical_solution(eq::Telegraphers, scenario::PlanarWaves)
    true
end


text:
using DataFrames
using Recommenders: leave_one_out_split, ratio_split
using Test

@testset "Leave one out split by timestamp" begin
    df = DataFrame(
        userid = [1, 1, 1, 2, 2, 3, 3, 3, 3],
        timestamp = [0, 2, 10, 3, 5, 4, 5, 11, -1],
    )
    df_train, df_test = leave_one_out_split(df, col_user = :userid, col_time = :timestamp)
    sort!(df_train, [:userid, :timestamp])
    sort!(df_test, [:userid, :timestamp])

    @test df_train ==
          DataFrame(userid = [1, 1, 2, 3, 3, 3], timestamp = [0, 2, 3, -1, 4, 5])
    @test df_test == DataFrame(userid = [1, 2, 3], timestamp = [10, 5, 11])
end

@testset "Ratio split" begin
    df = DataFrame((c = [1, 1, 1, 2, 2, 3, 3, 3, 3, 4],))
    df_train, df_test = ratio_split(df, 0.8)
    @test nrow(df_train) == 8
    @test nrow(df_test) == 2
end


text:
using Revise
using Random
using LinearAlgebra
using POMDPs
using POMDPModelTools
using POMDPPolicies
using POMDPSimulators
using IterTools
using AutomotiveDrivingModels
using AutomotivePOMDPs: polygon, overlap
using Parameters
using GridInterpolations
using StaticArrays
using AutoUrban
using AutoViz
using Printf
using Reel
using DiscreteValueIteration
using SparseArrays
using ProgressMeter
using JLD2
using FileIO

includet("environment.jl")
includet("mdp_definition.jl")
includet("helpers.jl")

mdp = MergingMDP()
@load "policies.jld2" policies
policy = policies[1] # [:wide, :cv, :ca, :cd]

function safe_actions(svec::Vector{Float64}, threshold::Float64=0.99)
    @assert length(svec) == length(first(states(mdp)))

    ## interpolate state
    sitp, witp = interpolants(mdp.state_grid, svec)
    normalize!(witp, 1)
    vals = zeros(n_actions(mdp))
    for (i, si) in enumerate(sitp)
        w = witp[i]
        s = ind2x(mdp.state_grid, si)
        if isapprox(svec[4], AbsentPos)
            s[4] = AbsentPos
        end
        vals += w .* actionvalues(policy, s)
    end
    return findall(vals .> threshold)
end

println("Loaded model checker")


text:
using DataFrames

A = readdlm("u1.base",'\t';header=false)
file = "u1.base"

function SpearmanCorrelation(x,y)
  @assert(length(x) == length(y))
  n = length(x)
  @assert(n>0)
	RankX    = Int64[]
	RankY    = Int64[]
  d        = Float64[]
  dsquared = Float64[]
	IndexX   = sortperm(x)
	IndexY   = sortperm(y)
  SumD2    = Float64{}
  SumD2    = 0.0
	for(i = 1: length(x))
		push!(RankX,0)
		push!(RankY,0)
    push!(d,0.0)
    push!(dsquared,0.0)
	end
	for(i = 1:length(x))
		RankX[IndexX[i]] = RankY[IndexY[i]] = length(x) - i + 1
	end
  for(i = 1:length(x))
    d[i]          = RankX[i] - RankY[i]
    dsquared[i]   = d * d
    SumD2        += dsquared[i]
  end
  spearCoeff = 1 - 6 * SumD2/(N^3 - N)
  return spearCoeff
end

function SpearmanCorrelationHelper(file)
     df = readtable(file,separator = '\t',header = false);
     rename!(df,:x1,:userId)
     rename!(df,:x2,:itemId)
     rename!(df,:x3,:rating)
     #rename!(df,:x4,:timestamp)
     M1 =Set(df[df[:userId] .== 1, :itemId])
     CorrelationArray = Float64[]
     (N,ind) = findmax(df[:,1])
     for(i = 2:N)
            M_i = Set(df[df[:userId] .== i, :itemId])
            if(length(intersect(M1,M_i)) == 0) push!( CorrelationArray , NA)
            else
                  X = Float64[]
                  Y = Float64[]
                  for(movie in intersect(M1,M_i))
                      push!(X , (df[(df[:userId]   .== 1)  & (df[:itemId] .== movie), :rating])[1])
                      push!(Y , (df[(df[:userId]   .== i)  & (df[:itemId] .== movie), :rating])[1])
                  end
                  push!( CorrelationArray , SpearmanCorrelation(X,Y))
            end
     end
end


text:
using Test
using ModelingToolkit
import ModelingToolkit: ExtraVariablesSystemException, ExtraEquationsSystemException

include("../examples/electrical_components.jl")

function UnderdefinedConstantVoltage(; name, V = 1.0)
    val = V
    @named p = Pin()
    @named n = Pin()
    @parameters V
    eqs = [
        V ~ p.v - n.v,
        # Remove equation
        # 0 ~ p.i + n.i
    ]
    ODESystem(eqs, t, [], [V], systems = [p, n], defaults = Dict(V => val), name = name)
end

function OverdefinedConstantVoltage(; name, V = 1.0, I = 1.0)
    val = V
    val2 = I
    @named p = Pin()
    @named n = Pin()
    @parameters V I
    eqs = [V ~ p.v - n.v
           # Overdefine p.i and n.i
           n.i ~ I
           p.i ~ I]
    ODESystem(eqs, t, [], [V], systems = [p, n], defaults = Dict(V => val, I => val2),
              name = name)
end

R = 1.0
C = 1.0
V = 1.0
@named resistor = Resistor(R = R)
@named capacitor = Capacitor(C = C)
@named source = UnderdefinedConstantVoltage(V = V)

rc_eqs = [connect(source.p, resistor.p)
          connect(resistor.n, capacitor.p)
          connect(capacitor.n, source.n)]

@named rc_model = ODESystem(rc_eqs, t, systems = [resistor, capacitor, source])
@test_throws ModelingToolkit.ExtraVariablesSystemException structural_simplify(rc_model)

@named source2 = OverdefinedConstantVoltage(V = V, I = V / R)
rc_eqs2 = [connect(source2.p, resistor.p)
           connect(resistor.n, capacitor.p)
           connect(capacitor.n, source2.n)]

@named rc_model2 = ODESystem(rc_eqs2, t, systems = [resistor, capacitor, source2])
@test_throws ModelingToolkit.ExtraEquationsSystemException structural_simplify(rc_model2)


text:
### I could use this file to set up the system for the test.

using DemeNet: DemeSpec, Deme, Signer, save, DemeID
using PeaceCypher
using PeaceFounder: BraiderConfig, RecorderConfig, CertifierConfig, BraidChainConfig, PeaceFounderConfig, AddressRecord, Port, certify, serialize, deserialize, config
using PeaceVote: AbstractChain

demespec = DemeSpec("PeaceDeme",:default,:PeaceCypher,:default,:PeaceCypher,:PeaceFounder)
save(demespec) ### Necessary to connect with Mixer
uuid = demespec.uuid
deme = Deme(demespec)

maintainer = Signer(deme,"maintainer")
server = Signer(deme,"server")

MIXER_ID = server.id
SERVER_ID = server.id
MAINTAINER_ID = maintainer.id

MIXER_PORT = Port(3001) # Self mixing
BRAIDER_PORT = Port(3000)
REGISTRATOR_PORT = Port(3002)
VOTING_PORT = Port(3003)
PROPOSAL_PORT = Port(3004)
SYNC_PORT = Port(3005)
TOOKEN_PORT = Port(3006)
CERTIFIER_PORT = Port(3007)

braiderconfig = BraiderConfig(BRAIDER_PORT,MIXER_PORT,UInt8(3),UInt8(64),SERVER_ID,DemeID(uuid,MIXER_ID))
recorderconfig = RecorderConfig([MAINTAINER_ID,SERVER_ID],server.id,REGISTRATOR_PORT,VOTING_PORT,PROPOSAL_PORT)
braidchainconfig = BraidChainConfig(SERVER_ID,MIXER_PORT,SYNC_PORT,braiderconfig,recorderconfig)
certifierconfig = CertifierConfig(MAINTAINER_ID,SERVER_ID,TOOKEN_PORT,CERTIFIER_PORT)

peacefounderconfig = PeaceFounderConfig(braidchainconfig,certifierconfig,AddressRecord[])

serialize(deme,peacefounderconfig)
certify(deme,maintainer)
deserialize(deme,PeaceFounderConfig)

braidchain = AbstractChain(deme)

max_stars_repo_name: akels/DemocracyCredit

max_issues_repo_name: akels/DemocracyCredit

max_forks_repo_name: akels/DemocracyCredit


text:
# The solver type
mutable struct ValueIterationSolver <: Solver
    max_iterations::Int64 # max number of iterations
    belres::Float64 # the Bellman Residual
    verbose::Bool
    include_Q::Bool
    init_util::Vector{Float64}
end
# Default constructor
function ValueIterationSolver(;max_iterations::Int64 = 100,
                               belres::Float64 = 1e-3,
                               verbose::Bool = false,
                               include_Q::Bool = true,
                               init_util::Vector{Float64}=Vector{Float64}(undef, 0))
    return ValueIterationSolver(max_iterations, belres, verbose, include_Q, init_util)
end

@POMDP_require solve(solver::ValueIterationSolver, mdp::Union{MDP,POMDP}) begin
    P = typeof(mdp)
    S = statetype(P)
    A = actiontype(P)
    @req discount(::P)
    @req n_states(::P)
    @req n_actions(::P)
    @subreq ordered_states(mdp)
    @subreq ordered_actions(mdp)
    @req transition(::P,::S,::A)
    @req reward(::P,::S,::A,::S)
    @req stateindex(::P,::S)
    @req actionindex(::P, ::A)
    @req actions(::P, ::S)
    as = actions(mdp)
    ss = states(mdp)
    a = first(as)
    s = first(ss)
    dist = transition(mdp, s, a)
    D = typeof(dist)
    @req support(::D)
    @req pdf(::D,::S)
end

#####################################################################
# Solve runs the value iteration algorithm.
# The policy input argument is either provided by the user or
# initialized during the function call.
# Verbose is a flag that triggers text output to the command line
# Example code for running the function:
# mdp = GridWorld(10, 10) # initialize a 10x10 grid world MDP (user written code)
# solver = ValueIterationSolver(max_iterations=40, belres=1e-3)
# policy = ValueIterationPolicy(mdp)
# solve(solver, mdp, policy, verbose=true)
#####################################################################
function solve(solver::ValueIterationSolver, mdp::Union{MDP,POMDP}; kwargs...)

    # deprecation warning - can be removed when Julia 1.0 is adopted
    if !isempty(kwargs)
        @warn("Keyword args for solve(::ValueIterationSolver, ::MDP) are no longer supported. For verbose output, use the verbose option in the ValueIterationSolver")
    end

    @warn_requirements solve(solver, mdp)

    # solver parameters
    max_iterations = solver.max_iterations
    belres = solver.belres
    discount_factor = discount(mdp)
    ns = n_states(mdp)
    na = n_actions(mdp)

    # intialize the utility and Q-matrix
    if !isempty(solver.init_util)
        @assert length(solver.init_util) == ns "Input utility dimension mismatch"
        util = solver.init_util
    else
        util = zeros(ns)
    end
    include_Q = solver.include_Q
    if include_Q
        qmat = zeros(ns, na)
    end
    pol = zeros(Int64, ns)

    total_time = 0.0
    iter_time = 0.0

    # create an ordered list of states for fast iteration
    states = ordered_states(mdp)

    # main loop
    for i = 1:max_iterations
        residual = 0.0
        iter_time = @elapsed begin
        # state loop
        for (istate,s) in enumerate(states)
            sub_aspace = actions(mdp, s)
            if isterminal(mdp, s)
                util[istate] = 0.0
                pol[istate] = 1
            else
                old_util = util[istate] # for residual
                max_util = -Inf
                # action loop
                # util(s) = max_a( R(s,a) + discount_factor * sum(T(s'|s,a)util(s') )
                for a in sub_aspace
                    iaction = actionindex(mdp, a)
                    dist = transition(mdp, s, a) # creates distribution over neighbors
                    u = 0.0
                    for (sp, p) in weighted_iterator(dist)
                        p == 0.0 ? continue : nothing # skip if zero prob
                        r = reward(mdp, s, a, sp)
                        u += p * r
                        if !isterminal(mdp, sp, a)
                            isp = stateindex(mdp, sp)
                            u += p * discount_factor * util[isp]
                        end
                    end
                    # BEFORE
                    # for (sp, p) in weighted_iterator(dist)
                    #     p == 0.0 ? continue : nothing # skip if zero prob
                    #     r = reward(mdp, s, a, sp)
                    #     isp = stateindex(mdp, sp)
                    #     u += p * (r + discount_factor * util[isp])
                    # end
                    new_util = u
                    if new_util > max_util
                        max_util = new_util
                        pol[istate] = iaction
                    end
                    include_Q ? (qmat[istate, iaction] = new_util) : nothing
                end # action
                # update the value array
                util[istate] = max_util
                diff = abs(max_util - old_util)
                diff > residual ? (residual = diff) : nothing
            end
        end # state
        end # time
        total_time += iter_time
        solver.verbose ? @printf("[Iteration %-4d] residual: %10.3G | iteration runtime: %10.3f ms, (%10.3G s total)\n", i, residual, iter_time*1000.0, total_time) : nothing
        residual < belres ? break : nothing
    end # main
    if include_Q
        return ValueIterationPolicy(mdp, qmat, util, pol)
    else
        return ValueIterationPolicy(mdp, utility=util, policy=pol, include_Q=false)
    end
end


text:
using StatsBase, Distances, Hungarian, Printf

export MatchingDistance, MatchDist
export FastMatchingDistance, FastMatchDist
export FixedPenaltyMatchingDistance, FixPenMatchDist
export MinDistMatchingDistance, MinDistMatchDist
export AvgSizeMatchingDistance, AvgSizeMatchDist
export get_cost_matrix_dynamic, get_cost_matrix_fixed

struct MatchingDistance{T<:Metric} <: Metric
    ground_dist::T
end

const MatchDist = MatchingDistance

function get_cost_matrix_dynamic(
    d::MatchingDistance,
    X::Vector{T}, Y::Vector{T}
    ) where {T}

    if length(X) < length(Y)
        return get_cost_matrix_dynamic(d,Y,X)
    elseif length(X)==length(Y)
        return pairwise_inbounds(d.ground_dist, X, Y)
    else
        C = pairwise_inbounds(d.ground_dist, X, Y)
        null_dists = [d.ground_dist(nothing, p) for p in X]
        size_diff = length(X)-length(Y)
        return [C [x for x∈null_dists, j=1:size_diff]]
    end
end

function get_cost_matrix_fixed(
    d::MatchingDistance,
    X::Vector{T}, Y::Vector{T}
    ) where {T}

    if length(X) < length(Y)
        C = pairwise_inbounds(d.ground_dist, X, Y)
        pentalty_vec = [d.ground_dist(nothing, p) for p in Y]
        size_diff = length(Y) - length(X)
        return [C ; [y for i=1:size_diff, y ∈ pentalty_vec]]
    elseif length(X)==length(Y)
        return pairwise_inbounds(d.ground_dist, X, Y)
    else
        C = pairwise_inbounds(d.ground_dist, X, Y)
        null_dists = [d.ground_dist(nothing, p) for p in X]
        size_diff = length(X)-length(Y)
        return [C [x for x∈null_dists, j=1:size_diff]]
    end

end

struct FastMatchingDistance{T<:Metric} <: Metric
    ground_dist::T
    C::Matrix{Float64}
    function FastMatchingDistance(ground_dist::S, K::Int) where {S<:Metric}
        new{S}(ground_dist, zeros(K,K))
    end
end

const FastMatchDist = FastMatchingDistance

Base.show(io::IO, d::FastMatchingDistance) = print(io, "$(typeof(d))(K=$(size(d.C,1)))")

function get_cost_matrix_dynamic(
    d::FastMatchingDistance,
    X::Vector{T}, Y::Vector{T}
    ) where {T}

    N,M = (length(X), length(Y))
    if N < M
        return get_cost_matrix_dynamic(d,Y,X)
    elseif N==M
        C = view(d.C, 1:N, 1:N)
        pairwise_inbounds!(C, d.ground_dist, X, Y)
    else
        C = view(d.C, 1:N, 1:N)
        pairwise_inbounds!(C, d.ground_dist, X, Y)
        for i in 1:N
            C[i,M+1] = d.ground_dist(nothing, X[i])
        end
        for j in (M+2):N
            for i in 1:N
                C[i,j] = C[i,j-1]
            end
        end
    end
    return C
end

function get_cost_matrix_fixed(
    d::FastMatchingDistance,
    X::Vector{T}, Y::Vector{T}
    ) where {T}

    N,M = (length(X), length(Y))
    if N < M
        C = view(d.C, 1:M, 1:M)
        pairwise_inbounds!(C, d.ground_dist, X, Y)
        for j in 1:M
            C[N+1,j] = d.ground_dist(nothing, Y[j])
        end
        for i in (N+2):N
            for j in 1:M
                C[i,j] = C[i-1,j]
            end
        end
    elseif N==M
        C = view(d.C, 1:N, 1:N)
        pairwise_inbounds!(C, d.ground_dist, X, Y)
    else
        C = view(d.C, 1:N, 1:N)
        pairwise_inbounds!(C, d.ground_dist, X, Y)
        for i in 1:N
            C[i,M+1] = d.ground_dist(nothing, X[i])
        end
        for j in (M+2):M
            for i in 1:N
                C[i,j] = C[i,j-1]
            end
        end
    end
    return C

end

function (d::Union{MatchDist,FastMatchDist})(X::Nothing, Y::Vector{T})::Float64 where {T}
    return sum(p->d.ground_dist(nothing,p), Y)
end
function (d::Union{MatchDist,FastMatchDist})(X::Vector{T}, Y::Nothing)::Float64 where {T}
    return d(Y, X)
end
function (d::Union{MatchDist,FastMatchDist})(X::Nothing, Y::Nothing)::Float64 where {T}
    return 0.0
end

struct FixedPenaltyMatchingDistance{T<:Metric} <: Metric
    ground_dist::T
    penalty::Float64
end

const FixPenMatchDist = FixedPenaltyMatchingDistance
# const FixPenMatchDist{T} = FixedPenaltyMatchingDistance{T}
# FixPenMatchDist(d_g::Metric, penalty::Float64) = FixedPenaltyMatchingDistance(d_g, penalty)

Base.show(io::IO, d::FixPenMatchDist) = print(io, "$(typeof(d))(ρ=$(d.penalty))")

function get_cost_matrix_dynamic(
    d::FixPenMatchDist,
    X::Vector{T}, Y::Vector{T}
    ) where {T}

    d_g = d.ground_dist
    N, M = (length(X), length(Y))
    if N < M
        return get_cost_matrix_dynamic(d, Y, X)
    elseif N == M
        return pairwise_inbounds(d_g, X, Y)
    else
        C = pairwise_inbounds(d_g, X, Y)
        size_diff = (N - M)
        return [C [d.penalty for i ∈ eachindex(X), j=1:size_diff]]
    end
end

function get_cost_matrix_fixed(
    d::FixPenMatchDist,
    X::Vector{T}, Y::Vector{T}
    ) where {T}

    d_g = d.ground_dist
    N, M = (length(X), length(Y))
    C = pairwise_inbounds(d_g, X, Y)
    if N < M
        size_diff = (M - N)
        return [C ; [d.penalty for i in 1:size_diff, j ∈ eachindex(Y)]]
    elseif N > M
        size_diff = (N - M)
        return [C [d.penalty for i ∈ eachindex(X), j in 1:size_diff]]
    else
        return pairwise_inbounds(d_g, X, Y)
    end

end

# function (d::FixPenMatchDist)(
#     S1::Vector{T}, S2::Vector{T}
#     ) where {T}

#     length(S1) < length(S2) ? (N_min, N_max)=(length(S1),length(S2)) : (N_min,N_max)=(length(S2),length(S1))

#     C = fill(d.penalty, N_max, N_max)
#     @views Distances.pairwise!(C, d.ground_dist, S1, S2)

#     # If we have max_dist > 2*ρ we must extend C to allow for interactions in BOTH observations being un-matched.
#     if maximum(view(C, 1:length(S1), 1:length(S2))) > 2*d.penalty
#         C = [
#             C fill(d.penalty/2, N_max, N_min);
#             fill(d.penalty/2, N_min, N_max) fill(0.0, N_min, N_min)
#             ]
#     end
#     return hungarian(C)[2]
# end

function (d::FixPenMatchDist)(X::Nothing, Y::Vector{T})::Float64 where {T}
    return d.penalty * length(Y)
end
function (d::FixPenMatchDist)(X::Vector{T}, Y::Nothing)::Float64 where {T}
    d(Y, X)
end
function (d::FixPenMatchDist)(X::Nothing, Y::Nothing)::Float64 where {T}
    return 0.0
end

struct AvgSizeMatchingDistance{T<:Metric} <: Metric
    ground_dist::T
    penalty::Float64
end

const AvgSizeMatchDist = AvgSizeMatchingDistance
# const AvgSizeMatchDist{T} = AvgSizeMatchingDistance{T}

Base.show(io::IO, d::AvgSizeMatchDist) = print(io, "$(typeof(d))(ρ=$(d.penalty))")

function get_cost_matrix_dynamic(
    d::AvgSizeMatchDist,
    X::Vector{T}, Y::Vector{T}
    ) where {T}

    d_g = d.ground_dist
    N, M = (length(X), length(Y))
    if N < M
        return get_cost_matrix_dynamic(d, Y, X)
    elseif N == M
        return pairwise_inbounds(d_g, X, Y)
    else
        C = pairwise_inbounds(d_g, X, Y)
        mean_size = mean(y -> d_g(y, nothing), Y) # Find mean size in smaller set
        pentalty_vec = [abs(d_g(x, nothing) - mean_size) + d.penalty for x in X]
        size_diff = (N - M)
        return [C [x for x ∈ pentalty_vec, j=1:size_diff]]
    end
end

function get_cost_matrix_fixed(
    d::AvgSizeMatchDist,
    X::Vector{T}, Y::Vector{T}
    ) where {T}

    d_g = d.ground_dist
    N, M = (length(X), length(Y))
    C = pairwise_inbounds(d_g, X, Y)
    if N < M
        mean_size = mean(x->d_g(x,nothing), X)
        penalty_vec = [abs(d_g(y,nothing) - mean_size) + d.penalty for y in Y]
        size_diff = (M - N)
        return [C ; [y for i in 1:size_diff, y ∈ penalty_vec]]
    elseif N > M
        mean_size = mean(y->d_g(y,nothing), Y)
        penalty_vec = [abs(d_g(x,nothing) - mean_size) + d.penalty for x in X]
        size_diff = (N - M)
        return [C [x for x ∈ penalty_vec, i in 1:size_diff]]
    else
        return pairwise_inbounds(d_g, X, Y)
    end

end

function (d::AvgSizeMatchDist)(X::Nothing, Y::Vector{T})::Float64 where {T}
    return (d.penalty * length(Y)) + sum(x->dist.ground_dist(x,nothing), Y)
end
function (d::AvgSizeMatchDist)(X::Vector{T}, Y::Nothing)::Float64 where {T}
    d(Y, X)
end
function (d::AvgSizeMatchDist)(X::Nothing, Y::Nothing)::Float64 where {T}
    return 0.0
end

struct MinDistMatchingDistance{T<:Metric} <: SemiMetric
    ground_dist::T
    penalty::Float64
end

const MinDistMatchDist{T} = MinDistMatchingDistance where {T}

Base.show(io::IO, d::MinDistMatchDist) = print(io, "$(typeof(d))(ρ=$(d.penalty))")

function get_cost_matrix_dynamic(
    d::MinDistMatchDist,
    X::Vector{T}, Y::Vector{T}
    ) where {T}

    d_g = d.ground_dist
    N, M = (length(X), length(Y))
    if N < M
        return get_cost_matrix_dynamic(d, Y, X)
    elseif N == M
        return pairwise_inbounds(d_g, X, Y)
    else
        C = pairwise_inbounds(d_g, X, Y)
        penalty_vec = map(minimum, eachrow(C))
        size_diff = (N - M)
        return [C [x for x ∈ penalty_vec, j=1:size_diff]]
    end
end

function get_cost_matrix_fixed(
    d::MinDistMatchDist,
    X::Vector{T}, Y::Vector{T}
    ) where {T}

    d_g = d.ground_dist
    N, M = (length(X), length(Y))
    C = pairwise_inbounds(d_g, X, Y)
    if N < M
        penalty_vec = map(minimum, eachcol(C))
        size_diff = (M - N)
        return [C ; [y for i in 1:size_diff, y ∈ penalty_vec]]
    elseif N > M
        penalty_vec = map(minimum, eachrow(C))
        size_diff = (N - M)
        return [C [x for x ∈ penalty_vec, i in 1:size_diff]]
    else
        return pairwise_inbounds(d_g, X, Y)
    end

end

# Define type union for all these complete matching distances

const CompleteMatchingDistance = Union{MatchingDistance, FastMatchingDistance, FixedPenaltyMatchingDistance, AvgSizeMatchingDistance, MinDistMatchingDistance}


text:
# Exponential decay at a point
#
# We set up a problem in which a tracer obeys
#
# ∂c/∂t = - c
#
# at a single point. This leads to exponential decay,
#
# c(t) = exp(-t)
#
# in the case that c(0) = 1. We use this problem to test
# the order of convergence of the time-stepper.

module PointExponentialDecay

using Printf
using Statistics
using Oceananigans

include("analysis.jl")

# Simple Exponential decay
c(t) = exp(-t)
Fᶜ = Relaxation(rate=1)

function run_test(; Δt, stop_iteration, timestepper, architecture = CPU())

    grid = RegularRectilinearGrid(size=(1, 1, 1), x=(0, 1), y=(0, 1), z=(0, 1))

    model = NonhydrostaticModel(architecture = architecture,
                                 timestepper = timestepper,
                                        grid = grid,
                                    coriolis = nothing,
                                    buoyancy = nothing,
                                     tracers = :c,
                                     forcing = (c=Fᶜ,))

    set!(model, c=1)
    simulation = Simulation(model, Δt=Δt, stop_iteration=stop_iteration, iteration_interval=stop_iteration)

    @info "Running point exponential decay test for Δt = $Δt, stop_iteration = $stop_iteration..."
    run!(simulation)

    # Calculate errors
    c_simulation = model.tracers.c[1, 1, 1]
    c_analytical = c(model.clock.time)
    L₁_error = abs(c_simulation - c_analytical)

    return (simulation = c_simulation,
            analytical = c_analytical,
                    L₁ = L₁_error)
end

end # module


text: println("Hello World! This is my very first program in one of the most anticipated programming language, Julia ")

max_stars_repo_name: kartike-sood/Hello-World

max_issues_repo_name: kartike-sood/Hello-World

max_forks_repo_name: kartike-sood/Hello-World


text:
using CUDA

export convert_cu, cuzeros, cuones

convert_cu(in_a, X) =  X isa CuArray ? cu(in_a) : in_a
cuzeros(X, args...) = X isa CuArray ? CUDA.zeros(args) : zeros(Float32, args)
cuones(X, args...) = X isa CuArray ? CUDA.ones(args) : ones(Float32, args)


text:
include("../learner/euclidean-1-NN.jl")
include("../algoritmos-busqueda/busqueda-local.jl")
include("../utils/funcion-objetivo.jl")
"""
    ILS_LearnerOneNN(data::Matrix{<:Real}, labels::Vector{<:Real},  iteraciones::Int)
Devuelve un clasificado 1-NN refinado por unos pesos en la distancia euclídea.
Ese vector de pesos se ha obtenido a partir de un algoritmo de búsqueda local
"""
function ILS_LearnerOneNN(data::Matrix{<:Real}, labels, iteraciones::Int=15)
    # Creamos función objetivo
    a = 0.5
    umbral_tasa_reduccion = 0.1
    F = CrearFuncionObjetivo(data, labels, a, umbral_tasa_reduccion )
    w_mejor = 0
    F_w_mejor = -100

    numero_atributos = size(data)[2] # los atributos son las columnas
    # Criterios de acorde a los requisitos
    numero_maximo_vecinos_sin_mejora = 20*numero_atributos
    nuero_maximo_evaluaciones = 1000
    w = 0
    # Buscamos pesos de acorde a parámetros calculados
    for i in 1:iteraciones
        w = PrimeroMejor(numero_maximo_vecinos_sin_mejora,
        nuero_maximo_evaluaciones,numero_atributos,umbral_tasa_reduccion, F, w)
        f_w = F(w)
        if f_w > F_w_mejor
            w_mejor = w
            F_w_mejor = F_w_mejor
        end
        w = GenNeighbourhoodIesimo(
            w ,1, 0.4, numero_atributos, 0.1)
    end
    return WeightedLearnerEuclideanOneNN(w_mejor, data, labels), F_w_mejor, w_mejor
end


text:
# LOB
mutable struct ObjectiveSurface
    param1_name::String
    param1_values::Array{Any, 1}
    param2_name::String
    param2_values::Array{Any, 1}
    surface_points::Int64
    replications::Int64
    parameters::Dict{String, Any}
end
function ObjectiveSurface(param1_name::String, param1_values::Array{Float64, 1},
    param2_name::String, param2_values::Array{Float64, 1},
    surface_points::Int64, parameters::Dict)
    return ObjectiveSurface(param1_name, param1_values, param2_name,
    param2_values, surface_points, parameters["num_paths"], parameters)
end

function (os::ObjectiveSurface)(seed, para=false)
    iterations = os.surface_points^2

    seeds = Int.(rand(MersenneTwister(seed), UInt32, iterations))
    if para==true
        sample_price_paths = SharedArray{Float64,2}((os.parameters["T"]+1, iterations*os.replications))
        @sync @distributed for i in 1:iterations
            os.parameters[os.param1_name] = os.param1_values[i]
            os.parameters[os.param2_name] = os.param2_values[i]
            rdpp = ReactionDiffusionPricePaths(os.parameters)
            index_start = (i-1) * os.replications + 1
            index_end = i * os.replications
            _, _, sample_price_paths[:, index_start:index_end], _, _ = rdpp(seeds[i])
        end
        return sample_price_paths
    else
        sample_price_paths = Array{Float64,2}(undef, os.parameters["T"]+1, iterations*os.replications)
        for i in 1:iterations
            os.parameters[os.param1_name] = os.param1_values[i]
            os.parameters[os.param2_name] = os.param2_values[i]
            rdpp = ReactionDiffusionPricePaths(os.parameters)
            index_start = (i-1) * os.replications + 1
            index_end = i * os.replications
            _, _, sample_price_paths[:, index_start:index_end], _, _ = rdpp(seeds[i])
        end
        return sample_price_paths
    end
end


text:
using Test
using DFTK
using LinearAlgebra

@testset "Hydrogen anion all electron" begin
    function energy_forces(;Ecut, tol)
        lattice = 10 * Matrix{Float64}(I, 3, 3)
        atoms = [ElementCoulomb(:H) => [zeros(3)]]
        model = model_DFT(lattice, atoms, [], n_electrons=2)
        basis = PlaneWaveBasis(model, Ecut, kgrid=(1, 1, 1))

        is_converged = DFTK.ScfConvergenceDensity(tol)
        scfres = self_consistent_field(basis, is_converged=is_converged)
        scfres.energies.total, DFTK.compute_forces(scfres)
    end

    E, forces = energy_forces(Ecut=5, tol=1e-10)
    @test E ≈ -0.28310157350203013
    @test norm(forces) < 1e-9
end


text:
# Issue #15
@testset "anonymous parameter" begin
    f(::Type{T}, n::Int) where T<:Unsigned = rand(T, n)

    patch = @patch f(::Type{UInt8}, n::Int) = collect(UnitRange{UInt8}(1:n))

    apply(patch) do
        @test (@mock f(UInt8, 2)) == [0x01, 0x02]
    end
end


text:
# redundant imports, so this file can be run standalone
using Test
import MD5
using Random

import Nettle

function test_equal_state(md5_value::MD5.MD5_CTX,nettle_value::Nettle.Hasher)
    @test nettle_value.state[1:16] ==
        reinterpret(UInt8, md5_value.state)
end

@testset "Against Nettle initialize, update!, digest!" begin
    for _ in 1:10, chunkcount in [0,1,2,3, 10, 100]
        md5_value = MD5.MD5_CTX()
        nettle_value = Nettle.Hasher("md5")
        test_equal_state(md5_value, nettle_value)
        for _ in 1:chunkcount
            chunksize = rand(0:10000)
            data = rand(UInt8, chunksize)
            MD5.update!(md5_value, data)
            Nettle.update!(nettle_value, data)
            test_equal_state(md5_value,nettle_value)
        end
        @test MD5.digest!(md5_value) == Nettle.digest!(nettle_value)
    end
end

@testset "Against Nettle end to end" begin
    for offset in [0,10^3, 10^4, 10^5]
        iter = offset:(offset+1000)
        for l in iter
            s = randstring(l)
            @test Nettle.hexdigest("md5", s) == bytes2hex(md5(s))
        end
    end
end


text:
# https://github.com/pmodels/mpich/blob/master/src/include/mpi.h.in

const MPI_Aint = Int
const MPI_Count = Int64
const MPI_Offset = Int64

const MPI_Status_size = 20
const MPI_Status_Source_offset = 8
const MPI_Status_Tag_offset = 12
const MPI_Status_Error_offset = 16

const MPI_Info = Cint
const MPI_INFO_NULL = Cint(469762048)

const MPI_Win = Cint
const MPI_WIN_NULL = Cint(536870912)

const MPI_Comm = Cint
const MPI_COMM_NULL = Cint(67108864)
const MPI_COMM_SELF = Cint(1140850689)
const MPI_COMM_WORLD = Cint(1140850688)

const MPI_Request = Cint
const MPI_REQUEST_NULL = Cint(738197504)

const MPI_File = Ptr{Cvoid}
MPI_File_f2c(c::Cint) = ccall((:MPI_File_f2c,libmpi),MPI_File,(Cint,),c)
const MPI_FILE_NULL = Cint(0)

const MPI_Op = Cint
const MPI_OP_NULL = Cint(402653184)
const MPI_BAND = Cint(1476395014)
const MPI_BOR = Cint(1476395016)
const MPI_BXOR = Cint(1476395018)
const MPI_LAND = Cint(1476395013)
const MPI_LOR = Cint(1476395015)
const MPI_LXOR = Cint(1476395017)
const MPI_MAX = Cint(1476395009)
const MPI_MIN = Cint(1476395010)
const MPI_SUM = Cint(1476395011)
const MPI_PROD = Cint(1476395012)
const MPI_MAXLOC = Cint(1476395020)
const MPI_MINLOC = Cint(1476395019)
const MPI_REPLACE = Cint(1476395021)
const MPI_NO_OP = Cint(1476395022)

const MPI_Datatype = Cint
const MPI_DATATYPE_NULL = Cint(201326592)
const MPI_PACKED = Cint(1275068687)
const MPI_CHAR = Cint(1275068673)
const MPI_SIGNED_CHAR = Cint(1275068696)
const MPI_UNSIGNED_CHAR = Cint(1275068674)
const MPI_WCHAR = Cint(1275069454)
const MPI_BYTE = Cint(1275068685)
const MPI_SHORT = Cint(1275068931)
const MPI_UNSIGNED_SHORT = Cint(1275068932)
const MPI_INT = Cint(1275069445)
const MPI_UNSIGNED = Cint(1275069446)
const MPI_LONG = Cint(1275070471)
const MPI_UNSIGNED_LONG = Cint(1275070472)
const MPI_FLOAT = Cint(1275069450)
const MPI_DOUBLE = Cint(1275070475)
const MPI_INT8_T = Cint(1275068727)
const MPI_INT16_T = Cint(1275068984)
const MPI_INT32_T = Cint(1275069497)
const MPI_INT64_T = Cint(1275070522)
const MPI_UINT8_T = Cint(1275068731)
const MPI_UINT16_T = Cint(1275068988)
const MPI_UINT32_T = Cint(1275069501)
const MPI_UINT64_T = Cint(1275070526)
const MPI_C_FLOAT_COMPLEX = Cint(1275070528)
const MPI_C_DOUBLE_COMPLEX = Cint(1275072577)

const MPI_THREAD_SINGLE = Cint(0)
const MPI_THREAD_FUNNELED = Cint(1)
const MPI_THREAD_SERIALIZED = Cint(2)
const MPI_THREAD_MULTIPLE = Cint(3)
const MPI_PROC_NULL = Cint(-1)
const MPI_ANY_SOURCE = Cint(-2)
const MPI_ANY_TAG = Cint(-1)
const MPI_SUCCESS = Cint(0)
const MPI_UNDEFINED = Cint(-32766)
const MPI_LOCK_EXCLUSIVE = Cint(234)
const MPI_LOCK_SHARED = Cint(235)
const MPI_MAX_INFO_KEY = Cint(255)
const MPI_MAX_INFO_VAL = Cint(1024)
const MPI_TAG_UB = Cint(1681915905)
const MPI_COMM_TYPE_SHARED = Cint(1)
const MPI_ORDER_C = Cint(56)
const MPI_ORDER_FORTRAN = Cint(57)
const MPI_UNIVERSE_SIZE = Cint(1681915913)
const MPI_MAX_ERROR_STRING = Cint(512)
const MPI_SUCCESS = Cint(0)
const MPI_MODE_RDONLY = Cint(2)
const MPI_MODE_RDWR = Cint(8)
const MPI_MODE_WRONLY = Cint(4)
const MPI_MODE_CREATE = Cint(1)
const MPI_MODE_EXCL = Cint(64)
const MPI_MODE_DELETE_ON_CLOSE = Cint(16)
const MPI_MODE_UNIQUE_OPEN = Cint(32)
const MPI_MODE_SEQUENTIAL = Cint(256)
const MPI_MODE_APPEND = Cint(128)

const MPI_BOTTOM = reinterpret(SentinelPtr, 0)
const MPI_IN_PLACE = reinterpret(SentinelPtr, -1)
const MPI_STATUS_IGNORE = reinterpret(SentinelPtr, 1)
const MPI_STATUSES_IGNORE = reinterpret(SentinelPtr, 1)

const MPI_IDENT = Cint(0)
const MPI_CONGRUENT = Cint(1)
const MPI_SIMILAR = Cint(2)
const MPI_UNEQUAL = Cint(3)


text:
using RecipesBase

import StaticArrays:
    SArray,
    MArray

import ..Simplices:
    Simplex,
    SSimplex,
    MutableSSimplex,
    connectvertices,
    splitaxes

##################################
# Plotting single triangulations
##################################
@recipe function f(pts::Vector{Vector{T}}, DT::DelaunayTriangulation) where {T}

    for i = 1:length(DT)
        s = Simplex(pts[DT[i]])
        @series begin
            seriestype := :path
            splitaxes(connectvertices(s))
        end
    end
end

@recipe function f(pts::AbstractArray{T, 2}, DT::DelaunayTriangulation) where {T}

    for i = 1:length(DT)
        s = Simplex(pts[:, DT[i]])
        @series begin
            seriestype := :path
            splitaxes(connectvertices(s))
        end
    end
end

@recipe function f(pts::Vector{MArray{Tuple{D},T,1,D}}, DT::DelaunayTriangulation) where {D, T}

    for vertexinds in DT
        s = MutableSSimplex(pts[vertexinds])
        @series begin
            seriestype := :path
            splitaxes(connectvertices(s))
        end
    end
end

@recipe function f(pts::Vector{SArray{Tuple{D},T,1,D}}, DT::DelaunayTriangulation) where {D, T}

    for vertexinds in DT
        s = SSimplex(pts[vertexinds])
        @series begin
            seriestype := :path
            splitaxes(connectvertices(s))
        end
    end
end


text:
# scalar
set = Setup(
    "scalar", # matter
    "advection", # case
    "1d0f0v", # space
    "gks", # flux
    "", # collision: for scalar conservation laws there are none
    1, # species
    1, # interpolation order
    "vanleer", # limiter
    "period", # boundary
    0.5, # cfl
    1.0, # simulation time
)
pSpace = PSpace1D(0.0, 1.0, 100, 1)
vSpace = nothing
property = Scalar(1.0, 1e-6)
w0 = 1.0
prim0 = conserve_prim(w0, property.a)
ib = IB(x -> sin(2π * x), property)
ks = SolverSet(set, pSpace, vSpace, property, ib)

ctr, face = init_fvm(ks, ks.ps)
for i in eachindex(ctr)
    ctr[i].w = sin(2π * ks.pSpace.x[i])
end

dt = 1e-3
reconstruct!(ks, ctr)
evolve!(ks, ctr, face, dt)
update!(ks, ctr, face, dt, 0.0)


text:
using Test, CSV, Dates, Tables, WeakRefStrings, CategoricalArrays, DataFrames

const dir = joinpath(dirname(pathof(CSV)), "../test/testfiles")

@eval macro $(:try)(ex)
    quote
        try $(esc(ex))
        catch
        end
    end
end

@testset "CSV" begin

include("files.jl")
include("iteration.jl")
include("write.jl")

@testset "CSV.validate" begin
    io = IOBuffer("""A,B,C
    1,1,10
    6,1""")

    @test_throws CSV.ExpectedMoreColumnsError CSV.validate(io)

    io = IOBuffer("""A;B;C
    1,1,10
    2,0,16""")
    @test_throws CSV.TooManyColumnsError CSV.validate(io)

    io = IOBuffer("""A;B;C
    1,1,10
    2,0,16""")
    @test_throws CSV.ExpectedMoreColumnsError CSV.validate(io; delim=';')

    io = IOBuffer("""a b c d e
    1 2  3 4 5
    1 2 3  4 5
    1  2 3  4 5""")
    @test_throws CSV.TooManyColumnsError CSV.validate(io; delim=' ')

    # 323
    CSV.validate(IOBuffer("a0001000\na0001000"), datarow=1, categorical=true)
end

@testset "transform" begin
    csv = """A,B,C
    1,1.0,"hey"
    ,2.0,"there"
    3,3.0,"sailor"
    """

    df = CSV.read(IOBuffer(csv), transforms=Dict{String, Base.Callable}("C" => Symbol))
    @test df.C == [:hey, :there, :sailor]

    df = CSV.read(IOBuffer(csv); transforms=Dict("A"=>x->x+1))
    @test isequal(df.A, [2, missing, 4])
    @test typeof(df.A) == Vector{Union{Missing, Int64}}
    @test size(df) == (3, 3)

    df = CSV.read(IOBuffer(csv); transforms=Dict("A"=>x->coalesce(x+1, 0)))
    @test df.A == [2, 0, 4]
    @test eltype(df.A) <: Signed

    df = CSV.read(IOBuffer(csv); transforms=Dict("A"=>x->coalesce(x+1, 0.0)))
    @test df.A == [2, 0.0, 4]
    @test eltype(df.A) <: Real

    df = CSV.read(IOBuffer(csv); transforms=Dict(2=>x->x==2.0 ? missing : x))
    @test size(df) == (3, 3)
    @test isequal(df.B, [1.0, missing, 3.0])
    @test typeof(df.B) == Vector{Union{Float64, Missing}}

    f = CSV.File(IOBuffer(csv); allowmissing=:auto)

    # transforms on CSV.Columns
    df = f |> transform(A=x->x+1) |> DataFrame
    @test size(df) == (3, 3)
    df = f |> transform(A=x->x+1) |> columntable
    @test length(df) == 3
    df = f |> transform(A=x->x+1) |> rowtable
    @test length(df) == 3

    # transforms on CSV.RowIterator
    f = CSV.File(joinpath(dir, "pandas_zeros.csv"), allowmissing=:none);
    df = f |> transform(_0=x->x+1) |> DataFrame;
    @test size(df) == (100000, 50)
    df = f |> transform(_0=x->x+1) |> columntable;
    @test length(df) == 50
    df = f |> transform(_0=x->x+1) |> rowtable;
    @test length(df) == 100000
end

@testset "CategoricalArray levels (including ordering)" begin
    f = CSV.File(IOBuffer("X\nb\nc\na\nc"), types=[CategoricalString{UInt32}])
    v = iterate(f, 1)[1].X
    @test v == "b"
    @test levels(v.pool) == ["b"]
    v = iterate(f, 2)[1].X
    @test v == "c"
    @test levels(v.pool) == ["b", "c"]
    v = iterate(f, 3)[1].X
    @test v == "a"
    @test levels(v.pool) == ["a", "b", "c"]
    v = iterate(f, 4)[1].X
    @test v == "c"
    @test levels(v.pool) == ["a", "b", "c"]
end

include("deprecated.jl")

end


text:
import Base.show

export Landscape, prob_neutral

@doc """ function makelinks(n::Int64, k::Int64, near::Bool)

Return an array of links (bit masks) that select the linked loci.
Thus, link[i] selects the loci that are linked to locus  `i`.
Each link contains `k + 1` one bits.  link[i] is guaranteed to have bit i
(counting from the right) set to one.  The positions of the remaining `k`
one bits are randomly selected from the set {1:n} \ {i} (where \ denotes
set difference).
"""
function makelinks(n::Int64, k::Int64, near::Bool)
  if k >= n
    error("k must be strictly less than n")
  end
  links = zeros(AlleleMask, n)
  if near
    # TODO: See if we can come up with a prettier algorithm here.
    if isodd(k)
      error("Periodic boundary requires an even-valued K.")
    end
    offset = div(k, 2)
    for i = 1:n
      mask = AlleleMask(0)
      for j = (i - offset):(i + offset)
        shift = if j < 1
          j + n - 1
        elseif j > n
          j - n - 1
        else
          j - 1
        end
        mask = mask | (AlleleMask(1) << shift)
      end
      links[i] = mask
    end
  else
    for i = 1:n
      linkedloci = sample(setdiff(0:(n - 1), [i - 1]), k, replace=false)
      mask = reduce((a, b) -> a | (AlleleMask(1) << b),
          AlleleMask(1) << (i - 1), linkedloci)
      links[i] = mask
    end
  end
  return links
end

# TODO: Provide a function to pre-fill contribs (for use with small k).

abstract type Landscape end

show(io::Base.IO, l::Landscape) = print(io, "$(typeof(l))($(l.n), $(l.k))")

include("landscape/nk.jl")
include("landscape/nkq.jl")
include("landscape/nkp.jl")


text:
# --------------------------------------------------------------------------
# ACE.jl and SHIPs.jl: Julia implementation of the Atomic Cluster Expansion
# Copyright (c) 2019 Christoph Ortner <christophortner0@gmail.com>
# Licensed under ASL - see ASL.md for terms and conditions.
# --------------------------------------------------------------------------

using Test
using ACE, JuLIP, JuLIP.Testing, QuadGK, LinearAlgebra, ACE.JacobiPolys,
      BenchmarkTools
using ACE: TransformedJacobi, transform, transform_d,
             alloc_B, alloc_temp, alloc_temp_d, alloc_dB, IntS

using JuLIP: evaluate!, evaluate_d!

function randR()
   R = rand(JVecF) .- 0.5
   return (0.9 + 2 * rand()) * R/norm(R)
end
randR(N) = [ randR() for n=1:N ], zeros(Int16, N), 0
randcoeffs(B) = rand(length(B)) .* (1:length(B)).^(-2)

##

Nmax = 4
rl, ru = 0.5, 3.0
fcut =  PolyCutoff2s(2, rl, ru)
trans = PolyTransform(2, 1.0)
spec = SparseSHIP(Nmax, 10)
aceB = SHIPBasis(spec, trans, fcut)

bgrp = aceB.bgrps[1]

##

Nr = 50
Rs, Zs = randR(Nr)
tmp = alloc_temp(aceB, Nr)
B = ACE.alloc_B(aceB)
tmp2 = alloc_temp(aceB2, Nr)
B2 = ACE.alloc_B(aceB2)

evaluate!(B, tmp, aceB, Rs, Zs, 0)
evaluate!(B2, tmp2, aceB2, Rs, Zs, 0)
@show B ≈ B2

@btime evaluate!($B, $tmp, $aceB, $Rs, $Zs, 0);
@btime evaluate!($B2, $tmp2, $aceB2, $Rs, $Zs, 0);

##

tmpd = alloc_temp_d(aceB, Nr)
dB = alloc_dB(aceB, Nr)
evaluate_d!(B, dB, tmpd, aceB, Rs, Zs, 0)

tmpd2 = alloc_temp_d(aceB2, Nr)
dB2 = alloc_dB(aceB2, Nr)
evaluate_d!(dB2, tmpd2, aceB2, Rs, Zs, 0)

@show dB2 ≈ dB

@btime evaluate_d!($B, $dB, $tmpd, $aceB, $Rs, $Zs, 0)
@btime evaluate_d!($dB2, $tmpd2, $aceB2, $Rs, $Zs, 0)

##

# ##
#
# f = let dB2 = dB2, tmpd2 = tmpd2, aceB2 = aceB2, Rs = Rs, Zs = Zs
#    () -> evaluate_d!(dB2, tmpd2, aceB2, Rs, Zs, 0)
# end
# runn(f, N) = (for n=1:N; f(); end)
# runn(f, 1)
# ##
#
# using Profile
# Profile.clear()
# @profile runn(f, 100)
# Profile.print()
#
# ##
#
#
# # debugging -> shows that ACE.grad_AA_Rj! is correct...
#
# _AA(Rs) = (
#    ACE.precompute_A!(tmp2, aceB2, Rs, Zs, 1);
#    ACE.precompute_AA!(tmp2, aceB2, 1);
#    return copy(tmp2.AA[1])
#    )
#
# _dAA(Rs, j) = (
#    ACE.precompute_dA!(tmpd2, aceB2, Rs, Zs, 1);
#    ACE.precompute_AA!(tmpd2, aceB2, 1);
#    ACE.grad_AA_Rj!(tmpd2, aceB2, j, Rs, Zs, 1);
#    return copy(tmpd2.dAAj[1])
#    )
#
#
# u1 = rand(JVecF); u1 /= norm(u1)
# jj = 4
#
# for p = 2:12
#    h = 0.1^p
#
#    A = _AA(Rs)
#    Rs_h = copy(Rs); Rs_h[jj] += h * u1
#    A_h = _AA(Rs_h)
#    dA_h = (A_h - A) / h
#
#    dAA = _dAA(Rs, jj)
#    dAA_u1 = [ dot(u1, da) for da in dAA ]
#
#    @show norm(real.(dAA_u1 - dA_h), Inf)
# end
#
#

# h = 1e-5
#
# A = _AA(Rs)
# Rs_h = copy(Rs); Rs_h[1] += h * u1
# A_h = _AA(Rs_h)
# dA_h = (A_h - A) / h
#
# dAA = _dAA1(Rs)
# dAA_u1 = [ dot(da, u1) for da in dAA ]
#
# @show norm(real.(dAA_u1 - dA_h), Inf)
#
# @show imag.(dAA_u1 - dA_h)
#
#
#
# _A(Rs) = (
#    ACE.precompute_A!(tmp2, aceB2, Rs, Zs, 1);
#    return copy(tmp2.A[1])
#    )
#
# _dA1(Rs) = begin
#       ACE.precompute_dA!(tmpd2, aceB2, Rs, Zs, 1);
#       alist = aceB2.alists[1]
#       dA = zeros(JVec{ComplexF64}, length(alist))
#       for n = 1:length(alist)
#          zklm = alist[n]
#          dA[n] = ACE.grad_phi_Rj(Rs[1], 1, zklm, tmpd2)
#       end
#       return dA
#    end
#
# _dA1(Rs)
#
# for p = 2:12
#    h = 0.1^p
#
#    A = _A(Rs)
#    Rs_h = copy(Rs); Rs_h[1] += h * u1
#    A_h = _A(Rs_h)
#    dA_h = (A_h - A) / h
#
#    dA = _dA1(Rs)
#    dA_u1 = [ dot(u1, da) for da in dA ]
#
#    @show norm(imag.(dA_u1 - dA_h), Inf)
# end
#
#
# R = Rs[1]
# _phi(R) = _A([R])
# _dphi(R) = begin
#       ACE.precompute_dA!(tmpd2, aceB2, [R], [0], 1);
#       alist = aceB2.alists[1]
#       dphi = zeros(JVec{ComplexF64}, length(alist))
#       for n = 1:length(alist)
#          zklm = alist[n]
#          dphi[n] = ACE.grad_phi_Rj(R, 1, zklm, tmpd2)
#       end
#       return dphi
#    end
#
# for p = 2:12
#    h = 0.1^p
#
#    phi = _phi(R)
#    phi_h = _phi(R + h * u1)
#    dphi_h = (phi_h - phi) / h
#
#    dphi = dot.(Ref(u1), _dphi(R))
#
#    @show norm(imag.(dphi - dphi_h), Inf)
# end
#
# _phi(R)
# dot.(_dphi(R), Ref(u1))
#
#
#
# tmpd.J ≈ tmpd2.JJ
# tmpd.dJ ≈ tmpd2.dJJ
# tmpd.Y ≈ tmpd2.YY
# tmpd.dY ≈ tmpd2.dYY
#
#
# SH = aceB2.SH
# Y = copy(tmpd2.Y)
# dY = copy(tmpd2.dY)
#
# _Y(R) = evaluate(SH, R)
# _dY(R, u) = dot.(Ref(u), evaluate_d(SH, R)[2])
#
# h = 1e-5
# dY0 = _dY(R, u1)
# dYh = (_Y(R+h*u1) - _Y(R)) / h
# norm(imag.(dY0 - dYh), Inf)
# imag.(dY0 - dYh)


text:
using ThreadsXBenchmarks
using Test

@testset "ThreadsXBenchmarks.jl" begin
    # Write your own tests here.
end


text:
using DrWatson
quickactivate(@__DIR__)
@show projectdir()
##
using StatsBase
using StatsPlots
using DataFrames
using CSV
using CovidSurvey
##
const mobility_names = [
    "retail_and_recreation_percent_change_from_baseline",
    "grocery_and_pharmacy_percent_change_from_baseline",
    "workplaces_percent_change_from_baseline",
    "residential_percent_change_from_baseline"
]

const new_mobility_names = [
    "retail",
    "grocery",
    "workplaces",
    "residential"
]
## ============================================================================
#                    national
#  ============================================================================
fname = "/home/and/data/google-mobility/Global_Mobility_Report.csv"
mobil_raw = CSV.File(fname) |> DataFrame

function selectdenmark(mobil_raw)
    mobil_raw = filter(:country_region => ==("Denmark"), mobil_raw)
    mobil_dk = filter(:sub_region_1 => ismissing, mobil_raw)
    mobil_dk[!, "sub_region_1"] = replace(mobil_dk.sub_region_1,
        "Capital Region of Denmark" => "capital"
        ,"Central Denmark Region" => "central"
        ,"North Denmark Region" => "north"
        ,"Region of Southern Denmark" => "south"
        ,"Region Zealand" => "zealand"
    )
    @assert issorted(mobil_dk.date)
    mobil_dk = select(mobil_dk, :date, mobility_names .=> new_mobility_names)
    mobil_dk[!,:nonresidential] = [ mean(r) for r in eachrow( select(mobil_dk, Not([:date, :residential])) ) ]
    disallowmissing!(mobil_dk)
    return mobil_dk
end
mobil_dk = selectdenmark(mobil_raw)
CSV.write(projectdir("data/google/", "mobility_denmark.csv"), mobil_dk)

## ============================================================================
#                    regional
#  ============================================================================
# select all regions
fname = "/home/and/data/google-mobility/Global_Mobility_Report.csv"
mobil_raw = CSV.File(fname) |> DataFrame
function selectregions(mobil_raw)
    mobil_raw = filter(:country_region => ==("Denmark"), mobil_raw)
    mobil_regions = filter(:sub_region_1 => x->!ismissing(x), mobil_raw)
    mobil_regions[!, "sub_region_1"] = replace(mobil_regions.sub_region_1,
        "Capital Region of Denmark" => "capital"
        ,"Central Denmark Region" => "central"
        ,"North Denmark Region" => "north"
        ,"Region of Southern Denmark" => "south"
        ,"Region Zealand" => "zealand"
    )
    return mobil_regions
end
mobil_regions = selectregions(mobil_raw)
##
function selectregion(mobil_regions, region)
    mobil_region = filter(:sub_region_1 => ==(region), mobil_regions)
    mobil_region = filter(:sub_region_2 => ismissing, mobil_region)
    @assert issorted(mobil_region.date)
    mobil_region = select(mobil_region, :date, mobility_names .=> new_mobility_names)
    mobil_region[!,:nonresidential] = [ mean(r) for r in eachrow( select(mobil_region, Not([:date, :residential])) ) ]
    disallowmissing!(mobil_region)
    return mobil_region
end

for region in Regional.regions
    mobil_region = selectregion(mobil_regions, region)
    CSV.write(projectdir("data/google/", "mobility_region=$region.csv"), mobil_region)
end

##           OLD
# function select_region(df, r)
#     mask = df.sub_region_1 .== r
#     mask = DataFrames.coalesce.(mask, Ref(false))
#     df = df[mask, :]
#     df = df[.! ismissing.(df.sub_region_2), :]
#     df.sub_region_2 = standardize_name.(df.sub_region_2)
#
#     # @assert any(.! ismissing.(Array(mobil[:,mobility_names])))
#     # for n in mobility_names
#     #     df[!,n] = convert(Vector{Float64}, df[:,n])
#     # end
#     return df
# end
#
# function weightedmean(xs, ws)
#     nom = 0.
#     denom = 0.
#     for (x,w) in zip(xs, ws)
#         if !ismissing(x)
#             nom += x*w
#             denom += w
#         end
#     end
#     return nom/denom
# end
#
# function cleanup(mobil)
#     gp = groupby(mobil, "date", sort=true)
#     mobil = combine(gp, [[n,"population"] for n in mobility_names] .=> weightedmean, renamecols=false)
#     mobil[!,"nonresidential"] = mean(eachcol(select(mobil, Not(["residential_percent_change_from_baseline_population", "date"]))))
#     rename!(mobil, [n*"_population" for n in mobility_names] .=> new_mobility_names)
#     @assert issorted(mobil.date)
#     return mobil
# end
#
# function readpopulation()
#     fname = projectdir("data/population/Denmark_2021Q2.csv")
#     pop = CSV.File(fname; header=false, select=[4,5]) |> DataFrame
#     rename!(pop, "Column4"=>"name", "Column5"=>"count")
#     pop.name = standardize_name.(pop.name)
#     pop
# end
#
# const population = readpopulation()
#
# function get_population(s)
#     if s == "Brondby"
#         s = "Br_ndby"
#     end
#     i = findfirst(x->contains(x,s), population.name)
#     isnothing(i) && (@error s; return missing)
#     return population.count[i]
# end
#
# function replace_special_characters(s)
#     s = replace(s, '�'=>'_')
#     s = replace(s, 'ø'=>'_')
#     s = replace(s, 'æ'=>'_')
#     s = replace(s, 'å'=>'_')
#     s = replace(s, '\ue6'=>'_')
#     s = replace(s, '\uf8'=>'_')
#     s = replace(s, '\ue5'=>'_')
#     return s
# end
#
# function standardize_name(s_raw)
#     s = if split(s_raw, ' ')|>last == "Municipality"
#             first(split(s_raw, ' '))
#         else
#             s_raw
#         end
#     replace_special_characters(s)
# end
#
# for region in Regional.regions
#     mobil = select_region(mobil_raw, region)
#     mobil[!,:population] = get_population.(mobil.sub_region_2)
#     @assert any(.! ismissing.(mobil.population))
#
#     mobil = cleanup(mobil)
#     CSV.write(projectdir("data/google/", "mobility_region=$region.csv"), mobil)
# end

max_stars_repo_name: andreaskoher/CovidSurvey

max_issues_repo_name: andreaskoher/CovidSurvey

max_forks_repo_name: andreaskoher/CovidSurvey


text:
module FileFormats

include("csv.jl")
include("json.jl")

end


text:
# # Detailed Look
#
#md # [![](https://mybinder.org/badge_logo.svg)](@__BINDER_ROOT_URL__/notebooks/detailed_look.ipynb)
#md # [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/notebooks/detailed_look.ipynb)
#
# A more detailed look at spatial interpolation, integration through time, and I/O.
# For additional documentation e.g. see
# [1](https://JuliaClimate.github.io/IndividualDisplacements.jl/dev/),
# [2](https://JuliaClimate.github.io/MeshArrays.jl/dev/),
# [3](https://docs.juliadiffeq.org/latest/solvers/ode_solve.html),
# [4](https://en.wikipedia.org/wiki/Displacement_(vector)).
# Here we illustrate a few things in more detail:
#
# - reading velocities from file.
#   - gridded velocity output (U*data, V*data)
#   - pre-computed trajectory output (`float_traj*data`)
# - interpolating `U,V` from gridded output to individual locations
#   - compared with `u,v` from `float_traj*data`
# - computing trajectories (location v time) using `OrdinaryDiffEq.jl`
#   - compared with `x(t),y(t)` from `float_traj*data`

# ## 1. Import Software

using IndividualDisplacements, OrdinaryDiffEq, DataFrames, MITgcmTools
p=dirname(pathof(IndividualDisplacements))
include(joinpath(p,"../examples/recipes_plots.jl"))
include(joinpath(p,"../examples/example123.jl"))
include(joinpath(p,"../examples/helper_functions.jl"));

# ## 2. Read Trajectory Output
#
# from `MITgcm/pkg/flt`

dirIn=IndividualDisplacements.flt_example
prec=Float32
df=read_flt(dirIn,prec);

#md plt=plot_paths(df,300,100000.0)

# ## 3. Read Gridded Variables
#
# using `MeshArrays.jl` and e.g. a NamedTyple

𝑃,Γ=example2_setup();

# ## 4. Visualize Velocity Fields
#
# ```
# plt=heatmap(Γ["mskW"][1,1].*𝑃.u0,title="U at the start")
# plt=heatmap(Γ["mskW"][1,1].*𝑃.u1-𝑃.u0,title="U end - U start")
# ```

# ## 5. Visualize Trajectories
#
# (select one trajectory)

tmp=df[df.ID .== 200, :]
tmp[1:4,:]

# Super-impose trajectory over velocity field (first for u ...)

x=Γ["XG"].f[1][:,1]
y=Γ["YC"].f[1][1,:]
z=transpose(Γ["mskW"][1].*𝑃.u0);

# plt=contourf(x,y,z,c=:delta)
# plot!(tmp[:,:lon],tmp[:,:lat],c=:red,w=4,leg=false)

# Super-impose trajectory over velocity field (... then for v)

x=Γ["XC"].f[1][:,1]
y=Γ["YG"].f[1][1,:]
z=transpose(Γ["mskW"][1].*𝑃.v0);

# plt=contourf(x,y,z,c=:delta)
# plot!(tmp[:,:lon],tmp[:,:lat],c=:red,w=4,leg=false)

# ## 6. Interpolate Velocities

dx=Γ["dx"]
uInit=[tmp[1,:lon];tmp[1,:lat]]./dx
nSteps=Int32(tmp[end,:time]/3600)-2
du=fill(0.0,2);

# Visualize and compare with actual grid point values -- jumps on the tangential component are expected with linear scheme:

tmpu=fill(0.0,100)
tmpv=fill(0.0,100)
tmpx=fill(0.0,100)
for i=1:100
    tmpx[i]=500.0 *i./dx
    dxy_dt(du,[tmpx[i];0.499./dx],𝑃,0.0)
    tmpu[i]=du[1]
    tmpv[i]=du[2]
end

#md plt=plot(tmpx,tmpu,label="u (interp)")
#md plot!(Γ["XG"].f[1][1:10,1]./dx,𝑃.u0[1:10,1],marker=:o,label="u (C-grid)")
#md plot!(tmpx,tmpv,label="v (interp)")
#md plot!(Γ["XG"].f[1][1:10,1]./dx,𝑃.v0[1:10,1],marker=:o,label="v (C-grid)")

# And similarly in the other direction

tmpu=fill(0.0,100)
tmpv=fill(0.0,100)
tmpy=fill(0.0,100)
for i=1:100
    tmpy[i]=500.0 *i./dx
    dxy_dt(du,[0.499./dx;tmpy[i]],𝑃,0.0)
    tmpu[i]=du[1]
    tmpv[i]=du[2]
end

# plt=plot(tmpx,tmpu,label="u (interp)")
# plot!(Γ["YG"].f[1][1,1:10]./dx,𝑃.u0[1,1:10],marker=:o,label="u (C-grid)")
# plot!(tmpx,tmpv,label="v (interp)")
# plot!(Γ["YG"].f[1][1,1:10]./dx,𝑃.v0[1,1:10],marker=:o,label="v (C-grid)")

# Compare recomputed velocities with those from `pkg/flt`

nSteps=2998
tmpu=fill(0.0,nSteps); tmpv=fill(0.0,nSteps);
tmpx=fill(0.0,nSteps); tmpy=fill(0.0,nSteps);
refu=fill(0.0,nSteps); refv=fill(0.0,nSteps);
for i=1:nSteps
    dxy_dt_replay(du,[tmp[i,:lon],tmp[i,:lat]],tmp,tmp[i,:time])
    refu[i]=du[1]./dx
    refv[i]=du[2]./dx
    dxy_dt(du,[tmp[i,:lon],tmp[i,:lat]]./dx,𝑃,Float64(tmp[i,:time]))
    tmpu[i]=du[1]
    tmpv[i]=du[2]
end

# plt=plot(tmpu,label="u")
# plot!(tmpv,label="v")
# plot!(refu,label="u (ref)")
# plot!(refv,label="v (ref)")

# ## 6. Compute Trajectories
#
# Solve through time using `OrdinaryDiffEq.jl` with
#
# - `dxy_dt` is the function computing `dxy/dt`
# - `uInit` is the initial condition `u @ tspan[1]`
# - `tspan` is the time interval
# - `uvetc` are parameters for `dxy_dt`
# - `Tsit5` is the time-stepping scheme
# - `reltol` and `abstol` are tolerance parameters

tspan = (0.0,nSteps*3600.0)
#prob = ODEProblem(dxy_dt_replay,uInit,tspan,tmp)
prob = ODEProblem(dxy_dt,uInit,tspan,𝑃)
sol = solve(prob,Tsit5(),reltol=1e-8,abstol=1e-8)
sol[1:4]

# Compare recomputed trajectories with originals from `MITgcm/pkg/flt`

ref=transpose([tmp[1:nSteps,:lon] tmp[1:nSteps,:lat]])
maxLon=80*5.e3
maxLat=42*5.e3
#show(size(ref))
for i=1:nSteps-1
    ref[1,i+1]-ref[1,i]>maxLon/2 ? ref[1,i+1:end]-=fill(maxLon,(nSteps-i)) : nothing
    ref[1,i+1]-ref[1,i]<-maxLon/2 ? ref[1,i+1:end]+=fill(maxLon,(nSteps-i)) : nothing
    ref[2,i+1]-ref[2,i]>maxLat/2 ? ref[2,i+1:end]-=fill(maxLat,(nSteps-i)) : nothing
    ref[2,i+1]-ref[2,i]<-maxLat/2 ? ref[2,i+1:end]+=fill(maxLat,(nSteps-i)) : nothing
end
ref=ref./dx;

#md plt=plot(sol[1,:],sol[2,:],linewidth=5,title="Using Recomputed Velocities",
#md      xaxis="lon",yaxis="lat",label="Julia Solution") # legend=false
#md plot!(ref[1,:],ref[2,:],lw=3,ls=:dash,label="MITgcm Solution")


text:
using DiffEqFlux
using OrdinaryDiffEq
using Optimisers
using Flux
using CUDA

CUDA.allowscalar(false)

using PyCall
using PyCallChainRules.Torch: TorchModuleWrapper, torch

u0 = Float32[2.; 0.]
datasize = 30
tspan = (0.0f0, 1.5f0)

function trueODEfunc(du,u,p,t)
    true_A = [-0.1 2.0; -2.0 -0.1]
    du .= ((u.^3)'true_A)'
end
t = range(tspan[1],tspan[2],length=datasize)
prob = ODEProblem(trueODEfunc,u0,tspan)
ode_data = (Array(solve(prob,Tsit5(),saveat=t)))
ode_data = CUDA.cu(ode_data)

torch_module = torch.nn.Sequential(
    torch.nn.Linear(2, 50), torch.nn.Tanh(),
    torch.nn.Linear(50, 2)
)

torch_module = torch_module.to(device=torch.device("cuda"))

jlmod = TorchModuleWrapper(torch_module)
p, re = Optimisers.destructure(jlmod)

u0 = CUDA.cu(u0)

dudt(u, p, t) = re(p)(u)
prob = ODEProblem(dudt, u0, tspan)

function predict_n_ode_gpu(p)
    Flux.gpu(solve(prob,Tsit5(),u0=u0,p=p,saveat=t))
end

function loss_n_ode_gpu(p)
    pred = predict_n_ode_gpu(p)
    loss = sum(abs2,ode_data .- pred)
    loss
end

@info "before" loss_n_ode_gpu(p)

function train(p;nsteps=1000)
    opt = Optimisers.ADAM(0.01)
    state = Optimisers.setup(opt, p)

    for i in 1:nsteps
        gs, = Flux.gradient(p) do ps
            loss_n_ode_gpu(ps)
        end
        state, p = Optimisers.update(state, p, gs)
    end
    return p
end

newp = train(p)

@info "after" loss_n_ode_gpu(newp)


text:
# Hypre.jl: Wrapper around the Hypre library (http://acts.nersc.gov/hypre/)
VERSION >= v"0.4.0-dev+6521" && __precompile__()
module Hypre
const depsfile = joinpath(dirname(dirname(@__FILE__)), "deps", "deps.jl")
if isfile(depsfile)
    include(depsfile)
else
    error("Hypre not properly installed. Please run Pkg.build(\"Hypre\")")
end

using MPI

include("types_and_consts.jl")
include("HYPRE_sstruct_ls.jl")
include("HYPRE_struct_ls.jl")
include("HYPRE_parcsr_ls.jl")
include("HYPRE_DistributedMatrixPilutSolver_protos.jl")
include("HYPRE_lobpcg.jl")

end # module