text:
using Test
using CUDAnative
using CuArrays

using CthulhuVision.Random
using CthulhuVision.Rendering: makeprng

function splitmix64_gpu(seed::UInt64, n::Int, result)

    x = threadIdx().x

    if x == 1
        sm = SplitMix64(seed)
        for i = 1:n
            @inbounds result[i] = next(sm)
        end
    end

    nothing
end

function xoshiro256pp_gpu(s0::UInt64, s1::UInt64, s2::UInt64, s3::UInt64, n::Int, result)

    x = threadIdx().x

    if x == 1
        r = Xoshiro256pp(s0, s1, s2, s3)
        for i = 1:n
            @inbounds result[i] = next(r)
        end
    end

    nothing
end

function uniform32_gpu(us, fs, n::Int)

    x = threadIdx().x

    if x == 1
        for i = 1:n
            @inbounds fs[i] = uniform32(us[i])
        end
    end

    nothing
end

# This kernel tests the factory function that creates a PRNG, with a seed value
# generated from the CUDA thread and block indices.
# Each thread generates a single Float32 value from the PRNG, and stores it in
# the `values` array. This allows us to verify that all Float32 values are different,
# so that we don't accidentally use a single seed, or a seed that repeats.
function prng_from_thread_and_block_gpu(values, n)
    idx = threadIdx().x + (blockIdx().x - 1) * blockDim().x
    if idx <= n
        rng = makeprng()
        @inbounds values[idx] = next(rng)
    end
    nothing
end

# These test values were generated from the C code available at
# http://prng.di.unimi.it/splitmix64.c
# 2020-04-16

@testset "CthulhuVision.Random GPU" begin
    @testset "uniform32 is 0.0f0 <= x < 1.0f0" begin
        n = 1000
        us = CuArray{UInt64}(rand(UInt64, n))
        fs = CuArray{Float32}(undef, n)

        CUDAnative.@sync begin
            @cuda threads=16 uniform32_gpu(us, fs, n)
        end

        floats_host = Vector{Float32}(fs)

        for k in floats_host
            @test 0.0f0 <= k < 1.0f0
        end
    end

    @testset "SplitMix64 validation: Seed 0" begin
        n = 100
        result = CuArray{UInt64}(undef, n)

        CUDAnative.@sync begin
            @cuda threads=16 splitmix64_gpu(UInt64(0), n, result)
        end

        # Seed 0x0000000000000000:
        fromseed0 = [
            0xe220a8397b1dcdaf, 0x6e789e6aa1b965f4, 0x06c45d188009454f, 0xf88bb8a8724c81ec, 0x1b39896a51a8749b,
            0x53cb9f0c747ea2ea, 0x2c829abe1f4532e1, 0xc584133ac916ab3c, 0x3ee5789041c98ac3, 0xf3b8488c368cb0a6,
            0x657eecdd3cb13d09, 0xc2d326e0055bdef6, 0x8621a03fe0bbdb7b, 0x8e1f7555983aa92f, 0xb54e0f1600cc4d19,
            0x84bb3f97971d80ab, 0x7d29825c75521255, 0xc3cf17102b7f7f86, 0x3466e9a083914f64, 0xd81a8d2b5a4485ac,
            0xdb01602b100b9ed7, 0xa9038a921825f10d, 0xedf5f1d90dca2f6a, 0x54496ad67bd2634c, 0xdd7c01d4f5407269,
            0x935e82f1db4c4f7b, 0x69b82ebc92233300, 0x40d29eb57de1d510, 0xa2f09dabb45c6316, 0xee521d7a0f4d3872,
            0xf16952ee72f3454f, 0x377d35dea8e40225, 0x0c7de8064963bab0, 0x05582d37111ac529, 0xd254741f599dc6f7,
            0x69630f7593d108c3, 0x417ef96181daa383, 0x3c3c41a3b43343a1, 0x6e19905dcbe531df, 0x4fa9fa7324851729,
            0x84eb4454a792922a, 0x134f7096918175ce, 0x07dc930b302278a8, 0x12c015a97019e937, 0xcc06c31652ebf438,
            0xecee65630a691e37, 0x3e84ecb1763e79ad, 0x690ed476743aae49, 0x774615d7b1a1f2e1, 0x22b353f04f4f52da,
            0xe3ddd86ba71a5eb1, 0xdf268adeb6513356, 0x2098eb73d4367d77, 0x03d6845323ce3c71, 0xc952c5620043c714,
            0x9b196bca844f1705, 0x30260345dd9e0ec1, 0xcf448a5882bb9698, 0xf4a578dccbc87656, 0xbfdeaed9a17b3c8f,
            0xed79402d1d5c5d7b, 0x55f070ab1cbbf170, 0x3e00a34929a88f1d, 0xe255b237b8bb18fb, 0x2a7b67af6c6ad50e,
            0x466d5e7f3e46f143, 0x42375cb399a4fc72, 0x8c8a1f148a8bb259, 0x32fcab5daed5bdfc, 0x9e60398c8d8553c0,
            0xee89cceb8c4064c0, 0xdb0215941d86a66f, 0x5ccde78203c367a8, 0xf1bcbc6a1ec11786, 0xef054fceee954551,
            0xdf82012d0555c6df, 0x292566ff72403c08, 0xc4dd302a1bfa1137, 0xd85f219db5c554e1, 0x6a27ff807441bcd2,
            0x96a573e9b48216e8, 0x46a9fdac40bf0048, 0x3dd12464a0ee15b4, 0x451e521296a7eea1, 0x56e4398a98f8a0fd,
            0x7b7dc2160e3335a7, 0xc679ee0bebcb1cca, 0x928d6f2d7453424e, 0x1b38994205234c6d, 0x8086d193a6f2b568,
            0x21c6e26639ac2c65, 0xd9dccac414d23c6f, 0x91cd642057e00235, 0x77fc607dc6589373, 0x05b8abe26dd3aee7,
            0x12f6436ac376cc66, 0x64952424897b2307, 0xee8c2baf6343e5c3, 0xdc4c613d9eba2304, 0x3505b7796bd1a506,
        ]

        result_host = Vector{UInt64}(result)

        for i = 1:n
            @test result_host[i] == fromseed0[i]
        end
    end

    @testset "SplitMix64 validation: Seed 0x1234567890123456" begin
        n = 100
        result = CuArray{UInt64}(undef, n)

        CUDAnative.@sync begin
            @cuda threads=16 splitmix64_gpu(0x1234567890123456, n, result)
        end

        # Seed 0x1234567890123456:
        fromseed123 = [
            0xa747f481346acb72, 0x8ad918349ab73966, 0xbec62b7cc97a0873, 0xf4975b3f04f272de, 0x309a2a499503b323,
            0x75f0d35e440c1380, 0xa0e250c4832f2146, 0x5000ad7279d4564d, 0x527f776bf010fc54, 0xb178837968cf7887,
            0x03a2b6f219eb2ac9, 0x0c292f6dac0b89ce, 0xf406a38407ab8e94, 0x675980f999c91458, 0x888499c0bc4770c8,
            0x410c7af835212035, 0xf1a2cce43237861f, 0x1b3988311d0850c3, 0x70dc6bfcab337fc3, 0x12cec9bb9600463d,
            0xbce1012111683484, 0x63d63227c1131aef, 0x935018dead98b912, 0x50427855dd872b27, 0x5e1701eb8aed1c32,
            0x3231ecfe36e059f7, 0x3ee6eae41b75b100, 0x8640df45d78097ac, 0x6bd68158bdc48a49, 0xaa6af115f13bed47,
            0x9680ce74a4a7e996, 0x4b7d22f3cbd00a25, 0xc71682b11cc96c07, 0x663a228deabdf2b1, 0x75018ae2ecab3af5,
            0x0e6c0c485496b5db, 0x4f5f6e4b5283597f, 0x51a2c78fa67277c2, 0x7d6587ff16a97e90, 0x0226d182b0dd9a48,
            0x61b5565f8c28bd10, 0x9201bb21659aefd1, 0x3971d699e9f55c64, 0x0a3d12bae04ab182, 0xcbaf187d97a60d1e,
            0x2615a8f09a9a6b7a, 0xef5b6a1cbcabf696, 0x4dd54df163973e33, 0xe4b483b1d6762fbf, 0x3bd5b1a09a569677,
            0x49b90c90ce82c302, 0xf193d2253c266d22, 0xf885d93f2e56a330, 0xc6fdf5a84b29694a, 0x84b09a60cd8f6d50,
            0x859bdff3c39724b9, 0x8f7bf534dfdb6f8b, 0x162ade8198d9f7dd, 0x56d5af032eb87ad9, 0x961bf7f92b4da6cf,
            0x189850a17f576b2e, 0x52661a31fc19528f, 0x4486fe7a2234d24b, 0xea4d0bbe2643c140, 0xd682459a551543b1,
            0xee1c33c957b2b3c2, 0x5d3c39bcd8abf73b, 0xcfb078345aae27ae, 0xdd8cae3594569323, 0xb5c53b2f99372232,
            0x358f6866529407ca, 0xadf308a600d090a2, 0xe16cb21fb6a309a7, 0x07f86d1786ef4264, 0x317f3628ac7ccc19,
            0x8d78766af85d7ce4, 0x8d8790adea039d5e, 0x17c8ff07ebccb417, 0x93be8e429f5d8b3e, 0x8893278971571734,
            0xce4f4aab2da09bde, 0xe75a99e8d6afd288, 0xae74da6eea491702, 0xce202d270eca927c, 0x040cf01b89c75d41,
            0x2c9608b9e1138f11, 0x5a23a33b63d4cdb1, 0x5f1064ca5b023240, 0x180712765e0c7518, 0xb409002087be28db,
            0x15582962b3b3dd73, 0x1fc1c3a0880262c2, 0x99dd5b63f583202f, 0x99da0a71d47b3176, 0x7cd6a4aa5b2760c8,
            0xea5d86aa35596e3e, 0x28c3164807d7201a, 0x7203b1f76c6e2866, 0x3911c8344d36bd8b, 0x580cc832a21b96c6,
        ]

        result_host = Vector{UInt64}(result)

        for i = 1:n
            @test result_host[i] == fromseed123[i]
        end
    end

    @testset "Xoshiro256++ validation: First seed" begin
        n = 100
        result = CuArray{UInt64}(undef, n)

        CUDAnative.@sync begin
            @cuda threads=16 xoshiro256pp_gpu(0xe220a8397b1dcdaf, 0x6e789e6aa1b965f4, 0x06c45d188009454f, 0xf88bb8a8724c81ec, n, result)
        end

        # The seeds are taken from the SplitMix64 test vectors.
        # Seed 0xe220a8397b1dcdaf 0x6e789e6aa1b965f4 0x06c45d188009454f 0xf88bb8a8724c81ec:
        xshr256_test_1 = [
            0x53175d61490b23df, 0x61da6f3dc380d507, 0x5c0fdf91ec9a7bfc, 0x02eebf8c3bbe5e1a, 0x7eca04ebaf4a5eea,
            0x0543c37757f08d9a, 0xdb7490c75ab5026e, 0xd87343e6464bc959, 0x4b7da0a02389f0ff, 0x1300fc58c0424c16,
            0x5084843206c19968, 0x10ea073de9aa4dfc, 0x1aae554343960cc1, 0x1804139f10fae720, 0x10d790e7b8ac10fa,
            0x667d2bffdd1496f7, 0xa04620d3d0fc04a8, 0x1d50881230af9cc3, 0x53be287ded35f698, 0x673235793f7908e1,
            0x46e91feb4535fbdc, 0x216c1524cbac57c0, 0x0a53eb08063a44df, 0x45f965b948778197, 0x6f2fa9d01ba03887,
            0x60c57eba69ed4e15, 0x22c65ce977dd39cb, 0xa5d1ce0c5a7c6abf, 0xe8e26337cde13268, 0x0b4a575fdb6f8160,
            0x400feb0bae786424, 0x633e0b621080bf50, 0x5a456e5a144e059b, 0xdc75548b5cd2e8cd, 0xdf9d76f766648113,
            0x342bf8b7aec0de41, 0x831593e6b50ae928, 0x29e12b2a1872d7db, 0xb6362d8b640aec49, 0x2e78698eb5bba4a9,
            0x9064494b8287afb9, 0x4c04974c6c1b4767, 0x5863b8685408be73, 0x0e8ca571066bc302, 0x088959d638956a37,
            0x2e9392dfd5c30e86, 0x36da000d696e9d9e, 0x2a839b60548c1044, 0x3ebbaffcc5f270ca, 0x6da02738c0f92ee5,
            0x962fd83157fe1682, 0x856dcc088cece014, 0xca8717351ab24cbd, 0x231527552d018184, 0x06793b14839607ec,
            0xc54f89a7e193e5c1, 0xbacc209dd739707c, 0x7dc7053580f1ff20, 0x4ee696659cc1be91, 0xa3cb5d7769921646,
            0x9c002aaa8a687ded, 0xc0c3a216563d9ae2, 0x035b6d98ee8a1b19, 0x68d89ab6ea60f57d, 0xec45ebfcc75a43fc,
            0xc5baecd6bbbc6a26, 0x699d89305420574e, 0x48dce98c29e291a9, 0xb48e120ad1f9f23e, 0xe7c6368bd2aa145f,
            0xcf6e7679e474e800, 0x3d55f9180cf7bd5a, 0xd07064c4f0125475, 0x9784173bc20bb2cf, 0xc3296df2c4be97be,
            0x7ef38a3feaea3dcc, 0xa93ae5c3e3ac91b3, 0xc173a3156ae8e099, 0xa44373eaf0aff364, 0xb2a0bfb5000b2197,
            0xa49a7367daa7357f, 0x6be7b09a41e8061b, 0xa78c4d974c96b3fb, 0x28f7fbded613a963, 0x294d2574459f8cca,
            0xf3b457be6842a7ca, 0x9ce67b8b49364712, 0xd1b547099732fb15, 0x3538f30d9343b292, 0xc30c6e4c0dbee4c4,
            0x90b83117f578872c, 0xee166f8ef6076653, 0x671a16f58b431a4d, 0xa272fd92ddc2c04e, 0xf2bf4fd243cb8c2d,
            0xd2090a5c8c865eef, 0x7fd93d8b5401a273, 0x1ca67e32b0e47584, 0xb30ff462f8392fec, 0xb16ca64963ede86f,
        ]

        result_host = Vector{UInt64}(result)

        for i = 1:n
            @test result_host[i] == xshr256_test_1[i]
        end
    end

    @testset "Xoshiro256++ validation: Second seed" begin
        n = 100
        result = CuArray{UInt64}(undef, n)

        CUDAnative.@sync begin
            @cuda threads=16 xoshiro256pp_gpu(0xa747f481346acb72, 0x8ad918349ab73966, 0xbec62b7cc97a0873, 0xf4975b3f04f272de, n, result)
        end

        # Seed 0xa747f481346acb72 0x8ad918349ab73966 0xbec62b7cc97a0873 0xf4975b3f04f272de:
        xshr256_test_2 = [
            0x8764a3205cb8bb19, 0x83430004f370c12d, 0x4b25f419b4c33798, 0x8a513a4e0cb13f7d, 0x838ec104e7520e51,
            0x9d4e0fbd9efbd410, 0x4699cac599252763, 0x9d12a1920aa92ea6, 0x16f58890f28f73f0, 0xb9a3d53227a00f2a,
            0x0a13b7dbe3e75ea8, 0x79ffcf5f0c8d1100, 0x99ece4fbeb7313dc, 0xc3f87b278b8ee4e1, 0x16a11441a19431e6,
            0x67513ad840dac4c9, 0xc1e06281c8136e1c, 0x45fff01d45e28278, 0xbd34ecbec4c2dfd1, 0x8e1473a6621058bb,
            0x5a4d7b7f47cd940d, 0xc1895931810f51e4, 0xe2da937313d1c62b, 0x2457af5f33bd5dc8, 0xeee84ffaf344e299,
            0xba15cf4e32bdea22, 0xf6bdb5b21f9ec218, 0x1f57e7082a0e65e1, 0x0b0308f317319866, 0x9d7d1dd370af8c34,
            0x6dcaa675fa0b21f7, 0x81f6098b1b7fd554, 0x10b95cf55f11840f, 0x2618e9ca09439e81, 0xe3365a6aad20217a,
            0x91d9055d0163d663, 0x0e65d218a7775add, 0xd044ac77cd43fea5, 0x617118ea17e0b314, 0xc01d7ae0b8fa3785,
            0xacbe5ba699d219e2, 0xe7e11d4508e32276, 0xbfb68436cc59c458, 0x0625ba7690a02a54, 0x94a9adb4ad8289c3,
            0xfd9b4e6396f0b240, 0xb29486a69153d373, 0x541afedc53063ca5, 0x84b59a52889dddb7, 0x869d3adffcd07c97,
            0x13a2a3363624c07f, 0x9a5db06194d2cd66, 0x588fc64b9737ed95, 0xe1c08df219eb64bb, 0x361074055893ede1,
            0x9ca5480c4e77c5a8, 0x49f4b91279da85b3, 0xbb455871f21e0c38, 0x5ca49b23e994955b, 0xdb9261adea4fb9b6,
            0x1f439d23db8369a5, 0x4b23e1a8a2b567f0, 0xf7f9b8afc8a989f4, 0x2eb334a70476a1d5, 0x473826fbd4e3a603,
            0xf09a459b1860dd91, 0xc82e3ccb001f5a52, 0xd6cbdf410e30ce1d, 0x9b3bf1f29c0c524b, 0x896c2533de1e7e8a,
            0x9cdad573ce10cde1, 0xbb084c48360588db, 0x6140567a93f1c7a8, 0x66f960c75ec1a817, 0x3e6de19b9f056c71,
            0x3f563c8f3d00f8f2, 0xfa0db14f62a2aa2c, 0x4f085fd50a99a966, 0x9b4c8efdf3c99fb5, 0xbc18d391b6930fdc,
            0x75cfbab414b7cb69, 0x4c0894009af8eea7, 0x945283b2b9824dc7, 0x9db7b027499d18ca, 0x37510bb9f0fefc6f,
            0xfa2aceb93f756441, 0x604623f851ca0192, 0xebeeebf959db2c26, 0xb482490e0bdc79aa, 0x6ec59aa11bb8b941,
            0xe000f1d2445876af, 0x3810a2f6e9ca3e49, 0x54434cd018715317, 0xd243f33a759bc4ca, 0xb7859923fe236fb3,
            0x20d77e1aced47630, 0x31f0224fc82805e7, 0xcc6b181f051b66d0, 0x3c10226521a38a2c, 0x7a740af5acd8d654,
        ]

        result_host = Vector{UInt64}(result)

        for i = 1:n
            @test result_host[i] == xshr256_test_2[i]
        end
    end

    @testset "Generate PRNG from thread and block indices" begin
        threads = 128
        blocks = 8

        n = threads * blocks
        values_dev = CuArray{Float32}(undef, n)

        CuArrays.@sync begin
            @cuda threads=threads blocks=blocks prng_from_thread_and_block_gpu(values_dev, n)
        end

        values_host = Vector{Float32}(values_dev)
        values_set = Set{Float32}(values_host)

        # This ensures that each element in `values_dev` is different from all
        # other elements.
        @test length(values_set) == n
    end
end


text:
include("./sympatricpopulation.jl")

using .SympatricPopulation

plot_population_histories(
    1000, 2000, [.01 .05 .1 .2 .3 .4 .5 .6 .7 .8 .9 1. .21 .67], n_chromosomes=62÷2,
    n_turnover=10, n_iterations=100, sample_function=crosspollination_sampling,
    variance_plot=true, mean_heatmap=true,fixation_plot=true,
    hybridization_content_plot=true, postfix="primula")
plot_population_histories(
    1000, 2000, [.01 .05 .1 .2 .3 .4 .5 .6 .7 .8 .9 1. .30 .42], n_chromosomes=26÷2,
    n_turnover=10, n_iterations=100, sample_function=crosspollination_sampling,
    variance_plot=true, mean_heatmap=true,fixation_plot=true,
    hybridization_content_plot=true, postfix="rhododendron")
plot_population_histories(
    1000, 2000, [.01 .05 .1 .2 .3 .4 .5 .6 .7 .8 .9 1. .39 .85 .75], n_chromosomes=36÷2,
    n_turnover=10, n_iterations=100, sample_function=crosspollination_sampling,
    variance_plot=true, mean_heatmap=true,fixation_plot=true,
    hybridization_content_plot=true, postfix="gentiana")

ext: jl

max_stars_repo_name: biothomme/Geosmin

max_issues_repo_name: biothomme/Geosmin

max_forks_repo_name: biothomme/Geosmin


text:
using Literate

for t in readdir(".")
    if occursin(r".*_tut.jl", t)
        println("\nTutorial $t in $(pwd())\n")
        Literate.markdown(t, "."; documenter=false);
        # Literate.notebook(t, "."; execute=false, documenter=false);
        # @show n, ext = splitext(t)
        # @show "../docs/tutorials/" * n * ".md"
        cp(t, "../../../src/" * t, force = true)
    end
end


text:
"""
The second part of the Boltzmann Wealth Distribution example.

In the first part of the Boltzmann example, the experiment is ran without
a spatial structure. In the second part, a spatial structure is added,
and agents are required to only give money to agents who are on the
same node.

This example can be ran by navigating to the examples/ folder, starting
a julia REPL session and running:

```
julia> include("boltzmann_wealth_distribution_with_grid.jl")
```

This will instantiate the model and create a `DataFrame` `data` that
contains the result of running the model 10 steps. After running the
model, the results can be visualized in DataVoyager like this:

```
julia> visualize_data(data);
```

...which should result in a pop-up window that displays graphs
depicting the results of the experiment.

"""

using Agents

"""
Defines the agent type.

The agent type must be a subtype of AbstractAgent.

"""
mutable struct MyAgent{T<:Integer} <: AbstractAgent
  "The identifier number of the agent."
  id::T
  "The agent's grid position."
  pos::Tuple{T, T}
  "The agent's wealth."
  wealth::T
end

"Define the model type."
mutable struct MyModel{A<:AbstractVector, S<:AbstractSpace} <: AbstractModel
  "A space dimension."
  space::S
  "An array of agents."
  agents::A
  "A field for the scheduler function."
  scheduler::Function
end

"The space type that serves as the grid for the model."
mutable struct MyGrid{T<:Integer, Y<:AbstractVector} <: AbstractSpace
  "The dimensions of the grid."
  dimensions::Tuple{T, T}
  "A field for the space type."
  space::SimpleGraph
  "An array of arrays for each grid node."
  agent_positions::Y
end

"Function to instantiate the model."
function instantiate_model(; numagents, griddims)
  # An array of arrays for each node of the space.
  agent_positions = [Int64[] for i in 1:gridsize(griddims)]
  # Instantiate the grid structure.
  mygrid = MyGrid(griddims, grid(griddims), agent_positions)
  # Create a list of agents, each with position (1,1) and one unit of
  # wealth.
  agents = [MyAgent(i, (1,1), 1) for i in 1:numagents]

  # Instantiate and return the model.
  model = MyModel(mygrid, MyAgent[], random_activation)

  # Use the `add_agent!` function to add agents to the model.
  for agent in agents
    add_agent!(agent, model)
  end

  return model
end

"""
Define the agent step function.

Defines what the agent should do at each step.

"""
function agent_step!(agent::AbstractAgent, model::AbstractModel)
  # If the agent's wealth is zero, then do nothing.
  if agent.wealth == 0
    return
  # Otherwise..
  else
    #...create a list of all agents on the same node and select a random agent.
    available_agents = get_node_contents(agent, model)
    random_neighbor_agent_id = rand(available_agents)
    random_neighbor_agent = [i for i in model.agents
							 if i.id == random_neighbor_agent_id][1]
    # Then decrement the current agent's wealth and increment the neighbor's wealth.
    agent.wealth -= 1
	random_neighbor_agent.wealth += 1

    # Now move the agent to a random node.
    # If the agent weren't moved, agents would merely trade wealth
    # amongst themselves on the same node.
    neighboring_nodes = node_neighbors(agent, model)
    move_agent!(agent, rand(neighboring_nodes), model)
  end
end

# Instantiate the model.
model = instantiate_model(numagents=100, griddims=(5,5))

# For each agent, move the agent to a random location on the grid by using the
# `move_agent!` function.
for agent in model.agents
  move_agent!(agent, model)
end

# Run the model multiple steps and collect data.
# An array of Symbols for the agent fields that are to be collected, in
# this case wealth is the only variable to be collected.
agent_properties = [:wealth]
# Specifies at which steps data should be collected.
steps_to_collect_data = collect(1:10)
# Use the step function to run the model 10 times and collect data at
# each step.
data = step!(agent_step!, model, 10, agent_properties, steps_to_collect_data)


text:
module OpenStreetMapParser
    import LibExpat # XPCallbacks, XPStreamHandler, parsefile)
    import DataFrames # DataFrame, NA
    import LightGraphs

    export parseOSM, osm2dataframe, osmnetwork

    include("osm_classification.jl")
    include("osm_elements.jl")
    include("osm_parser.jl")
    include("osm_dataframes.jl")
    include("osm_rules.jl")
    include("osm_network.jl")
end # module


text:
include("partition.jl")
include("io.jl")

"""
    convergence(finame::String, foname::String; metric::String="snr", skipoverlap::Bool=false)

Compute the RMS error or SNR for progressively longer stacks of cross-correlations and (for RMS) a reference cross-correlation

# Arguments
- `finame::String,`    : Input file name for unstacked cross-correlations e.g. "./inputData/BPnetwork.jld2"
- `refname::String,`    : Input file name for reference cross-corrleations e.g. "./inputData/BPnetwork.jld2"
- `foname::String,`    : Output file name for fully stacked xcorrs e.g. "referenceXcorr.jld2"
- `metric::String`    : Method of tracking convergence: "snr" for peak signal-to-noise ratio, "rms" for root mean square misfit, or "cc" for correlation coefficient
- `skipoverlap::Bool`    : Whether or not to skip overlap in cross-correlation windows. If cross-correlations were computed
                           with no overlap, set this to false to ensure no time is skipped.
- `slice::Union{Bool, Float64, Array{Float64,1}}`    : whether to slice the cross-correlations before computing convergence. If so, a single Float64 specifies the
                                                       lag at which to slice to keep just the ballistic wave and an array of Float64 specifies the start and end lag
                                                       to keep only the coda.

# Output
- `foname.jld2`    : contains arrays of RMS errors for progressively longer stacks of cross-correlations

"""
function convergence(basefiname::String, refname::String, foname::String; ntimes::Int64=-1, metric::String="cc", skipoverlap::Bool=false, slice::Union{Bool, Float64, Array{Float64,1}}=false, thresh::Float64=-1)
    # input file contains metadata (stationlist, timestamplist)
    f = jldopen(basefiname*".jld2")
    tslist = f["info/timestamplist"] # time stamps
    close(f) # input file

    # read reference xcorrs
    f_ref = jldopen(refname)

    # dictionary to keep track of convergence vectors for each station pair
    conv = Dict()
    # dictionary to keep track of stacked data
    stackData = Dict{String, Array{Float64,2}}()

    if ntimes==-1
        ntimes=length(tslist)
    end
    nWins = Dict()
    nGood = Dict()
    # iterate over timestamps
    for tstamp in tslist[1:ntimes]
        println("Processing $tstamp")
        # read xcorrs for each time step
        f_cur = jldopen(basefiname*".$tstamp.jld2")
        stnpairs = sort(keys(f_cur[tstamp]))
        stnpairs = ["BP.CCRB..BP1.BP.SMNB..BP1"]

        # iterate over station pairs
        for stnpair in stnpairs
            # load unstacked cross-correlations
            data = try f_cur["$tstamp/$stnpair"] catch; continue end

            # load reference cross-correlation for the current station pair
            reference = f_ref[stnpair]

            # optionally slice reference to keep only coda or only ballistic
            ref = corr_slice(reference, 1:size(reference.corr, 2), slice)

            if haskey(nWins, stnpair)
                push!(nWins[stnpair], size(data.corr, 2))
            else
                nWins[stnpair] = [size(data.corr, 2)]
            end

            datacopy, cList = selective_stacking(data, reference, threshold=thresh)
            good_inds = findall(x->(x>thresh), cList)
            data.corr = data.corr[:, good_inds]

            if haskey(nGood, stnpair)
                push!(nGood[stnpair], length(good_inds))
            else
                nGood[stnpair] = [length(good_inds)]
            end

            # iterate over windowed cross-correlations
            for i=1:size(data.corr)[2]
                if skipoverlap i=2*i-1 end # assumes cc_step is 1/2 cc_len

                # slice current cross-correlation if necessary
                data_cur = corr_slice(data, i:i, slice)

                # build input to convergence functions
                # here we either create a progressively longer stack of cross-correlation
                # functions, or concatenate cross-correlation functions for snr
                if haskey(stackData, stnpair) && (metric=="rms" || metric=="cc")
                    # keep running sum of cross-correlations, stack if key exists, assign if not
                    stackData[stnpair] .+= data_cur
                elseif haskey(stackData, stnpair) && metric=="snr"
                    # create matrix of xcorrs
                    stackData[stnpair] = hcat(stackData[stnpair], data_cur)
                else
                    # reshape into 2D array (necessary for snr input)
                    if typeof(slice) == Array{Float64,1}
                        stackData[stnpair] = data_cur
                    else
                        stackData[stnpair] = reshape(data_cur, length(data_cur), 1)
                    end
                end

                # compute convergence data point and append to vector
                if metric=="rms"
                    # compute RMS error of the running sum with the reference
                    if haskey(conv, stnpair)
                        push!(conv[stnpair], rms(normalize(ref[:,1]), normalize(stackData[stnpair][:,1])))
                    else
                        conv[stnpair] = [rms(normalize(ref[:,1]), normalize(stackData[stnpair][:,1]))]
                    end
                elseif metric=="cc"
                    # compute Pearson correlation-coefficient between current and reference cross-correlation
                    if haskey(conv, stnpair)
                        if typeof(slice)==Array{Float64,1}
                            # compute correlation-coefficient for both positive and negative coda
                            # store in 2d array
                            neg_coda_cor = cor(ref[:, 1], stackData[stnpair][:, 1])
                            pos_coda_cor = cor(ref[:, 2], stackData[stnpair][:, 2])
                            push!(conv[stnpair], [neg_coda_cor, pos_coda_cor])
                        else
                            push!(conv[stnpair], cor(ref[:,1], stackData[stnpair][:,1]))
                        end
                    else
                        if typeof(slice)==Array{Float64,1}
                            # compute correlation-coefficient for both positive and negative coda
                            # store in 2d array
                            neg_coda_cor = cor(ref[:, 1], stackData[stnpair][:, 1])
                            pos_coda_cor = cor(ref[:, 2], stackData[stnpair][:, 2])
                            conv[stnpair] = [[neg_coda_cor, pos_coda_cor]]
                        else
                            conv[stnpair] = [cor(ref[:,1], stackData[stnpair][:,1])]
                        end
                    end
                elseif metric=="snr"
                    # compute peak signal-to-noise ratio of increasingly long stacks
                    if haskey(conv, stnpair)
                        # compute psnr for both positive and negative coda
                        if typeof(slice)==Array{Float64,1}
                            neg_coda = findall(x->(x%2!=0), collect(1:size(stackData[stnpair], 2)))
                            pos_coda = findall(x->(x%2==0), collect(1:size(stackData[stnpair], 2)))
                            neg_coda_snr = maximum(snr(stackData[stnpair][:, neg_coda], reference.fs))
                            pos_coda_snr = maximum(snr(stackData[stnpair][:, pos_coda], reference.fs))

                            push!(conv[stnpair], [neg_coda_snr, pos_coda_snr])
                        else
                            push!(conv[stnpair], maximum(snr(stackData[stnpair], reference.fs)))
                        end
                    else
                        if typeof(slice)==Array{Float64,1}
                            neg_coda = findall(x->(x%2!=0), collect(1:size(stackData[stnpair], 2)))
                            pos_coda = findall(x->(x%2==0), collect(1:size(stackData[stnpair], 2)))
                            neg_coda_snr = maximum(snr(stackData[stnpair][:, neg_coda], reference.fs))
                            pos_coda_snr = maximum(snr(stackData[stnpair][:, pos_coda], reference.fs))

                            conv[stnpair] = [[neg_coda_snr, pos_coda_snr]]
                        else
                            conv[stnpair] = [maximum(snr(stackData[stnpair], reference.fs))]
                        end
                    end
                end
            end
        end
        close(f_cur) # xcorrs for each time step
    end
    close(f_ref) # reference xcorrs

    # write convergence vectors to disk
    f_conv = jldopen(foname, "a+")
    for key in keys(conv)
        f_conv["$metric/$key"] = conv[key]
        f_conv["$metric/nWin/$key"] = nWins[key]
        f_conv["$metric/nGood/$key"] = nGood[key]
    end

    close(f_conv) # convergence file
end

function corr_slice(data::CorrData, range::UnitRange{Int64}, slice::Union{Bool, Float64, Array{Float64,1}}=false)
    if typeof(slice)==Float64
        # time vector
        tvec = -data.maxlag:1/data.fs:data.maxlag
        # find all times within [-slice,slice] time
        t_inds = findall(x->(x>=-slice && x<=slice), tvec)
        # slice reference (maintain 2d array)
        d = data.corr[t_inds, range]
    elseif typeof(slice)==Array{Float64,1}
        # convert startlag/endlag[s] to startlag/windowlength[samples]
        win_len = Int(diff(slice)[1] * data.fs)
        startlag = Int(slice[1] * data.fs)
        # partition reference to keep only the coda
        d = partition(data.corr[:, range], startlag, win_len)
    else
        println("Proceeding without slicing.")
        d = data.corr
    end
    return d
end

max_stars_repo_name: kura-okubo/SeisXcorrelation

max_issues_repo_name: kura-okubo/SeisXcorrelation

max_forks_repo_name: kura-okubo/SeisXcorrelation


text:
"""
```
augment_variables!(m::BruSan, stategrid::StateGrid, funcvar::OrderedDict{Symbol, Vector{S}},
    endo::OrderedDict{Symbol, Vector{S}}) where {S <: Real}

augment_variables_nojump!(m::BruSan, stategrid::StateGrid, ode_f::Function, funcvar::OrderedDict{Symbol, Vector{S}},
    endo::OrderedDict{Symbol, Vector{S}}, odesol::ODESolution) where {S <: Real}
```

calculates derivatives and additional endogenous variables that can be computed
after solving the equilibrium system of functional equations.
"""
function augment_variables!(m::BruSan, stategrid::StateGrid, funcvar::OrderedDict{Symbol, Vector{S}},
                            derivs::OrderedDict{Symbol, Vector{S}},
                            endo::OrderedDict{Symbol, Vector{S}}) where {S <: Real}
end

"""
```
augment_variables_nojump!(m::BruSan, stategrid::StateGrid, ode_f::Function, funcvar::OrderedDict{Symbol, Vector{S}},
    endo::OrderedDict{Symbol, Vector{S}}, odesol::ODESolution) where {S <: Real}
```

calculates derivatives and additional endogenous variables that can be computed
after solving the equilibrium system of functional equations. The first derivative
of `p` is directly calculated using `ode_f`
"""
function augment_variables_nojump!(m::BruSan, stategrid::StateGrid, ode_f::Function, funcvar::OrderedDict{Symbol, Vector{S}},
                                   derivs::OrderedDict{Symbol, Vector{S}},
                                   endo::OrderedDict{Symbol, Vector{S}}, odesol::ODESolution) where {S <: Real}

    if m[:γₑ].value == 1. && m[:γₕ].value == 1. && m[:ψₑ].value == 1. && m[:ψₕ].value == 1.
        return augment_variables_nojump_log!(m, stategrid, ode_f, funcvar, derivs, endo, odesol)
    elseif m[:ψₑ].value == 1. && m[:ψₕ].value == 1.
        return augment_variables_nojump_nonunit_gamma!(m, stategrid, ode_f, funcvar, derivs, endo, odesol)
    end
end

function augment_variables_nojump_log!(m::BruSan, stategrid::StateGrid, ode_f::Function, funcvar::OrderedDict{Symbol, Vector{S}},
                                       derivs::OrderedDict{Symbol, Vector{S}},
                                       endo::OrderedDict{Symbol, Vector{S}}, odesol::ODESolution) where {S <: Real}

    # Unpack equilibrium endogenous variables and stategrid
    q       = funcvar[:q]
    ∂q_∂η   = derivs[:∂q_∂η]
    ∂²q_∂η² = derivs[:∂²q_∂η²]
    η       = stategrid[:η]
    @unpack φ_e, φ_h, μ_q, σ_q, ςₑ, ςₕ, μ_η, σ_η, ι, Φ, μ_K = endo
    θ = parameters_to_named_tuple(m.parameters)
    @unpack ρₑ, ρₕ, νₑ, νₕ, aₑ, aₕ, σ, δ, χ₁, χ₂, τ = θ
    cₑ_per_nₑ = ρₑ + νₑ
    cₕ_per_nₕ = ρₕ + νₕ

    # Interpolate remainder of the solution for q from odesol
    i = findfirst(η .>= odesol.t[end])
    q[1:i] .= odesol(η[1:i])
    ∂q_∂η[1:i] .= map(j -> ode_f(q[j], θ, η[j]), 1:i)
    if cₑ_per_nₑ == cₕ_per_nₕ
        q[(i + 1):end] .= boundary_conditions(m)[:q][2]
        ∂q_∂η[(i + 1):end] .= 0.
    else
        q[(i + 1):end] .= (χ₂ * aₑ + 1.) ./ (χ₂ .* (η[(i + 1):end] .* cₑ_per_nₑ +
                                                    (1. .- η[(i + 1):end]) .* cₕ_per_nₕ) .+ χ₁)
        ∂q_∂η[(i + 1):end] .= -(χ₂ * aₑ + 1) / (χ₂ * (η  * cₑ_per_nₑ + (1 - η) * cₕ_per_nₕ) + χ₁)^2 *
            (χ₂ * (cₑ_per_nₑ - cₕ_per_nₕ))
    end
    dx = vcat(η[1], diff(η), 1. - η[end])
    ∂²q_∂η² .= CenteredDifference(2, 2, dx, length(η)) * RobinBC((0., 1., 0.), (0., 1., 0.), dx) * q

    # Solve for values not calculated during the calculation of equilibrium
    φ_e[1:i] .= ((q[1:i] .* (χ₂ .* (cₑ_per_nₑ .* η[1:i] + cₕ_per_nₕ .* (1. .- η[1:i])) .+ χ₁) .- 1.) ./ χ₂ .- aₕ) ./
        (aₑ - aₕ) ./ η[1:i]
    φ_e[(i + 1):end] .= 1. ./ η[(i + 1):end]
    φ_h .= (1. .- φ_e .* η) ./ (1. .- η)
    is1 = (i + 1):length(η) # experts hold all capital
    no1 = 1:i
    σ_q[no1] .= sqrt.((aₑ - aₕ) ./ q[no1] ./ (φ_e[no1] - φ_h[no1])) .- σ
    σ_q[is1] .= ∂q_∂η[is1] .* η[is1] .* (φ_e[is1] .- 1.) .* σ ./
        (q[is1] .- ∂q_∂η[is1] .* η[is1] .* (φ_e[is1] .- 1.))
    ςₑ .= φ_e .* (σ .+ σ_q).^2
    ςₕ .= φ_h .* (σ .+ σ_q).^2

    ι   .= (χ₁ .* q .- 1.) ./ χ₂
    Φ   .= (χ₁ / χ₂) .* log.(χ₁ .* q)
    μ_K .= Φ .- δ

    σ_η .= (φ_e .- 1.) .* (σ .+ σ_q)
    μ_η .= (aₑ .- ι) ./ q .- cₑ_per_nₑ .- τ .+ σ_η.^2
    μ_q .= ∂q_∂η .* μ_η .* η + ∂²q_∂η² ./ 2 .* (σ_η .* η).^2
end


text:
import Base.split

export split

"""
   c = split(s::SeisCon, inds)

Creates a new SeisCon object `c` without copying by referencing `s.blocks[inds]`

`inds` can be an array of indicies, a range, or a scalar.

# Example
```julia-repl
julia> b = split(s, 1:10);
julia> c = split(b, [1; 7; 9]);
julia> d = split(c, 1);
```
And we can confirm that `d` and `s` reference the same place in memory.

```julia-repl
julia> s.blocks[1] === d.blocks[1]
true
```
"""
function split(s::SeisCon, inds::Union{Vector{Ti}, AbstractRange{Ti}}) where {Ti<:Integer}
    c = SeisCon(s.ns, s.dsf, view(s.blocks, inds), s.fileheader)
end

split(s::SeisCon, inds::Integer) = split(s, [inds])

"""
    split(s::SeisBlock, inds::Union{Vector{Ti}, AbstractRange{Ti}};
               consume::Bool = false) where {Ti<:Integer}

Split the 'ind' traces of `s` into a sepretate `SeisBlock` object.

The default behaviour does not change the input object `s`. The returned `SeisBlock` is
created by referencing the traceheaders of `s`, and copying the subset of the data.

If memory use is a concern, the `consume` keyword changes the behaviour of `merge`.
If `consume` is true, then the `ind` traces of `s` are REMOVED and placed into the returned
`SeisBlock` object. In this case, `merge` will operate in place on `s`, while still returning a
 subset of `s`.

# Examples

    julia> using SegyIO

    julia> s = segy_scan(joinpath(dirname(pathof(SegyIO)),"data/"), "overthrust", ["GroupX"; "GroupY"], verbosity = 0);

    julia> d = s[1:length(s)]; @time b = split(d, 1:10000);
      0.005683 seconds (23 allocations: 28.649 MiB, 17.99% gc time)

    julia> println("d: ", size(d)); println("b: ", size(b))
    d: (751, 20013)
    b: (751, 10000)

    julia> data_in_memory = Base.summarysize(d.data) + Base.summarysize(b.data)
    90159052

Note that the default behaviour of `merge` is to create `b` by copying `d`.

    julia> d = s[1:length(s)]; @time b = split(d, 1:10000, consume=true);
      0.062364 seconds (37 allocations: 116.537 MiB, 19.44% gc time)

    julia> println("d: ", size(d)); println("b: ", size(b))
    d: (751, 10013)
    b: (751, 10000)

    julia> data_in_memory = Base.summarysize(d.data) + Base.summarysize(b.data)
    60119052

In this case, because `consume` was set to true, `merge` created `b` by removing traces
from `d`. `consume` prevents the duplication of data in memory at the cost of performance,
memory allocation, and risk.
"""
function split(s::SeisBlock, inds::Union{Vector{Ti}, AbstractRange{Ti}};
               consume::Bool = false) where {Ti<:Integer}
    if consume
        c = SeisBlock(s.fileheader, s.traceheaders[inds], s.data[:, inds])
        ii = BitArray(size(s.data))
        ii[:, inds] = true
        deleteat!(vec(s.data), vec(ii))
        deleteat!(s.traceheaders, inds)
        s.data = s.data[:, 1:end-length(inds)]
    else
        c = SeisBlock(s.fileheader, view(s.traceheaders, inds), s.data[:, inds])
    end

    return c
end

split(s::SeisBlock, inds::Integer) = split(s, [inds])


text:
#==============================================================================#
# CUR.jl
#
# This file is generated from:
# https://github.com/aws/aws-sdk-js/blob/master/apis/cur-2017-01-06.normal.json
#==============================================================================#

__precompile__()

module CUR

using AWSCore

"""
    using AWSSDK.CUR.delete_report_definition
    delete_report_definition([::AWSConfig], arguments::Dict)
    delete_report_definition([::AWSConfig]; <keyword arguments>)

    using AWSCore.Services.cur
    cur([::AWSConfig], "DeleteReportDefinition", arguments::Dict)
    cur([::AWSConfig], "DeleteReportDefinition", <keyword arguments>)

# DeleteReportDefinition Operation

Delete a specified report definition

# Arguments

## `ReportName = ::String`

# Returns

`DeleteReportDefinitionResponse`

# Exceptions

`InternalErrorException` or `ValidationException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/cur-2017-01-06/DeleteReportDefinition)
"""
@inline delete_report_definition(aws::AWSConfig=default_aws_config(); args...) = delete_report_definition(aws, args)

@inline delete_report_definition(aws::AWSConfig, args) = AWSCore.Services.cur(aws, "DeleteReportDefinition", args)

@inline delete_report_definition(args) = delete_report_definition(default_aws_config(), args)

"""
    using AWSSDK.CUR.describe_report_definitions
    describe_report_definitions([::AWSConfig], arguments::Dict)
    describe_report_definitions([::AWSConfig]; <keyword arguments>)

    using AWSCore.Services.cur
    cur([::AWSConfig], "DescribeReportDefinitions", arguments::Dict)
    cur([::AWSConfig], "DescribeReportDefinitions", <keyword arguments>)

# DescribeReportDefinitions Operation

Describe a list of report definitions owned by the account

# Arguments

## `MaxResults = ::Int`

## `NextToken = ::String`

# Returns

`DescribeReportDefinitionsResponse`

# Exceptions

`InternalErrorException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/cur-2017-01-06/DescribeReportDefinitions)
"""
@inline describe_report_definitions(aws::AWSConfig=default_aws_config(); args...) = describe_report_definitions(aws, args)

@inline describe_report_definitions(aws::AWSConfig, args) = AWSCore.Services.cur(aws, "DescribeReportDefinitions", args)

@inline describe_report_definitions(args) = describe_report_definitions(default_aws_config(), args)

"""
    using AWSSDK.CUR.put_report_definition
    put_report_definition([::AWSConfig], arguments::Dict)
    put_report_definition([::AWSConfig]; ReportDefinition=)

    using AWSCore.Services.cur
    cur([::AWSConfig], "PutReportDefinition", arguments::Dict)
    cur([::AWSConfig], "PutReportDefinition", ReportDefinition=)

# PutReportDefinition Operation

Create a new report definition

# Arguments

## `ReportDefinition = [ ... ]` -- *Required*

```
 ReportDefinition = [
        "ReportName" => <required> ::String,
        "TimeUnit" => <required> "HOURLY" or "DAILY",
        "Format" => <required> "textORcsv",
        "Compression" => <required> "ZIP" or "GZIP",
        "AdditionalSchemaElements" => <required> ["RESOURCES", ...],
        "S3Bucket" => <required> ::String,
        "S3Prefix" => <required> ::String,
        "S3Region" => <required> "us-east-1", "us-west-1", "us-west-2", "eu-central-1", "eu-west-1", "ap-southeast-1", "ap-southeast-2" or "ap-northeast-1",
        "AdditionalArtifacts" =>  ["REDSHIFT" or "QUICKSIGHT", ...]
    ]
```

# Returns

`PutReportDefinitionResponse`

# Exceptions

`DuplicateReportNameException`, `ReportLimitReachedException`, `InternalErrorException` or `ValidationException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/cur-2017-01-06/PutReportDefinition)
"""
@inline put_report_definition(aws::AWSConfig=default_aws_config(); args...) = put_report_definition(aws, args)

@inline put_report_definition(aws::AWSConfig, args) = AWSCore.Services.cur(aws, "PutReportDefinition", args)

@inline put_report_definition(args) = put_report_definition(default_aws_config(), args)

end # module CUR

#==============================================================================#
# End of file
#==============================================================================#


text:
# function record(data)
#   push!(data)
# end

# struct MemoryChannel{T}
#   x::Vector{T}
# end

# function (x::MemoryChannel{T})(c::Channel) where T
#   while true
#     v = rand(T)
#     push!(x.x, v)
#     put!(c, v)
#   end
# end

# mem = MemoryChannel([31])

max_stars_repo_name: zenna/expect

max_issues_repo_name: zenna/expect

max_forks_repo_name: zenna/expect


text:
# Second iteration of explicitly exloring the marginalization of discrete variables onto the continuous state space.
# Although this code is still excessive and messy, this is a significant feature expansion; and future versions will
# generalize the marginalization process to allow for implicit hypothesis exploration.  The messy explicit version of code is
# intended to help develop the required generalistic unit tests.  The unit tests will then be validated, frozen and used to
# confirm that future "algebraic" marginalization (implicit) versions operate correctly.
# FYI, the complexity of general multihypothesis convolutions can be deceiving, however, note that the coding
# complexity is contained for each indivual factor at a time.  Global Bayes tree inference then creates the symphony
# of non-Gaussian (multimodal) posterior beliefs from the entire factor graph.
#
# 2018/6/01 @dehann

"""
    $SIGNATURES

Return common vectors `(allmhp, certainidx,uncertnidx)` used for dealing with multihypo cases.
"""
function getHypothesesVectors(mhp::Vector{Float64})::Tuple{Vector{Int},Vector{Int},Vector{Int}}
  allmhp = 1:length(mhp)
  certainidx = allmhp[mhp .== 0.0]  # TODO remove after gwp removed
  uncertnidx = allmhp[0.0 .< mhp]
  return (allmhp,certainidx,uncertnidx)
end

"""
    $(SIGNATURES)

This function explicitly encodes the marginalization of a discrete categorical selection variable
for ambiguous data association situations.  This function populates `allelements` with particle
indices associated with particular multihypothesis selection while `activehypo` simultaneously
contains the hypothesis index and factor graph variables associated with that hypothesis selection.  The return value `certainidx` are the hypotheses that are not in question.

This function does not consider whether a unroll hypo lambda can actually exist, it just generates
a recipe of which options should be considered.  Whoever solves the recipe is responsible for
building the right lambda or doing something else.

Input:
- `maxlen` is the max number of samples across all variables

Output:
- `certainidx`:   non fractional variables
- `allelements`:  list of which particles go with which hypothesis selection
- `activehypo`:   list of which hypotheses go with which certain + fractional variables
- `mhidx`:        multihypothesis selection per particle idx
- `sfidx`:        Solve for idx

Example:
```julia
idx=(1,2,3)
multihypo=[1.0;0.5;0.5]    # X,La,Lb

# specfic example
# this is important -- e.g. `pts = approxConv(fg, :XLaLbf1, :La)`
sfidx=2

certainidx=
1-element Array{Int64,1}:
 1                         # X
allelements=
3-element Array{Any,1}:
 Int64[]
 [1, 2, 11, ...]
 [3, 4, 5, ...]
activehypo=
3-element Array{Any,1}:
 (0, [2])     # nullhypo -- forced afterward, might be deprecated if better solution is found
 (1, [1, 2])  # unroll hypo lambda for X,La
 (2, [1, 2])  # unroll hypo lambda for X,La
 (3, [2, 3])  # would be (but cannot) unroll hypo La,Lb # almost our nullhypo # wont build a lambda

# now select which lambdas to build based on converted `activehypo`  rand(  Categorical( [0; 0.5;0.5] ) )
 mhidx=
100-element Array{Int64,1}:
2, 2, 3, 3, 3,...
```

Another example on what is `certainidx`
```julia
multihypo=[1;1;0.5;0.5]
# results in: `certainidx = [1;2]`
```

Notes:
- Issue 427, race condition during initialization since n-ary variables not resolvable without other init.

DevNotes
- FIXME convert into some kind of `HypoRecipe` struct, to improve readibility and code maintainability
- TODO add nullhypo cases to returning result
- FIXME make type-stable `activehypo` and others
- Improved implementations should implicitly induce the same behaviour through summation (integration) when marginalizing any number of discrete variables.

```
# `allelements` example BearingRange [:x1, 0.5:l1a, 0.5:l1b]
# sfidx = (1=:x1,2=:l1a,3=:l1b)
if solvefor :x1, then allelem = [mhidx.==:l1a; mhidx.==l1b]
if solvefor :l1a, then allelem = [mhidx.==:l1a] and ARR[solvefor][:,mhidx.==:l1b]=ARR[:l1b][:,mhidx.==:l1b]
if solvefor :l1b, then allelem = [mhidx.==:l1b] and ARR[solvefor][:,mhidx.==:l1a]=ARR[:l1a][:,mhidx.==:l1a]
if solvefor 1, then allelem = [mhidx.==2; mhidx.==3]
if solvefor 2, then allelem = [mhidx.==2] and ARR[solvefor][:,mhidx.==3]=ARR[3][:,mhidx.==3]
if solvefor 3, then allelem = [mhidx.==3] and ARR[solvefor][:,mhidx.==2]=ARR[2][:,mhidx.==2]

# `activehypo` in example mh=[1.0;0.5;0.5]
sfidx=1, mhidx=2:  ah = [1;2]
sfidx=1, mhidx=3:  ah = [1;3]
sfidx=2, mhidx=2:  ah = [1;2]
sfidx=2, mhidx=3:  2 should take a value from 3 or nullhypo
sfidx=3, mhidx=2:  3 should take a value from 2 or nullhypo
sfidx=3, mhidx=3:  ah = [1;3]

# `activehypo` in example mh=[1.0;0.33;0.33;0.34]
sfidx=1, mhidx=2:  ah = [1;2]
sfidx=1, mhidx=3:  ah = [1;3]
sfidx=1, mhidx=4:  ah = [1;4]

sfidx=2, mhidx=2:  ah = [1;2]
sfidx=2, mhidx=3:  2 should take a value from 3 or nullhypo
sfidx=2, mhidx=4:  2 should take a value from 4 or nullhypo

sfidx=3, mhidx=2:  3 should take a value from 2 or nullhypo
sfidx=3, mhidx=3:  ah = [1;3]
sfidx=3, mhidx=4:  3 should take a value from 4 or nullhypo

sfidx=4, mhidx=2:  4 should take a value from 2 or nullhypo
sfidx=4, mhidx=3:  4 should take a value from 3 or nullhypo
sfidx=4, mhidx=4:  ah = [1;4]
```

Also keeping the default case documented:
```
# the default case where mh==nothing
# equivalent to mh=[1;1;1] # assuming 3 variables
sfidx=1, allelements=allidx[nhidx.==0], activehypo=(0,[1;])
sfidx=2, allelements=allidx[nhidx.==0], activehypo=(0,[2;])
sfidx=3, allelements=allidx[nhidx.==0], activehypo=(0,[3;])
```

TODO still need to compensate multihypo case for user nullhypo addition.
"""
function _prepareHypoRecipe!( mh::Categorical,
                                      maxlen::Int,
                                      sfidx::Int,
                                      lenXi::Int,
                                      isinit::Vector{Bool}=ones(Bool, lenXi),
                                      nullhypo::Real=0  )
  #
  allelements = []
  activehypo = []
  mhidx = Int[]

  allidx = 1:maxlen
  allmhp, certainidx, uncertnidx = getHypothesesVectors(mh.p)

  # select only hypotheses that can be used (ie variables have been initialized)
  @assert !(sum(isinit) == 0 && sfidx == certainidx) # cannot init from nothing for any hypothesis

  mhh = if sum(isinit) < lenXi - 1
    @assert isLeastOneHypoAvailable(sfidx, certainidx, uncertnidx, isinit)
    @info "not all hypotheses initialized, but at least one available -- see #427"
    mhp = deepcopy(mh.p)
    suppressmask = isinit .== false
    suppressmask[sfidx] = false
    mhp[suppressmask] .= 0.0
    mhp ./= sum(mhp)
    Categorical(mhp)
  else
    mh
  end

  # FIXME consolidate with addEntropyOnManifolds approach in `computeAcrossHypothesis!`
  # prepend for the mhidx=0, bad-init-null-hypothesis case (if solving a fractional variable)
  mhh = if sfidx in uncertnidx
    nhw = (length(uncertnidx)+1)
    nmhw = [1/nhw; length(uncertnidx)/nhw*mhh.p]
    nmhw ./= sum(nmhw) # renormalize (should not be necessary)
    Categorical(nmhw)
  else
    mhh
  end

  # prep mm-nultihypothesis selection values
  mhidx = rand(mhh, maxlen)  # selection of which hypothesis is correct
  pidx = 0
  if sfidx in uncertnidx
    # shift down to get mhidx=0 case
    mhidx .-= 1
    pidx = -1
  end

  sfincer = sfidx in certainidx
  for pval in mhh.p
    pidx += 1
    pidxincer = pidx in certainidx # ??
    # permutation vectors for later computation
    iterarr = allidx[mhidx .== pidx]
    iterah = Int[]
    if !pidxincer && sfincer && pidx != 0 # 1e-15 <= pval && mh.p[sfidx] < 1e-10  # proxy for sfidx in certainidx
      # solve for one of the certain variables containing uncertain hypotheses in others
      iterah = sort(union(certainidx, pidx)) # sort([sfidx;pidx])
      # DONE -- supports n-ary factors in multihypo mode
    elseif (pidxincer && !sfincer || sfidx == pidx) && pidx != 0 # pval < 1e-15 && mh.p[sfidx] >= 1e-10
      # solve for one of the uncertain variables
      iterah = sort(union(certainidx, sfidx)) # sort([sfidx;pidx])
      # EXPERIMENTAL -- support more than binary factors in multihypo mode
    elseif pidxincer && sfincer && pidx != 0 # pval < 1e-15 && mh.p[sfidx] < 1e-10
      iterarr = Int[]
      iterah = Int[] # may be moot anyway, but double check first
    elseif !pidxincer && !sfincer && pidx != 0 # pval >= 1e-15 && mh.p[sfidx] >= 1e-10
      iterah = uncertnidx #allmhp[mh.p .> 1e-15]
    elseif pidx == 0
      # nullhypo only take values from self sfidx, might add entropy later (for bad init case)
      iterah = [sfidx;]
    else
      error("Unknown hypothesis case, got sfidx=$(sfidx) with mh.p=$(mh.p), pidx=$(pidx)")
    end
    push!(allelements, iterarr)
    push!(activehypo, (pidx,iterah))
  end

  # # retroactively add nullhypo case (the 0 case)
  # if sfidx in uncertnidx
  #   #
  # end

  # hyporecipe::NamedTuple
  return (;certainidx, allelements, activehypo, mhidx)
end

function _prepareHypoRecipe!( mh::Nothing,
                                      maxlen::Int,
                                      sfidx::Int,
                                      lenXi::Int,
                                      isinit::Vector{Bool}=ones(Bool, lenXi),
                                      nullhypo::Real=0  )
  #
  # FIXME, consolidate with the general multihypo case

  # the default case where mh==nothing
  # equivalent to mh=[1;1;1] # assuming 3 variables
  # sfidx=1, allelements=allidx[nhidx.==0], activehypo=(0,[1;])

  #
  allelements = []
  activehypo = []

  # TODO add cases where nullhypo occurs, see DFG #536, and IIF #237
  nmhw = [nullhypo; (1-nullhypo)]
  nhh = Categorical(nmhw)

  # prep mmultihypothesis selection values
  # mhidx = Int[]
  # NOTE, must do something special to get around Categorical([0;10]) error
    # selection of which hypothesis is correct
  mhidx = nullhypo == 0 ? ones(Int, maxlen) : (rand(nhh, maxlen) .- 1)

  allidx = 1:maxlen
  certainidx = 1:lenXi
  # zero is nullhypo case, 1 is first sfidx variable
  nullarr = allidx[mhidx .== 0]
  # mhidx == 1 case is regular -- this will be all elements if nullhypo=0.0
  reguarr = allidx[mhidx .!= 0]
  pidxAll = [0;certainidx]
  for pidx in pidxAll
    if pidx == 0
      # elements that occur during nullhypo active
      push!(allelements, nullarr)
      push!(activehypo, (pidx,[sfidx;]))
    elseif pidx == 1
      # elements that occur during regular hypothesis true
      push!(allelements, reguarr)
      push!(activehypo, (pidx,certainidx))
    else
      # all remaining collections are empty (part of multihypo support)
      push!(allelements, Int[])
      push!(activehypo, (pidx,Int[]))
    end
  end

  # return hyporecipe::NamedTuple
  return (;certainidx, allelements, activehypo, mhidx)
end

#


text:
println("Installing all my favorite libraries...")
using Pkg
Pkg.add("Agents")
Pkg.add("BandedMatrices")
Pkg.add("BenchmarkTools")
Pkg.add("Cbc")
Pkg.add("CDCS")
Pkg.add("Conda")
Pkg.add("CSV")
Pkg.add("DataFrames")
Pkg.add("DataFramesMeta")
Pkg.add("DataStructures")
Pkg.add("Distributions")
Pkg.add("Flux")
Pkg.add("FredData")
Pkg.add("FreqTables")
Pkg.add("Genie")
Pkg.add("GLM")
Pkg.add("GLPK")
Pkg.add("GR")
Pkg.add("GraphPlot")
Pkg.add("HTTP")
Pkg.add("IJulia")
Pkg.add("Ipopt")
Pkg.add("InteractiveUtils")
Pkg.add("JuMP")
Pkg.add("Juno")
Pkg.add("LightGraphs")
Pkg.add("Logging")
Pkg.add("OpenStreetMapX")
Pkg.add("Parameters")
Pkg.add("PGFPlots")
Pkg.add("Plots")
Pkg.add("Pluto")
Pkg.add("PlutoUI")
Pkg.add("PyCall")
Pkg.add("PyPlot")
Pkg.add("RCall")
Pkg.add("ScikitLearn")
Pkg.add("SimpleHypergraphs")
Pkg.add("SparseArrays")
Pkg.add("StatsBase")
Pkg.add("StatsKit")
Pkg.add("Turing")
Pkg.add("XLSX")

using Conda
println("Installing some libraries with conda...")
Conda.runconda(`install jupyter --yes`)
Conda.runconda(`install matplotlib --yes`)
Conda.runconda(`install folium -c conda-forge --yes`)

println("Precompiling all packages...")
Pkg.precompile()


text:
using Base: Int64
using XLSX: length
#using Pkg
#Pkg.add("CSV")
#Pkg.add("DataFrames")
#Pkg.add("XLSX")
using XLSX, DataFrames, Tables, SymPy;
IO = XLSX.readdata("C:\\Users\\jaber\\OneDrive\\Documents\\AIBE\\io-table\\IO.xlsx", "io-table-5!A1:DV130");

#indexing vectors for initial data import groups
intermediaryCol = findall(x -> occursin("0", x), string.(IO[1,:]));
intermediaryRow = findall(x -> occursin("0", x), string.(IO[:,1]));
finalDemandCol = findall(x -> occursin('Q', x), string.(IO[3,:]));
factorRow = findall(x -> occursin('P', x), string.(IO[:,1]));
#initialising SAM
SAMDims = length(intermediaryCol) + length(finalDemandCol) + length(factorRow) + 1;
SAM = zeros(SAMDims,SAMDims);
#import numerical data into SAM
SAM[1:length(intermediaryRow), 1:length(intermediaryCol)] = IO[intermediaryRow, intermediaryCol];
SAM[1:length(intermediaryRow), length(intermediaryCol)+length(factorRow)+1:length(intermediaryCol)+length(factorRow)+ length(finalDemandCol)] = IO[intermediaryRow, finalDemandCol];
SAM[length(intermediaryRow) + 1 : length(intermediaryCol) + length(factorRow), intermediaryCol] = IO[factorRow, intermediaryCol];
#creating vectors of titles for SAM
IOCode = append!(string.(IO[1, intermediaryCol]), string.(IO[factorRow, 1]), string.(IO[3,finalDemandCol]), string.(["Total"]));
IOName = append!(string.(IO[2, intermediaryCol]), string.(IO[factorRow, 2]), string.(IO[2,finalDemandCol]), string.(["Total"]));
#=code to sum all investment final demand into one column, comment out if not needed
#convert to dataframe, make 1st cap formation column into the sum collumn, delete other collumns, delete other rows
investment = findall(x -> occursin("Capital Formation", x), string.(IO[2,:]));
SAM = DataFrame(SAM, :auto);
names!(SAM, Symbol.(IOCode));
SAM[:,investment[1]]=sum(eachcol(SAM[:,investment]));
SAM = SAM[Not(investment[Not(investment[1])]),Not(investment[Not(investment[1])])];
#alter title vectors accordingly (include Q in total investment collumn in IOcode)
IOCode[investment[1]]
IOCode = IOCode[Not(investment[Not(investment[1])])];
IONames = IONames[Not(investment[Not(investment[1])])];
#convert back to matrix
=#
#combining capital formation collumns into relevant collumns - dynamic kinda stuff, leave for another time
#=
#lambdaPub is the government owned share of public corporations (min 0.5)
#lambdaPriv is the government owned share of private corporations (max 0.5)
@vars lambdaPub
@vars lambdaPriv
#SAM = Any[SAM]
SAM = DataFrame(SAM, :auto);
privInv = findall(x -> occursin("Private", x), IOName)[1];
pubInv = findall(x -> occursin("Public", x), IOName)[3];
govInv = findall(x -> occursin("Government", x), IOName)[2];
householdsExpend = findall(x -> occursin("Households", x), IOName)[1];
govExpend = findall(x -> occursin("Government", x), IOName)[1];
HOH = []
GOV = []
for i in 1:length(IOName);
    push!(GOV, SAM[i, govExpend]+lambdaPub*SAM[i,pubInv]+lambdaPriv*SAM[i,privInv]+SAM[i, govInv]);
    push!(HOH, SAM[i, householdsExpend]+(1-lambdaPub)*pubInv+(1-lambdaPriv)*privInv);
end
df[!, "HOH"]=HOH
=#

#SAM = SAM[Not(investment[Not(investment[1])]),Not(investment[Not(investment[1])])];
#alter title vectors accordingly (include Q in total investment collumn in IOcode)
IOCode = IOCode[Not(investment[Not(investment[1])])];
IONames = IONames[Not(investment[Not(investment[1])])];
#creating vectors for groups within SAM
intermediary = findall(x -> occursin("0", x), IOcode);
finalDemand = findall(x -> occursin("Q", x), IOcode);
factor = findall(x -> occursin("P", x), IOcode);
#steps to compute missing values
#for ring in 1:127;
 #   SAM[ring,128] = sum(SAM[ring,1:127]);
  #  SAM[128,ring] = sum(SAM[1:127,ring]);
#end
#assuming compensation and gross operating surplus all goes into households;
#SAM[122, 116:117] .= SAM[116:117, 129];
#assuming all tax goes directly to the government;
#SAM[12];

#=
#convert to dictionary
function increment!( d::Dict{S, T}, k::S, i::T) where {T<:Real, S<:Any}
    if haskey(d, k)
        d[k] += i
    else
        d[k] = i
    end
end
increment!(d::Dict{S, T}, k::S ) where {T<:Real, S<:Any} = increment!( d, k, one(T))

function df2dict( df::DataFrame, key_col::Symbol, val_col::Symbol=:null)
    keytype = typeof(df[1,key_col])
    if val_col == :null
        valtype = Int
    else
        valtype = typeof(df[1,val_col])
    end
    D = Dict{keytype, valtype}()
    for i=1:size(df,1)
        if !ismissing(df[i,key_col])
            if val_col == :null
                increment!( D, df[i,key_col] )
            elseif valtype <: Real
                increment!( D, df[i,key_col], df[i,val_col] )
            else
                if haskey(D, df[i,key_col])
                    @warn("non-unique entry: $(df[i,key_col])")
                else
                    D[df[i,key_col]] = df[i,val_col]
                end
            end
        end
    end
    return D
end
df[!, "IOCode"]=IOcode
D = df2dict(df, :IOCode, :x3)
=#


text:
module StaticTypeAnalysis

macro lintpragma(ex); end

function __init__()
    global const EQ_METHOD_FALSE = which(==, Tuple{Void, Int})
    global const CONSTRUCTOR_FALLBACK = which(Void, Tuple{Void})
end

"""
    StaticTypeAnalysis.canequal(S::Type, T::Type) :: Nullable{Bool}

Given types `S` and `T`, return `Nullable(false)` if it is not possible for
`s::S == t::T`. Return `Nullable(true)` if it is possible, and
`Nullable{Bool}()` if it cannot be determined.

```jldoctest
julia> StaticTypeAnalysis.canequal(Int, Float64)
Nullable(true)

julia> StaticTypeAnalysis.canequal(Int, String)
Nullable(false)
```
"""
function canequal(S::Type, T::Type)
    if S == Union{} || T == Union{}
        return Nullable(false)
    elseif typeintersect(S, T) ≠ Union{}
        # TODO: this is not fully correct; some types are not Union{} but still
        # not instantiated
        return Nullable{Bool}(true)
    elseif isleaftype(S) && isleaftype(T) &&
           EQ_METHOD_FALSE == which(==, Tuple{S, T})
        # == falls back to === here, but we saw earlier that the intersection
        # is empty
        return Nullable(false)
    elseif try zero(S) == zero(T) catch false end
        return Nullable{Bool}(true)
    else
        return Nullable{Bool}()
    end
end

"""
    StaticTypeAnalysis.isknownerror(f, argtypes)

Return `true` if it is known, without inference, that calling `f` with
signature `argtypes` will result in an error.

Consumers of this package are advised to use `infertype` and check the result
against `Union{}`, which covers more cases.
"""
isknownerror(_f, _argtypes) = false
function isknownerror(::typeof(Base.getindex), argtypes)
    if isempty(argtypes)
        true
    elseif argtypes[1] <: Associative
        if Base.length(argtypes) ≠ 2
            true
        else
            try
                K = keytype(argtypes[1])
                ce = canequal(K, argtypes[2])
                !isnull(ce) && !get(ce)
            catch
                false
            end
        end
    else
        false
    end
end

"""
    StaticTypeAnalysis.infertype(f, argtypes)

Given a function `f` and a list of types `argtypes`, use inference and other
static type checking techniques to figure out a type `S` such that the result
of applying `f` to `argtypes` is always of type `S`.
"""
infertype(f, argtypes) = infertype(f, (argtypes...))
function infertype(f, argtypes::Tuple)
    if isknownerror(f, argtypes)
        Union{}
    elseif f === Base.getindex && Base.length(argtypes) == 2 &&
           argtypes[1] == UnitRange && argtypes[2] <: UnitRange &&
           eltype(argtypes[2]) <: Integer
        # TODO: would be nice to get rid of this odd special case
        UnitRange
    elseif isa(f, Type) && Base.length(argtypes) == 1 &&
           isleaftype(argtypes[1]) &&
           which(f, Tuple{argtypes[1]}) === CONSTRUCTOR_FALLBACK
        # we can infer better code for the constructor `convert` fallback by
        # inferring the convert itself
        Core.Inference.return_type(convert, Tuple{Type{f}, argtypes[1]})
    else
        try
            typejoin(Base.return_types(f, Tuple{argtypes...})...)
        catch  # error might be thrown if generic function, try using inference
            if all(isleaftype, argtypes)
                Core.Inference.return_type(f, Tuple{argtypes...})
            else
                Any
            end
        end
    end
end

"""
    StaticTypeAnalysis.getindexable(T::Type)

Return `true` if all objects of type `T` support `getindex`, and the `getindex`
operation on numbers is consistent with iteration order.

Note that, in particular, this is not true for `String` and `Dict`.
"""
getindexable{T<:Union{Tuple,Pair,Array,Number}}(::Type{T}) = true
getindexable(::Type) = false

"""
    StaticTypeAnalysis.length(T::Type) :: Nullable{Int}

If it can be determined that all objects of type `T` have length `n`, then
return `Nullable(n)`. Otherwise, return `Nullable{Int}()`.
"""
length(::Type{Union{}}) = Nullable(0)
length(::Type) = Nullable{Int}()
length{T<:Pair}(::Type{T}) = Nullable(2)

if VERSION < v"0.6.0-dev.2123" # where syntax introduced by julia PR #18457
    length{T<:Tuple}(::Type{T}) = if !isa(T, DataType) || Core.Inference.isvatuple(T)
        Nullable{Int}()
    else
        Nullable{Int}(Base.length(T.types))
    end
else
    include_string("""
    length(::Type{T}) where T <: NTuple{N, Any} where N = Nullable{Int}(N)
    """)
end

"""
    StaticTypeAnalysis.eltype(T::Type)

Return `S` as specific as possible such that all objects of type `T` have
element type `S`.
"""
eltype(::Type{Union{}}) = Union{}
eltype(T::Type) = Base.eltype(T)

_getindex_nth{n}(xs::Any, ::Type{Val{n}}) = xs[n]
_typeof_nth_getindex{T}(::Type{T}, n::Integer) =
    infertype(_getindex_nth, Any[T, Type{Val{Int(n)}}])

"""
    StaticTypeAnalysis.typeof_nth(T::Type)

Return `S` as specific as possible such that all objects of type `T`, when
iterated over, have `n`th element type `S`.
"""
typeof_nth(T::Type, n::Integer) =
    if getindexable(T)
        typeintersect(eltype(T), _typeof_nth_getindex(T, n))
    else
        eltype(T)
    end
typeof_nth{K,V}(::Type{Pair{K,V}}, n::Integer) =
    n == 1 ? K : n == 2 ? V : Union{}
typeof_nth(::Type{Union{}}, ::Integer) = Union{}

end


text:
function __init__polblogs()
    LINK = "https://netset.telecom-paris.fr/datasets/polblogs.tar.gz"
    DEPNAME = "PolBlogs"

    register(DataDep(DEPNAME,
    """
    Dataset : The $DEPNAME dataset
    Website : $LINK
    """,
    LINK,
    post_fetch_method = unpack
    ))
end

"""
    PolBlogs(; dir=nothing)

The Political Blogs dataset from the [The Political Blogosphere and
the 2004 US Election: Divided they Blog](https://dl.acm.org/doi/10.1145/1134271.1134277) paper.

`PolBlogs` is a graph with 1,490 vertices (representing political blogs) and 19,025 edges (links between blogs).

The links are automatically extracted from a crawl of the front page of the blog.

Each vertex receives a label indicating the political leaning of the blog: liberal or conservative.
"""
struct PolBlogs <: AbstractDataset
    metadata::Dict{String, Any}
    graphs::Vector{Graph}
end

function PolBlogs(; dir=nothing)
    DEPNAME = "PolBlogs"
    DATA = ["adjacency.csv", "labels.csv"]

    path = datafile(DEPNAME, DATA[1], dir)
    adj = 1 .+ Matrix{Int64}(readdlm(path, ','))
    s, t = adj[:,1], adj[:,2]

    path = datafile(DEPNAME, DATA[2], dir)
    labels = Matrix{Int64}(readdlm(path, ',')) |> vec

    metadata = Dict{String, Any}()
    g = Graph(; num_nodes=1490, num_edges=19025,
            edge_index = (s, t),
            node_data = (; labels))

    return PolBlogs(metadata, [g])
end

Base.length(d::PolBlogs) = length(d.graphs)
Base.getindex(d::PolBlogs) = d.graphs
Base.getindex(d::PolBlogs, i) = getindex(d.graphs, i)


text:
fname = path*"/SampleFiles/fdsn.conf"
sac_pz_file = path*"/SampleFiles/SAC/test_sac.pz"
hood_reg = Float64[44.8, 46.0, -122.4, -121.0]
rainier_rad = Float64[46.852886, -121.760374, 0.0, 0.1]

printstyled("  FDSN web requests\n", color=:light_green)

# FDSNsta
printstyled("    FDSNsta\n", color=:light_green)
S = FDSNsta("CC.VALT..,PB.B001..BS?,PB.B001..E??")
for i in ("PB.B001.T0.BS1", "PB.B001..EHZ", "CC.VALT..BHZ")
  j = findid(S, i)
  try
    @test j > 0
  catch
    @warn(string("No data from ", i, "; check connection!"))
  end
end

# FDSNsta with MultiStageResp
S = FDSNsta("CC.VALT..,PB.B001..BS?,PB.B001..E??", msr=true)

# With autoname
printstyled("      get_data(\"FDSN\", ..., autoname=true)\n", color=:light_green)
req_f = "2019.001.00.00.00.000.UW.VLL..EHZ.R.mseed"
req_ok = (try
    S = get_data("FDSN", "UW.VLL..EHZ", src="IRIS", s="2019-01-01", t=3600, autoname=true)
    true
  catch
    @warn("Station VLL appears to be offline; test skipped.")
    false
  end)
if req_ok
  @test safe_isfile(req_f)
  rm(req_f)

  printstyled("        test match to IRIS filename convention\n", color=:light_green)
  S = get_data("IRIS", "UW.VLL..EHZ", s="2019-01-01", t=3600, autoname=true)
  @test safe_isfile(req_f)
end

printstyled("      radius search (rad=)\n", color=:light_green)
rad = Float64[45.373514, -121.695919, 0.0, 0.1]
S = FDSNsta(rad=rainier_rad)
try
  @test S.n > 0 # Test will break if everything around Mt. Rainier is offline
catch
  @warn("Stations around Mt. Rainier appear to be offline.")
end

printstyled("      rectangular search (reg=)\n", color=:light_green)
S = FDSNsta(reg=hood_reg)
try
  @test S.n > 0 # Test will break if everything around Mt. Hood is offline
catch
  @warn("Stations around Mt. Rainier appear to be offline.")
end
printstyled("    get_data\n", color=:light_green)

printstyled("      GeoCSV output\n", color=:light_green)
S = get_data("FDSN", "CC.JRO..BHZ,IU.COLA.00.*", src="IRIS", s=-600, t=0, fmt="geocsv", w=true)
S = get_data("FDSN", "CC.JRO..BHZ,CC.VALT.*", src="IRIS", s=-300, t=0, fmt="geocsv.slist")

printstyled("      config file for channel spec\n", color=:light_green)

S = SeisData()
get_data!(S, "FDSN", fname, src="IRIS", s=-600, t=0, w=true)
deleteat!(S, findall(S.fs.<25.0))
filtfilt!(S, fl=0.01, fh=10.0)

# Ensure station headers are set
ids = ["UW.HOOD..ENE", "CC.VALT..BHZ", "UW.TDH..EHZ", "UW.VLL.EHZ"]
fss = [100.0, 50.0, 100.0, 100.0]
codes = ['N', 'H', 'H', 'H']
for i = 1:4
  j = findid(S, ids[i])
  if j > 0
    @test ≈(S.fs[j], fss[i])
    @test inst_code(S, j) == codes[i]
    break
  else
    @warn(string("No data from ", ids[i], "; check connection!"))
  end
end

# Check that headers get overwritten with SACPZ info when we use read_sacpz
Nc = S.n
read_sacpz!(S, sac_pz_file)
@test S.n > Nc
i = findid("CC.VALT..BHZ", S)
if i > 0
  @test S.misc[i]["OUTPUT UNIT"] == "COUNTS"
end
i = findid("UW.HOOD..ENE", S)
if i > 0
  @test S.misc[i]["INSTTYPE"] == "ES-T-3339=Q330S+-6410"
end

# Ensure we got data
L = [length(x) for x in S.x]
if isempty(L) == false
  @test (maximum(L) > 0)
end

# Check that msr=true works
S = SeisData()
get_data!(S, "FDSN", fname, src="IRIS", msr=true, s=-600, t=0)
if isempty(S)
  warn("Empty request; check connectivity!")
else
  for i in 1:S.n
    @test typeof(S.resp[i]) == MultiStageResp
  end
end

# Try a string array for input
printstyled("      string array for channel spec\n", color=:light_green)
S = SeisData()
get_data!(S, "FDSN", ["UW.HOOD..E??", "CC.VALT..???", "UW.XNXNX.99.QQQ"], src="IRIS", s=-600, t=0, opts="szsrecs=true")

# Try a single string
printstyled("      string for channel spec\n", color=:light_green)
S = get_data("FDSN", "CC.JRO..BHZ,IU.COLA.00.*", src="IRIS", s=-600, t=0, v=1,
  demean=true,
  detrend=true,
  rr=true,
  taper=true,
  ungap=true,
  unscale=true)

# This should return exactly 4 days of data, which we know IRIS' FDSN server has
printstyled("      multi-day request\n", color=:light_green)
ts = "2018-01-31T00:00:00"
te = "2018-02-02T00:00:00"
S = get_data("FDSN","CI.ADO..BH?", s=ts, t=te)
id = "CI.ADO..BHE"
i = findid(S, id)
if i == 0
  @warn(string("No data from ", id, "; check connection!"))
elseif isempty(S)
  @warn(string("Request empty; rest of test skipped!"))
else
  # Check that we have two complete days of data with no gaps

  if (length(S.x[i]) / (86400*S.fs[i])) != 2.0
    @warn(string("Partial outage; missing data from ", id, "; check connection!"))
  end

  printstyled("        are data written identically?\n", color=:green)

  # write ASDF first, since this modifies the first sample start time in S
  write_hdf5("sacreq.h5", S, add=true, ovr=true, len=Day(2))
  S1 = read_hdf5("sacreq.h5", ts, te, msr=false)

  for f in (:id, :loc, :fs, :gain, :resp, :units, :misc)
    @test isequal(getfield(S1,f), getfield(S,f))
  end
  for i in 1:S.n
    x1 = S1.x[i]
    x2 = S.x[i]
    if (any(isnan.(x1)) == false) && (any(isnan.(x2)) == false)
      @test x1 ≈ x2
    end
  end

  # Check that these data can be written and read faithfully in ASDF, SAC, and SeisIO
  writesac(S)
  wseis("sacreq.seis", S)
  S2 = rseis("sacreq.seis")[1]
  for f in (:id, :loc, :fs, :gain, :resp, :units, :misc)
    @test isequal(getfield(S2,f), getfield(S,f))
  end
  for i in 1:S.n
    x1 = S.x[i]
    x2 = S2.x[i]
    if (any(isnan.(x1)) == false) && (any(isnan.(x2)) == false)
      @test x1 ≈ x2
    end
  end

  # These are the only fields preserved; :loc is preserved to Float32 precision
  d0 = DateTime(ts)
  y0 = Year(d0).value
  j0 = md2j(y0, Month(d0).value, Day(d0).value)
  sac_str = string(y0, ".", lpad(j0, 3, '0'), "*CI.ADO..BH*SAC")
  S1 = read_data("sac", sac_str)
  for f in (:id, :fs, :gain, :t)
    @test getfield(S, f) == getfield(S1, f)
    @test getfield(S, f) == getfield(S2, f)
  end
  for f in (:lat, :lon, :el, :dep, :az, :inc)
    @test isapprox(getfield(S.loc[i], f), getfield(S1.loc[i], f), atol=1.0e-3)
  end
  for i in 1:S.n
    x = S.x[i]
    x1 = S1.x[i]
    x2 = S2.x[i]
    if (any(isnan.(x)) == false) && (any(isnan.(x1)) == false) && (any(isnan.(x2)) == false)
      @test x ≈ x1
      @test x1 ≈ x2
    end
  end

  # clean up
  rm("sacreq.seis")
  rm("sacreq.h5")
end

# A bad data format should produce a warning
printstyled("      request an unparseable format (sac.zip)\n", color=:light_green)

redirect_stdout(out) do
  get_data!(S, "FDSN", "UW.LON.."; src="IRIS", s=-600, t=0, v=3, fmt="sac.zip")
end

# Potsdam test
printstyled("      request from GFZ\n", color=:light_green)
S = get_data("FDSN", "GE.BKB..BH?", src="GFZ", s="2011-03-11T06:00:00", t="2011-03-11T06:05:00", v=0, y=false)
if isempty(S)
  @warn(string("No data from GFZ request; check connection!"))
end

# ❄❄❄❄❄❄❄❄❄❄❄❄❄❄❄❄❄❄ (oh, California...)
printstyled("    servers with special headers:\n", color=:light_green)
ds = now()-Day(1)
ds -= Millisecond(ds)
s = string(ds)
t = string(ds+Hour(1))

rubric = [
  "NCEDC" "BK.MOD..BHE"
  "SCEDC" "CI.SDD..BHZ"
]
❄ = size(rubric, 1)

for i = 1:❄
  printstyled("      ", rubric[i,1], ":\n", color=:light_green)
  printstyled("        station info\n", color=:light_green)
  S = FDSNsta(rubric[i,2], s=s, t=t, msr=true, src=rubric[i,1])
  if isempty(S)
    printstyled("        No data; check headers & connection!\n", color=:red)
  end

  if rubric[i] != "SCEDC"
    printstyled("        trace data\n", color=:light_green)
    S = get_data("FDSN", rubric[i,2], src=rubric[i,1], s=s, t=t, msr=true, w=true)
    if isempty(S)
      printstyled("        No data; check headers & connection!\n", color=:red)
    end
  end
end


text:
using Word2Vec
using Base.Test

data_dir = joinpath(Pkg.dir("Word2Vec"), "test", "data")
test_file = joinpath(data_dir, "text8_tiny")

printover(i) = print("\r$i\e[K")

function test_word_window()
    open(test_file) do fp
        test_word_window(fp)
    end
end

function test_word_window(fp)
    println("Testing words_of...")

    for i in Word2Vec.words_of(fp)
        printover(i)
    end

    seekstart(fp)
    for i in map(uppercase, Word2Vec.words_of(fp))
        printover(i)
    end

    println("\rDone.\e[K")

    println("Testing sliding_window...")
    seekstart(fp)
    for i in Word2Vec.sliding_window(Word2Vec.words_of(fp))
        printover(join(i, ','))
    end
    println("\nDone.")
end

test_word_window()


text:
function additional_prices(mdl::Union{Model,ModelInit})

    R,Rᴸ = mdl.Prices
    ω = (Lₛ * ( 1 + (1-ρ) * R) )^-1
    τ = (ϐ / ω) * (1 - Lₛ) / Lₛ

    return([ω,τ])

end


text:
n = 10
A = Diagonal(range(1, stop=2, length=n))
f(x) = dot(x,A*x)/2
g(x) = A*x
g!(stor,x) = copyto!(stor,g(x))
x0 = randn(n)

manif = Optim.Sphere()
Optim.optimize(f, g!, x0, Optim.ConjugateGradient(manifold=manif))

n = 10
A = Diagonal(range(1, stop=2, length=n))
f(x) = dot(x[1:n],A*x[1:n])/2 + dot(x[n+1:2n],A*x[n+1:2n])/2
x0 = randn(2n)

manif = Optim.ProductManifold(Optim.Sphere(), Optim.Sphere(), (n,), (n,))
res = Optim.optimize(f, x0, Optim.ConjugateGradient(manifold=manif))

v =Optim.minimizer(res)

norm(v[1:n])
norm(v[n+1:end])

norm(v)

#### iterative manifold ####

# r column vectors noramlized of dimension n

n = 10
r = n

manif = Optim.Sphere()
for layer = 1:r-1
    manif = Optim.ProductManifold(Optim.Sphere(), manif, (n,), (layer * n,))
end

function array_to_matrix(array, r, n)

    """gives a (r,n) matrix from the r, n-sized vectors concatenated in array"""

    mat = Matrix{eltype(array)}(undef, r, n)

    for line = 1 : r
        for col = 1 : n
            mat[line, col] = array[(line-1)*n+col]
        end
    end

    mat

end

function objective(array, r, n)

    M = array_to_matrix(array, r, n)

    gram_matrix = M'*M

    -abs(permanent_ryser(gram_matrix))

end

function gradient_permanent!(G, array, n, r)

    A = array_to_matrix(array, r, n)

    for i = 1:n
        for j = 1:n
            grad[i,j] = permanent_ryser(remove_row_col(A, i, j))
        end
    end

end

function foo(array)
    objective(array,r,n)
end

function g!(G, array)
    gradient_permanent!(G, array, n, r)
end

res = Optim.optimize(foo,rand_gram_matrix(n), Optim.GradientDescent(manifold=manif), Optim.Options(x_tol = 1e-16, iterations = 100000))
Float64(factorial(n))
v = Optim.minimizer(res)

println("###############")
for i = 1 : n
    println(real(norm(v[i,:])))
end


text:
# Horizontal and Vertical Recurrance Relationship of Head-Gordon and Pople
#
# A method for two-electron Gaussian integral and integral derivative evaluation using recurrence relations
# Martin Head-Gordon and John A. Pople
# Department of Chemistry, Carnegie-Mellon University, Pittsburgh, Pennsylvania 15213 (Received 26 May 1988; accepted 20 July 1988)

# An efficient method is presented for evaluating two-electron Cartesian
# Gaussian integrals, and their first derivatives with respect to nuclear
# coordinates. It is based on the recurrence relation (RR) of Obara and Saika
# [J. Chem. Phys. 84, 3963 (1986)], and an additional new RR, which are
# combined together in a general algorithm applicable to any angular momenta.
# This algorithm exploits the fact that the new RR can be applied outside
# contraction loops. It is shown, by floating point operation counts and
# comparative timings, to be generally superior to existing methods,
# particularly for basis sets containing d functions.
#  https://doi.org/10.1063/1.455553

function hrr(aexpn::Float64,ax::Float64,ay::Float64,az::Float64,aI::Int64,aJ::Int64,aK::Int64,
    bexpn::Float64,bx::Float64,by::Float64,bz::Float64,bI::Int64,bJ::Int64,bK::Int64,
    cexpn::Float64,cx::Float64,cy::Float64,cz::Float64,cI::Int64,cJ::Int64,cK::Int64,
    dexpn::Float64,dx::Float64,dy::Float64,dz::Float64,dI::Int64,dJ::Int64,dK::Int64)
    if bI > 0
        return hrr(aexpn,ax,ay,az,aI+1,aJ,aK,bexpn,bx,by,bz,bI-1,bJ,bK,
            cexpn,cx,cy,cz,cI,cJ,cK,dexpn,dx,dy,dz,dI,dJ,dK) +
        (ax-bx)*hrr(aexpn,ax,ay,az,aI,aJ,aK,bexpn,bx,by,bz,bI-1,bJ,bK,
                cexpn,cx,cy,cz,cI,cJ,cK,dexpn,dx,dy,dz,dI,dJ,dK)
    elseif bJ > 0
        return hrr(aexpn,ax,ay,az,aI,aJ+1,aK,bexpn,bx,by,bz,bI,bJ-1,bK,
            cexpn,cx,cy,cz,cI,cJ,cK,dexpn,dx,dy,dz,dI,dJ,dK) +
            (ay-by)*hrr(aexpn,ax,ay,az,aI,aJ,aK,bexpn,bx,by,bz,bI,bJ-1,bK,
                cexpn,cx,cy,cz,cI,cJ,cK,dexpn,dx,dy,dz,dI,dJ,dK)
    elseif bK > 0
        return hrr(aexpn,ax,ay,az,aI,aJ,aK+1,bexpn,bx,by,bz,bI,bJ,bK-1,
            cexpn,cx,cy,cz,cI,cJ,cK,dexpn,dx,dy,dz,dI,dJ,dK) +
            (az-bz)*hrr(aexpn,ax,ay,az,aI,aJ,aK,bexpn,bx,by,bz,bI,bJ,bK-1,
                cexpn,cx,cy,cz,cI,cJ,cK,dexpn,dx,dy,dz,dI,dJ,dK)
    elseif dI > 0
        return hrr(aexpn,ax,ay,az,aI,aJ,aK,bexpn,bx,by,bz,bI,bJ,bK,
            cexpn,cx,cy,cz,cI+1,cJ,cK,dexpn,dx,dy,dz,dI-1,dJ,dK) +
            (cx-dx)*hrr(aexpn,ax,ay,az,aI,aJ,aK,bexpn,bx,by,bz,bI,bJ,bK,
                cexpn,cx,cy,cz,cI,cJ,cK,dexpn,dx,dy,dz,dI-1,dJ,dK)
    elseif dJ > 0
        return hrr(aexpn,ax,ay,az,aI,aJ,aK,bexpn,bx,by,bz,bI,bJ,bK,
            cexpn,cx,cy,cz,cI,cJ+1,cK,dexpn,dx,dy,dz,dI,dJ-1,dK) +
            (cy-dy)*hrr(aexpn,ax,ay,az,aI,aJ,aK,bexpn,bx,by,bz,bI,bJ,bK,
                cexpn,cx,cy,cz,cI,cJ,cK,dexpn,dx,dy,dz,dI,dJ-1,dK)
    elseif dK > 0
        return hrr(aexpn,ax,ay,az,aI,aJ,aK,bexpn,bx,by,bz,bI,bJ,bK,
            cexpn,cx,cy,cz,cI,cJ,cK+1,dexpn,dx,dy,dz,dI,dJ,dK-1) +
            (cz-dz)*hrr(aexpn,ax,ay,az,aI,aJ,aK,bexpn,bx,by,bz,bI,bJ,bK,
                cexpn,cx,cy,cz,cI,cJ,cK,dexpn,dx,dy,dz,dI,dJ,dK-1)
    end
    return vrr(aexpn,ax,ay,az,aI,aJ,aK,bexpn,bx,by,bz,cexpn,cx,cy,cz,cI,cJ,cK,dexpn,dx,dy,dz,0)
end

"""
Vertical Recurrence Relationship (VRR)

"""
function vrr(aexpn::Float64,ax::Float64,ay::Float64,az::Float64,aI::Int64,aJ::Int64,aK::Int64,
        bexpn::Float64,bx::Float64,by::Float64,bz::Float64,
        cexpn::Float64,cx::Float64,cy::Float64,cz::Float64,cI::Int64,cJ::Int64,cK::Int64,
        dexpn::Float64,dx::Float64,dy::Float64,dz::Float64,m::Int64)
    px,py,pz = gaussian_product_center(aexpn,ax,ay,az,bexpn,bx,by,bz)
    qx,qy,qz = gaussian_product_center(cexpn,cx,cy,cz,dexpn,dx,dy,dz)
    zeta,eta = aexpn+bexpn,cexpn+dexpn
    wx,wy,wz = gaussian_product_center(zeta,px,py,pz,eta,qx,qy,qz)
    #println("P: $px,$py,$pz, Q: $qx,$qy,$qz, W: $wx,$wy,$wz, $zeta,$eta")

    val::Float64 = 0.0
    if cK>0
        val = (qz-cz)*vrr(aexpn,ax,ay,az,aI,aJ,aK,bexpn,bx,by,bz,
            cexpn,cx,cy,cz,cI,cJ,cK-1,dexpn,dx,dy,dz,m) +
            (wz-qz)*vrr(aexpn,ax,ay,az,aI,aJ,aK,bexpn,bx,by,bz,
                cexpn,cx,cy,cz,cI,cJ,cK-1,dexpn,dx,dy,dz,m+1)
        #println("val1=$val")
        if cK>1
            val += 0.5*(cK-1)/eta*(
                vrr(aexpn,ax,ay,az,aI,aJ,aK,bexpn,bx,by,bz,
                    cexpn,cx,cy,cz,cI,cJ,cK-2,dexpn,dx,dy,dz,m) -
            zeta/(zeta+eta)*
                vrr(aexpn,ax,ay,az,aI,aJ,aK,bexpn,bx,by,bz,
                    cexpn,cx,cy,cz,cI,cJ,cK-2,dexpn,dx,dy,dz,m+1) )
        #println("val2=$val")
        end
        if aK>0
            val += 0.5*aK/(zeta+eta)*
                vrr(aexpn,ax,ay,az,aI,aJ,aK-1,bexpn,bx,by,bz,
                    cexpn,cx,cy,cz,cI,cJ,cK-1,dexpn,dx,dy,dz,m+1)
        end
        #println("val3=$val")
        return val
    elseif cJ>0
        val = (qy-cy)*vrr(aexpn,ax,ay,az,aI,aJ,aK,bexpn,bx,by,bz,
            cexpn,cx,cy,cz,cI,cJ-1,cK,dexpn,dx,dy,dz,m) +
        (wy-qy)*vrr(aexpn,ax,ay,az,aI,aJ,aK,bexpn,bx,by,bz,
                cexpn,cx,cy,cz,cI,cJ-1,cK-1,dexpn,dx,dy,dz,m+1)
        #println("val4=$val")
        if cJ>1
            val += 0.5*(cJ-1)/eta*(
            vrr(aexpn,ax,ay,az,aI,aJ,aK,bexpn,bx,by,bz,
                cexpn,cx,cy,cz,cI,cJ-2,cK,dexpn,dx,dy,dz,m) -
            zeta/(zeta+eta)*
            vrr(aexpn,ax,ay,az,aI,aJ,aK,bexpn,bx,by,bz,
                cexpn,cx,cy,cz,cI,cJ-2,cK,dexpn,dx,dy,dz,m+1)
            )
        #println("val5=$val")
        end
        if aJ>0
            val += 0.5*aJ/(zeta+eta)*
            vrr(aexpn,ax,ay,az,aI,aJ-1,aK,bexpn,bx,by,bz,
                cexpn,cx,cy,cz,cI,cJ-1,cK,dexpn,dx,dy,dz,m+1)
        end
        #println("val6=$val")
        return val
    elseif cI>0
        val = (qx-cx)*vrr(aexpn,ax,ay,az,aI,aJ,aK,bexpn,bx,by,bz,
            cexpn,cx,cy,cz,cI-1,cJ,cK,dexpn,dx,dy,dz,m) +
        (wx-qx)*vrr(aexpn,ax,ay,az,aI,aJ,aK,bexpn,bx,by,bz,
                cexpn,cx,cy,cz,cI-1,cJ,cK-1,dexpn,dx,dy,dz,m+1)
        #println("val7=$val")
        if cI>1
            val += 0.5*(cI-1)/eta*(
            vrr(aexpn,ax,ay,az,aI,aJ,aK,bexpn,bx,by,bz,
                cexpn,cx,cy,cz,cI-2,cJ,cK,dexpn,dx,dy,dz,m) -
            zeta/(zeta+eta)*
            vrr(aexpn,ax,ay,az,aI,aJ,aK,bexpn,bx,by,bz,
                cexpn,cx,cy,cz,cI-2,cJ,cK,dexpn,dx,dy,dz,m+1)
            )
        end
        #println("val8=$val")
        if aI>0
            val += 0.5*aI/(zeta+eta)*
            vrr(aexpn,ax,ay,az,aI-1,aJ,aK,bexpn,bx,by,bz,
                cexpn,cx,cy,cz,cI-1,cJ,cK,dexpn,dx,dy,dz,m+1)
        end
        #println("val9=$val")
        return val
    elseif aK>0
        val = (pz-az)*vrr(aexpn,ax,ay,az,aI,aJ,aK-1,bexpn,bx,by,bz,
                cexpn,cx,cy,cz,cI,cJ,cK,dexpn,dx,dy,dz,m) +
        (wz-pz)*vrr(aexpn,ax,ay,az,aI,aJ,aK-1,bexpn,bx,by,bz,
                cexpn,cx,cy,cz,cI,cJ,cK,dexpn,dx,dy,dz,m+1)
        #println("val10=$val")
        if aK>1
            val += 0.5*(aK-1)/zeta*(
            vrr(aexpn,ax,ay,az,aI,aJ,aK-2,bexpn,bx,by,bz,
                cexpn,cx,cy,cz,cI-1,cJ,cK,dexpn,dx,dy,dz,m) -
            eta/(zeta+eta)*
            vrr(aexpn,ax,ay,az,aI,aJ,aK-2,bexpn,bx,by,bz,
                cexpn,cx,cy,cz,cI-1,cJ,cK,dexpn,dx,dy,dz,m+1)
            )
        end
        #println("val11=$val")
        return val
    elseif aJ>0
        val = (py-ay)*vrr(aexpn,ax,ay,az,aI,aJ-1,aK,bexpn,bx,by,bz,
                cexpn,cx,cy,cz,cI,cJ,cK,dexpn,dx,dy,dz,m)+
        (wy-py)*vrr(aexpn,ax,ay,az,aI,aJ-1,aK,bexpn,bx,by,bz,
                cexpn,cx,cy,cz,cI,cJ,cK,dexpn,dx,dy,dz,m+1)
        #println("val12=$val")
        if aJ>1
            val += 0.5*(aJ-1)/zeta*(
            vrr(aexpn,ax,ay,az,aI,aJ-2,aK,bexpn,bx,by,bz,
                cexpn,cx,cy,cz,cI,cJ,cK,dexpn,dx,dy,dz,m) -
            eta/(zeta+eta)*
            vrr(aexpn,ax,ay,az,aI,aJ-2,aK,bexpn,bx,by,bz,
                cexpn,cx,cy,cz,cI,cJ,cK,dexpn,dx,dy,dz,m+1)
            )
        end
        #println("val13=$val")
        return val
    elseif aI>0
        val = (px-ax)*vrr(aexpn,ax,ay,az,aI-1,aJ,aK,bexpn,bx,by,bz,
                cexpn,cx,cy,cz,cI,cJ,cK,dexpn,dx,dy,dz,m) +
        (wx-px)*vrr(aexpn,ax,ay,az,aI-1,aJ,aK,bexpn,bx,by,bz,
                cexpn,cx,cy,cz,cI,cJ,cK,dexpn,dx,dy,dz,m+1)
        #println("val14=$val")
        if aI>1
            val += 0.5*(aI-1)/zeta*(
            vrr(aexpn,ax,ay,az,aI-2,aJ,aK,bexpn,bx,by,bz,
                cexpn,cx,cy,cz,cI,cJ,cK,dexpn,dx,dy,dz,m) -
            eta/(zeta+eta)*
            vrr(aexpn,ax,ay,az,aI-2,aJ,aK,bexpn,bx,by,bz,
                cexpn,cx,cy,cz,cI,cJ,cK,dexpn,dx,dy,dz,m+1)
            )
        end
        #println("val15=$val")
        return val
    end

    rab2 = dist2(ax-bx,ay-by,az-bz)
    rcd2 = dist2(cx-dx,cy-dy,cz-dz)
    rpq2 = dist2(px-qx,py-qy,pz-qz)
    T = zeta*eta/(zeta+eta)*rpq2
    Kab = sqrt(2)pi^1.25/zeta*exp(-aexpn*bexpn*rab2/zeta)
    Kcd = sqrt(2)pi^1.25/eta*exp(-cexpn*dexpn*rcd2/eta)
    #println("rab2=$rab2,rcd2=$rcd2,rpq2=$rpq2,T=$T,Kab=$Kab,Kcd=$Kcd")
    return Kab*Kcd/sqrt(zeta+eta)*Fgamma(m,T)
end

function vrr_iter(aexpn::Float64,ax::Float64,ay::Float64,az::Float64,aI::Int64,aJ::Int64,aK::Int64,
        bexpn::Float64,bx::Float64,by::Float64,bz::Float64,
        cexpn::Float64,cx::Float64,cy::Float64,cz::Float64,cI::Int64,cJ::Int64,cK::Int64,
        dexpn::Float64,dx::Float64,dy::Float64,dz::Float64,M::Int64)
    px,py,pz = gaussian_product_center(aexpn,ax,ay,az,bexpn,bx,by,bz)
    qx,qy,qz = gaussian_product_center(cexpn,cx,cy,cz,dexpn,dx,dy,dz)
    zeta,eta = aexpn+bexpn,cexpn+dexpn
    wx,wy,wz = gaussian_product_center(zeta,px,py,pz,eta,qx,qy,qz)
    rab2 = dist2(ax-bx,ay-by,az-bz)
    rcd2 = dist2(cx-dx,cy-dy,cz-dz)
    rpq2 = dist2(px-qx,py-qy,pz-qz)
    T = zeta*eta/(zeta+eta)*rpq2
    Kab = sqrt(2)pi^1.25/zeta*exp(-aexpn*bexpn*rab2/zeta)
    Kcd = sqrt(2)pi^1.25/eta*exp(-cexpn*dexpn*rcd2/eta)
    mtot = aI+aJ+aK+cI+cJ+cK+M

    vrr_terms = zeros(Float64,(aI+1,aJ+1,aK+1,cI+1,cJ+1,cK+1,mtot+1))

    for m in 0:mtot
        vrr_terms[1,1,1, 1,1,1, m+1] = Fgamma(m,T)*Kab*Kcd/sqrt(zeta+eta)
    end

    for i in 0:(aI-1)
        for m in 0:(mtot-i-1)
            vrr_terms[i+2,1,1, 1,1,1, m+1] = (
                 (px-ax)*vrr_terms[i+1,1,1, 1,1,1, m+1] +
                 (wx-px)*vrr_terms[i+1,1,1, 1,1,1, m+2])

            if i>0
                vrr_terms[i+2,1,1, 1,1,1, m+1] += i/2/zeta*(
                    vrr_terms[i,1,1, 1,1,1, m+1] -
                    eta/(zeta+eta)*vrr_terms[i,1,1, 1,1,1, m+2])
            end
        end
    end
#=
    for j in 0:(aJ-1)
        for i in 0:aI
            for m in 0:(mtot-i-j-1)
                println(("b",i,j,m))
                vrr_terms[i+1,j+2,1, 1,1,1, m+1] = (
                (py-ay)*vrr_terms[i+1,j+1,1, 1,1,1, m+1] +
                (wy-py)*vrr_terms[i+1,j+1,1, 1,1,1, m+2])
                if j>0
                    vrr_terms[i+1,j+2,1, 1,1,1, m+1] += j/2/zeta*(
                        vrr_terms[i+1,j,1, 1,1,1, m+1] -
                        eta/(zeta+eta)*vrr_terms[i+1,j,1, 1,1,1, m+2])
                end
            end
        end
    end
    for k in 0:(aK-1)
        for j in 0:aJ
            for i in 0:aI
                for m in 0:(mtot-i-j-k-1)
                    println(("c",i,j,k,m))
                    vrr_terms[i+1,j+1,k+2, 1,1,1, m+1] = (
                    (pz-az)*vrr_terms[i+1,j+1,k+1, 1,1,1, m+1] +
                    (wz-pz)*vrr_terms[i+1,j+1,k+1, 1,1,1, m+2])
                    if k>0
                        vrr_terms[i+1,j+1,k+2, 1,1,1, m+1] += k/2/zeta*(
                        vrr_terms[i+1,j+1,k, 1,1,1, m+1] -
                        eta/(zeta+eta)*vrr_terms[i+1,j+1,k, 1,1,1, m+2])
                    end
                end
            end
        end
    end
    for q in 0:(cI-1)
        for k in 0:aK
            for j in 0:aJ
                for i in 0:aI
                    for m in 0:(mtot-i-j-k-q-1)
                        println(("d",i,j,k,q,m))
                        vrr_terms[i+1,j+1,k+1, q+2,1,1, m+1] = (
                        (qx-cx)*vrr_terms[i+1,j+1,k+1, q+1,1,1, m+1] +
                        (wx-qx)*vrr_terms[i+1,j+1,k+1, q+1,1,1, m+2])
                        if q>0
                            vrr_terms[i+1,j+1,k+1, q+2,1,1, m+1] += q/2/eta*(
                            vrr_terms[i+1,j+1,k+1, q,1,1, m+1] -
                            eta/(zeta+eta)*vrr_terms[i+1,j+1,k+1, q,1,1, m+2])
                        end
                        if i>0
                            vrr_terms[i+1,j+1,k+1, q+2,1,1, m+1] += (
                            i/2/(zeta+eta)*vrr_terms[i,j+1,j+1, q+1,1,1, m+2])
                        end
                    end
                end
            end
        end
    end
    for r in 0:(cJ-1)
        for q in 0:cI
            for k in 0:aK
                for j in 0:aJ
                    for i in 0:aI
                        for m in 0:(mtot-i-j-k-q-r-1)
                            println(("e",i,j,k,q,r,m))
                            vrr_terms[i+1,j+1,k+1, q+1,r+2,1, m+1] = (
                            (qy-cy)*vrr_terms[i+1,j+1,k+1, q+1,r+1,1, m+1] +
                            (wy-qy)*vrr_terms[i+1,j+1,k+1, q+1,r+1,1, m+2])
                            if r>0
                                vrr_terms[i+1,j+1,k+1, q+1,r+2,1, m+1] += r/2/eta*(
                                vrr_terms[i+1,j+1,k+1, q+1,r+1,1, m+1] -
                                zeta/(zeta+eta)*vrr_terms[i+1,j+1,k+1, q+1,r,1, m+2])
                            end
                            if j>0
                                vrr_terms[i+1,j+1,k+1, q+1,r+2,1, m+1] += (
                                j/2/(zeta+eta)*vrr_terms[i+1,j,k+1, q+1,r+1,1, m+2])
                            end
                        end
                    end
                end
            end
        end
    end
    for s in 0:(cK-1)
        for r in 0:cJ
            for q in 0:cI
                for k in 0:aK
                    for j in 0:aJ
                        for i in 0:aI
                            for m in 0:(mtot-i-j-k-q-r-s-1)
                                println(("f",i,j,k,q,r,s,m))
                                vrr_terms[i+1,j+1,k+1, q+1,r+1,s+2, m+1] = (
                                (qz-cz)*vrr_terms[i+1,j+1,k+1, q+1,r+1,s+1, m+1]+
                                (wz-qz)*vrr_terms[i+1,j+1,k+1, q+1,r+1,s+1, m+2])
                                if s>0
                                    vrr_terms[i+1,j+1,k+1, q+1,r+1,s+2,m+1] += s/2/eta*(
                                    vrr_terms[i+1,j+1,k+1,q+1,r+1,s, m+1] -
                                    zeta/(zeta+eta)*vrr_terms[i+1,j+1,k+1,q+1,r+1,s,m+2])
                                end
                                if k>0
                                    vrr_terms[i+1,j+1,k+1, q+1,r+1,s+2,m+1] += (
                                    k/2/(zeta+eta)*vrr_terms[i+1,j+1,k, q+1,r+1,s+1,m+2])
                                end
                            end
                        end
                    end
                end
            end
        end
    end
    println("before return")
    =#
    vrr_terms[aI+1,aJ+1,aK+1,cI+1,cJ+1,cK+1,M+1]
end

function test_vrr()
    @testset "test_vrr" begin
    ax=ay=az=bx=by=bz=cx=cy=cz=dx=dy=dz=0.0
    aexpn=bexpn=cexpn=dexpn=1.0
    aI=aJ=aK=0
    cI=cJ=cK=0
    M=0

    for (ax,ay,az, aI,aJ,aK, cI,cJ,cK, result) in [
            (0.,0.,0., 0,0,0, 0,0,0, 4.37335456733),
            (0.,0.,0., 1,0,0, 1,0,0, 0.182223107579),
            (0.,0.,0., 0,1,0, 0,1,0, 0.182223107579),
            (0.,0.,0., 0,0,1, 0,0,1, 0.182223107579),

            (0.,0.,0., 2,0,0, 2,0,0,  0.223223306785),
            (0.,0.,0., 0,2,0, 0,2,0,  0.223223306785),
            (0.,0.,0., 0,0,2, 0,0,2,  0.223223306785),

            (1.,2.,3., 1,0,0, 1,0,0, -5.63387712455e-06),
            (1.,2.,3., 0,1,0, 0,1,0, -0.000116463120359),
            (1.,2.,3., 0,0,1, 0,0,1, -0.000301178525749),

            (1.,2.,3., 2,0,0, 2,0,0, 0.00022503308545040895),
            (1.,2.,3., 0,2,0, 0,2,0, 0.0006102470883881907),
            (1.,2.,3., 0,0,2, 0,0,2, 0.0013427831014563411),

            (0.,0.,0., 1,1,0, 1,1,0, 0.0136667330685),
            (0.,0.,0., 0,1,1, 0,1,1, 0.0136667330685),
            (0.,0.,0., 1,0,1, 1,0,1, 0.0136667330685),

            (3.,2.,1., 1,1,0, 1,1,0, 5.976771621486971e-5),
            (3.,2.,1., 0,1,1, 0,1,1, 1.5742904443905067e-6),
            (3.,2.,1., 1,0,1, 1,0,1, 4.00292848649699e-6)
        ]

        val1 = vrr(aexpn,ax,ay,az,aI,aJ,aK,bexpn,bx,by,bz,
            cexpn,cx,cy,cz,cI,cJ,cK,dexpn,dx,dy,dz,M)
        val2 = vrr(cexpn,cx,cy,cz,cI,cJ,cK,dexpn,dx,dy,dz,
            aexpn,ax,ay,az,aI,aJ,aK,bexpn,bx,by,bz,M)
        @test val1 ≈ val2
        @test val1 ≈ result
        val3 = vrr_iter(aexpn,ax,ay,az,aI,aJ,aK,bexpn,bx,by,bz,
            cexpn,cx,cy,cz,cI,cJ,cK,dexpn,dx,dy,dz,M)
        val4 = vrr_iter(cexpn,cx,cy,cz,cI,cJ,cK,dexpn,dx,dy,dz,
            aexpn,ax,ay,az,aI,aJ,aK,bexpn,bx,by,bz,M)
    end
    end# testset
end

function test_hrr()
    @testset "test_hrr" begin
    ax=ay=az=bx=by=bz=cx=cy=cz=dx=dy=dz=0.0
    aexpn=bexpn=cexpn=dexpn=1.0
    aI=aJ=aK=0
    bI,bJ,bK = 1,0,1
    cI=cJ=cK=0
    dI,dJ,dK = 1,0,1

    for (ax,ay,az, aI,aJ,aK, cI,cJ,cK, result) in [
            (0.,0.,0., 0,0,0, 0,0,0, 0.0136667330685),
            (0.,0.,0., 1,0,0, 1,0,0, 0.00821630976139),
            (0.,0.,0., 0,1,0, 0,1,0, 0.00122024402397),
            (0.,0.,0., 0,0,1, 0,0,1, 0.00821630976139),

            (0.,0.,0., 2,0,0, 2,0,0,   0.0039759617781),
            (0.,0.,0., 0,2,0, 0,2,0,   0.000599953311785),
            (0.,0.,0., 0,0,2, 0,0,2,  0.0039759617781),

            (1.,2.,3., 1,0,0, 1,0,0, -1.1851316496333975e-6),
            (1.,2.,3., 0,1,0, 0,1,0,  -4.669991667384835e-6),
            (1.,2.,3., 0,0,1, 0,0,1, -3.474373852654044e-5),

            (1.,2.,3., 2,0,0, 2,0,0, 2.81002247462e-6),
            (1.,2.,3., 0,2,0, 0,2,0, 7.09856891538e-6),
            (1.,2.,3., 0,0,2, 0,0,2, 3.62153023224e-5),

            (0.,0.,0., 1,1,0, 1,1,0, 0.000599953311785),
            (0.,0.,0., 0,1,1, 0,1,1, 0.000599953311785),
            (0.,0.,0., 1,0,1, 1,0,1, 0.0116431617287),

            (3.,2.,1., 1,1,0, 1,1,0, 7.37307761485e-6),
            (3.,2.,1., 0,1,1, 0,1,1, 2.5333243119843164e-7),
            (3.,2.,1., 1,0,1, 1,0,1, 2.452115184675799e-6)
        ]
        #println("hrr($aexpn,$ax,$ay,$az,$aI,$aJ,$aK,$bexpn,$bx,$by,$bz,$bI,$bJ,$bK,")
        #println("    $cexpn,$cx,$cy,$cz,$cI,$cJ,$cK,$dexpn,$dx,$dy,$dz,$dI,$dJ,$dK)")
        val1 = hrr(aexpn,ax,ay,az,aI,aJ,aK,bexpn,bx,by,bz,bI,bJ,bK,
            cexpn,cx,cy,cz,cI,cJ,cK,dexpn,dx,dy,dz,dI,dJ,dK)
        val2 = hrr(cexpn,cx,cy,cz,cI,cJ,cK,dexpn,dx,dy,dz,dI,dJ,dK,
            aexpn,ax,ay,az,aI,aJ,aK,bexpn,bx,by,bz,bI,bJ,bK)
        @test val1 ≈ val2
        @test val1 ≈ result
    end
    end #testset
end


text:
module Trixi2Img

# Include other packages
using EllipsisNotation
using Glob: glob
using HDF5: h5open, attributes, haskey
using Plots: plot, plot!, gr, savefig, contourf!
using TimerOutputs
import GR

# Number of spatial dimensions
"""
    ndims

Number of spatial dimensions (= 2).
"""
const ndim = 2

# Include all source files
include("interpolation.jl")
include("interpolate.jl")
include("io.jl")

# Include top-level conversion method
include("convert.jl")

# export types/functions that define the public API of Trixi2Img
export trixi2img

end # module Trixi2Img


text:
module Utils_Model
using JuMP
using DataFrames

include(dirname(dirname(@__FILE__))*"\\Model\\R_sep_Model.jl")
include(dirname(dirname(@__FILE__))*"\\Model\\C_sep_Model.jl")
using .R_sep_Model
using .C_sep_Model

export export_model,
       batch_model,
       read_model

""" Exmports the JuMP model to SDPA format """
export_model(model,file_name::String = pwd()*"\\Defualt_sep_mod.dat-s") =
    JuMP.write_to_file(model, file_name, format=MOI.FileFormats.FORMAT_SDPA)

""" Reads a JuMP model from SDPA format text document """
function read_model(file_name::String)
    @assert splitext(file_name)[end] == ".dat-s" "Can only load models from .dat-s (SDPA) files!"
    return JuMP.read_from_file(file_name, format=MOI.FileFormats.FORMAT_SDPA)
end

end

max_stars_repo_name: JAndriesJ/sep-rank-bounding

max_issues_repo_name: JAndriesJ/sep-rank-bounding

max_forks_repo_name: JAndriesJ/sep-rank-bounding