text:
using Gamess
using Test

import Gamess: Configuration, substitutions

load_file(name) = Gamess.load(joinpath(dirname(@__FILE__), "data", "$(name).out"))

@testset "Gamess.jl" begin
    @testset "Water CIS" begin
        water = load_file("water")
        atoms = water.molecule.atoms
        @test size(atoms) == (3, 5)
        @test atoms.Atom == ["OXYGEN", "HYDROGEN", "HYDROGEN"]

        ci = water.ci
        states = ci.excited_states
        @test size(states) == (100,3)

        dipoles = ci.length_gauge_dipole
        @test size(dipoles) == (2601,10)
        @test dipoles[1, :z] ≈ 0.958273779741262
        @test dipoles[2599,:z] ≈ 0.044333
        @test dipoles[2599,:A] ≈ -5.4312e-41
        @test dipoles[2599,:B] ≈ 797650.0

        expansions = ci.state_expansions
        @test length(expansions) == 100

        e = first(expansions)
        reference = e.reference
        N = 13
        @test length(reference.α) == length(reference.β) == N
        @test reference == Configuration(vcat(trues(5), falses(N-5)), vcat(trues(5), falses(N-5)))

        @test substitutions(reference, e.configurations[1]) == ([5 => 6], [])
        @test substitutions(reference, e.configurations[2]) == ([5 => 9], [])
        @test substitutions(reference, e.configurations[3]) == ([],[5 => 6])
        @test substitutions(reference, e.configurations[4]) == ([],[5 => 9])
    end

    @testset "Neon CIS" begin
        @testset "Ne I" begin
            neon = load_file("ne-cis")
            atoms = neon.molecule.atoms
            @test size(atoms) == (1, 5)
            @test atoms.Atom == ["NE"]

            ci = neon.ci
            states = ci.excited_states
            @test size(states) == (100,3)

            dipoles = ci.length_gauge_dipole
            @test size(dipoles) == (2603,10)
            @test dipoles[4, :z] ≈ -0.250121

            expansions = ci.state_expansions
            @test length(expansions) == 100

            e = first(expansions)
            reference = e.reference
            N = 121
            @test length(reference.α) == length(reference.β) == N
            @test reference == Configuration(vcat(trues(5), falses(N-5)), vcat(trues(5), falses(N-5)))

            @test substitutions(reference, e.configurations[1]) == ([4 => 10], [])
            @test substitutions(reference, e.configurations[2]) == ([4 => 19], [])
            @test substitutions(reference, e.configurations[13]) == ([],[4 => 10])
            @test substitutions(reference, e.configurations[14]) == ([],[4 => 19])
        end

        @testset "Ne II" begin
            neon = load_file("ne+-guga")
            atoms = neon.molecule.atoms
            @test size(atoms) == (1, 5)
            @test atoms.Atom == ["NE"]

            ci = neon.ci
            states = ci.excited_states
            @test size(states) == (5,2)

            dipoles = ci.length_gauge_dipole
            @test size(dipoles) == (15,10)
            @test dipoles[11, :z] ≈ -0.521326

            dipoles = ci.velocity_gauge_dipole
            @test size(dipoles) == (10,8)
        end
    end
end


text:
struct Serializable{T}
    val::T

    function Serializable(val::T) where {T}
        @assert !isbits(T) "BSONSerializer cannot serialize a bitstype to BSON."
        return new{T}(val)
    end
end

ext: jl


text:
using Combinatorics, Statistics, DataFrames, CSV

data = CSV.read("../data/fertilizer.csv")
control = data.Control
fertilizer = data.FertilizerX

subGroups = collect(combinations([control;fertilizer],10))

meanFert = mean(fertilizer)
pVal = sum([mean(i) >= meanFert for i in subGroups])/length(subGroups)
println("p-value = ", pVal)

max_stars_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_issues_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_forks_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook


text:
using Knet: cuda11

function cuda11src(f, j=f, ex="$f(xi,yi)"; BLK=256, THR=256)
"""
__global__ void _$(f)_32_11(int n, float *x, float *y, float *z) {
  int i = threadIdx.x + blockIdx.x * blockDim.x;
  while (i < n) {
    float xi=x[i];
    float yi=y[i];
    z[i] = $ex;
    i += blockDim.x * gridDim.x;
  }
}
extern "C" {
  void $(f)_32_11(int n, float *x, float *y, float *z) {
    _$(f)_32_11<<<$BLK,$THR>>>(n,x,y,z);
  }
}
__global__ void _$(f)_64_11(int n, double *x, double *y, double *z) {
  int i = threadIdx.x + blockIdx.x * blockDim.x;
  while (i < n) {
    double xi=x[i];
    double yi=y[i];
    z[i] = $ex;
    i += blockDim.x * gridDim.x;
  }
}
extern "C" {
  void $(f)_64_11(int n, double *x, double *y, double *z) {
    _$(f)_64_11<<<$BLK,$THR>>>(n,x,y,z);
  }
}
"""
end

for a in cuda11
    isa(a,Tuple) || (a=(a,))
    print(cuda11src(a...))
end


text:
using Documenter
using EnergyExpressions
using AtomicLevels

makedocs(
    modules = [EnergyExpressions],
    sitename = "EnergyExpressions",
    pages = [
        "Home" => "index.md",
        "Theory" => [
            "Notation" => "notation.md",
            "Energy Expressions" => "energy_expressions.md",
            "Calculus of Variations" => "calculus_of_variations.md"
        ],
        "Implementation" => [
            "Conjugate orbitals" => "conjugate_orbitals.md",
            "Slater determinants" => "slater_determinants.md",
            "N-body operators" => "nbody_operators.md",
            "N-body matrix elements" => "nbody_matrix_elements.md",
            "Common N-body operators" => "common_operators.md",
            "N-body equations" => "equations.md",
            "Variation" => "variations.md",
            "System of equations" => "system_of_equations.md"
        ]
    ],
    format = Documenter.HTML(assets = ["assets/latex.js"]),
    doctest = false
)

deploydocs(repo = "github.com/JuliaAtoms/EnergyExpressions.jl.git")


text:
"""
    RandomDotOrgAPI

Use various functions from https://random.org (q.v.).

# Currently available
- get_usage(): obtain current bit quota for your IP.
- check_usage(): check if quota is non-zero.
- get_result(): get a stored signed result from RDO (stored >= 24 hours)
- verify_signature(): obtain signature verification for a previously signed result from RDO
- generate_integers(): obtain integers.
- generate_integer_sequences(): obtain randomized sequences of integers 1..N
- generate_strings(): obtain random strings of characters (upper/lower case letters, digits)
- generate_gaussians(): obtain random Gaussian numbers
- generate_decimal_fractions(): obtain random numbers on the interval (0,1)
- generate_uuids(): obtain random bytes in various formats
- generate_blobs(): generate bit blobs of sizes in multiples of 8 in base64 or hex
- pull_data(): extract data from RDO response (Dictionary)

Github repository at: https://github.com/edwcarney/RandomDotOrgAPI

"""
module RandomDotOrgAPI

using HTTP, JSON
# using OrderedCollections

const url = "https://api.random.org/json-rpc/2/invoke"
<<<<<<< HEAD
include("RandomDotOrg.key")
=======
const myapikey = "00000000-0000-0000-0000-000000000000"
>>>>>>> 02f29c07e06d8f6448742aada08e69066b19e03f
const myapikeyzeros = "00000000-0000-0000-0000-000000000000"

export  get_usage, check_usage, get_result, verify_signature, generate_integers, generate_integer_sequences, generate_strings, generate_gaussians,
        generate_decimal_fractions, generate_uuids, generate_blobs, pull_data

"""
    Get the current bit quota from Random.org
"""
function get_usage(; apiType = "basic")
    get_usage1 = Dict(
        "jsonrpc" => "2.0",
        "method" => "getUsage",
        "params" => Dict(
            "apiKey" => myapikeyzeros
            ),
        "id" => 22407
    )

    if apiType == "signed"
         get_usage1["params"]["apiKey"] = myapikey
    end

    r = HTTP.request("POST", url,
                ["Content-Type" => "application/json"],
                JSON.json(get_usage1))
    out = JSON.parse(String(r.body))
    return out
end;

"""

    check_usage(minimum = 500::Number)

    Test for sufficient quota to insure response. This should be set to match
    user's needs.
"""
function check_usage(minimum = 500::Number)

    return (get_usage()["result"]["bitsLeft"] >= minimum);
end;

"""
    get_result(serialnumber::Int)

    get_result(): get a stored signed result from RDO (stored >= 24 hours)

"""
function get_result(serialNumber::Int)
    get_result = Dict(
        "jsonrpc" => "2.0",
        "method" => "getResult",
        "params" => Dict(
            "apiKey" => myapikey,
            "serialNumber" => serialNumber
        ),
        "id" => 22407
    )

    r = HTTP.request("POST", url,
                ["Content-Type" => "application/json"],
                JSON.json(get_result))
    out = JSON.parse(String(r.body))
    return out
end;

"""
    verify_signature(result::Dict)

Verify a signature from a previous signed request.

# Arguments
- `result` : Dict with signed result from RDO for previous request

# Example
```
r1 = get_result(672)
verify_signature(r1)
"result" => Dict{String,Any} with 1 entry:
    "authenticity" => true
```
"""
function verify_signature(result::Dict)
    verifSignature = Dict(
        "jsonrpc" => "2.0",
        "method" => "verifySignature",
        "params" => Dict(
            "random" => result["result"]["random"],
            "signature" => result["result"]["signature"]
        ),
        "id" => 22407
    )

    r = HTTP.request("POST", url,
                ["Content-Type" => "application/json"],
                JSON.json(verifSignature))
    out = JSON.parse(String(r.body))
    return out
end

"""
    generate_integers(n = 100::Number, min = 1, max = 20, base = 10, parse=true, check = true, replace = true, apitype = "basic")

Get `n` random integers on the interval `[min, max]` as strings
in one of 4 `base` values--[binary (2), octal (8), decimal (10), or hexadecimal (16)]
All numbers, except base 10, are returned as strings.

# Arguments
- `max`,`min` : [-1e9, 1e9]
- `base::Integer`: retrieved Integer format [2, 8, 10, 16]
- `numeric`: return numbers instead of strings (base = 10, only)
- `check::Bool`: perform a call to `checkQuota` before making request
- `replace::Bool`: use sampling with replacement
- `apiType::String`: "basic" or "signed"; key of zeros if "basic"

# Examples
```
generate_integers(5, max=50, base=16)
"result" => Dict{String,Any} with 5 entries:
    "bitsLeft"      => 3361463
        "random"        => Dict{String,Any} with 2 entries:
            "data"           => Any["1d", "0f", "0b", "11", "1e"]
            "completionTime" => "2020-08-01 22:54:49Z"
    "advisoryDelay" => 3150
    "bitsUsed"      => 28
    "requestsLeft"  => 773473

julia> generate_integers(5, max=4096, base=2, min=2048)
"result" => Dict{String,Any} with 5 entries:
    "bitsLeft"      => 3363817
        "random"        => Dict{String,Any} with 2 entries:
            "data"           => Any["110110110111", "100111010010", "111001100010", "111100100011", "110100101111"]
            "completionTime" => "2020-08-01 22:50:06Z"
    "advisoryDelay" => 3670
    "bitsUsed"      => 55
    "requestsLeft"  => 773569

julia> generate_integers(5, max=200)
"result" => Dict{String,Any} with 5 entries:
    "bitsLeft"      => 3362373
        "random"        => Dict{String,Any} with 2 entries:
            "data"           => Any[168, 168, 4, 101, 175]
            "completionTime" => "2020-08-01 22:52:27Z"
    "advisoryDelay" => 2970
    "bitsUsed"      => 38
    "requestsLeft"  => 773510
```
"""
function generate_integers(n = 100::Number; min = 1, max = 20, base = 10, numeric = true, check = true, replace=true, apiType="basic")
    (n < 1 || n > 10000) && return "Requests must be between 1 and 10,000 numbers"

    (min < -1f+09 || max > 1f+09 || min > max) && return "Range must be between -1e9 and 1e9"

    (!(base in [2, 8, 10, 16])) && return "Base has to be one of 2, 8, 10 or 16"

    (check && !check_usage()) && return "random.org suggests to wait until tomorrow"

#     urlbase = "https://www.random.org/integers/"
#     urltxt = @sprintf("%s?num=%d&min=%d&max=%d&col=%d&base=%d&format=plain&rnd=new",
#                     urlbase, n, Int(min), Int(max), col, base)
    genIntegers = Dict{String, Any}(
                "jsonrpc" => "2.0",
                "method" => "generate_integers",
                "params" => Dict{String, Any}(
                        "n" => n,
                        "min" => min,
                        "max" => max,
                        "base" => base,
                        "replacement" => replace
                ),
                "id" => 22407
    )

    apiType=="basic" ? push!(genIntegers, "method" => "generateIntegers") : push!(genIntegers, "method" => "generateSignedIntegers")
    apiType=="basic" ? push!(genIntegers["params"], "apiKey" =>   myapikeyzeros) : push!(genIntegers["params"], "apiKey" =>  myapikey)

    r = HTTP.request("POST", url,
                ["Content-Type" => "application/json"],
                JSON.json(genIntegers))
    out = JSON.parse(String(r.body))

    return out
end

"""
    generate_integer_sequences(n = 1; min = 1::Number, max = 20::Number, col = 1, check = true, replace=true, apiType="basic")

Get a randomized interval `[min, max]`. Returns (max - min + 1) randomized integers
All numbers are returned as strings (as Random.org sends them).

# Arguments
- `length` : size of the returned sequence; may be a tuple
- `min` : no less than 1; may be a tuple, as [1, 1]
- `max` : must be [-1e9, 1e9]; may be a tuple, as [69, 26]
- `base::Integer`: retrieved Integer format [2, 8, 10, 16]
- `check::Bool`: perform a call to `checkQuota` before making request
- `replace::Bool`: use sampling with replacement
- `apiType::String`: "basic" or "signed"; key of zeros, if "basic"

Returns a Dict with results in "result", data are nested in ["random"]["data"]. Multiple sequences
must be teased out with deeper slicing. For example, the first sequence in a two sequence return
will require the following: ["result"]["random"]["data"][1]

Additional information, such as the nuber of bits left in your quota will be provided. See this page for
further details on this function: https://api.random.org/json-rpc/2/basic

# Example of a simple sequence of random order of 1:10 with replace = false
```
julia> generate_integer_sequences(1, max=10)
"result"  => Dict{String,Any} with 5 entries:
    "bitsLeft"      => 3830058
        "random"        => Dict{String,Any} with 2 entries:
            "data"=>Any[Any[4, 6, 3, 2, 8, 1, 7, 10, 9, 5]],
            "completionTime"=>"2020-07-31 04:55:14Z"
    "advisoryDelay" => 2850 # recommended wait time for next request in ms
    "bitsUsed"      => 33
    "requestsLeft"  => 792944

 # Example of a multiple sequence from 1 to 10 and 1 to 5, with 10 samples from each with replace = true
 generate_integer_sequences(2, min = [1, 1],  max=[10, 5], replace=true)
"result"   => Dict{String,Any} with 5 entries:
    "bitsLeft"=>3825370,
        "random"  => Dict{String,Any} with 2 entries:
            "data"=>Any[Any[7, 2, 6, 1, 5, 3, 5, 2, 1, 9], Any[5, 3, 4, 2, 4, 3, 5, 1, 1, 4]],
            "completionTime"=>"2020-07-31 05:01:10Z"
    "advisoryDelay" => 2480
    "bitsUsed"      => 56
    "requestsLeft"  => 792747
```
"""
function generate_integer_sequences(n = 10::Number, length = 10::Number; min = 1::Number, max = 20::Number, base = 10, check = true, replace = true, apiType = "basic")
    # (min < -1f+09 || max > 1f+09 || min > max) && return "Range must be between -1e9 and 1e9"

    (check && !check_usage()) && return "random.org suggests to wait until tomorrow"

    genIntegerSequences = Dict{String, Any}(
        "jsonrpc" => "2.0",
        "method" => "genIntegerSequences",
        "params" => Dict{String, Any}(
                "n" => n,
                "length" => length,
                "min" => min,
                "max" => max,
                "base" => base,
                "replacement" => replace
        ),
        "id" => 22407
    )

    apiType=="basic" ? push!(genIntegerSequences, "method" => "generate_integer_sequences") : push!(genIntegerSequences, "method" => "generateSignedIntegerSequences")
    apiType=="basic" ? push!(genIntegerSequences["params"], "apiKey" =>   myapikeyzeros) : push!(genIntegerSequences["params"], "apiKey" =>  myapikey)

    r = HTTP.request("POST", url,
                ["Content-Type" => "application/json"],
                JSON.json(genIntegerSequences))
    out = JSON.parse(String(r.body))

    return out;
end

"""

    generate_strings(n=10::Number, length=5, characters="abcdefghijklmnopqrstuvwxyz"; check=true)

Get `n` random strings of length `len`

# Arguments
- `n::number` : number of strings
- `length::Number` : length of strings [1, 32]
- `characters::String`: strings formed from this set of characters
- `check::Bool`: perform a call to `checkQuota` before making request
- `replace::Bool`: use sampling with replacement (strings might be duplicated)
- `apiType::String`: "basic" or "signed"; "basic" sends zeros as  myapikey

# Examples
```
# Generate 2 strings of 10 characters each from the lower case alphabet with replacement

generate_strings(2, 10, "abcdefghijklmnopqrstuvwxyz", apiType="basic")
"result"  => Dict{String,Any} with 5 entries:
    "bitsLeft"  => 3812394
        "random"    => Dict{String,Any} with 2 entries:
            "data"           => Any["crarvpazcz", "iagtbmseus"]
            "completionTime" => "2020-07-31 05:27:21Z"
    "advisoryDelay" => 1560
    "bitsUsed"      => 56
    "requestsLeft"  => 792218

# Generate a single string of 15 characters from the upper and lower case alphabets + numerals
# without replacement

generate_strings(1, 15, "abcdefghijklmnopqrstuvwxyzABCDEFGHIJKLMNOPQRSTUVWXYZ0123456789"; replace=false, apiType="basic")
"result"  => Dict{String,Any} with 5 entries:
    "bitsLeft"      => 3356382
    "random"        => Dict{String,Any} with 2 entries:
        "data"           => Any["BEuIiiYud34tGa6"]
        "completionTime" => "2020-08-01 23:00:07Z"
    "advisoryDelay" => 2740
    "bitsUsed"      => 89
    "requestsLeft"  => 773264
```
"""
function generate_strings(n = 10::Number, length = 5::Number, characters = "abcdefghijklmnopqrstuvwxyz"::String; check = true, replace = true, apiType="basic")
    (n < 1 || n > 10000) && return "1 to 10,000 requests only"

    (length < 1 || length > 32) && return "Length must be between 1 and 32"

    (check && !check_usage()) && return "random.org suggests to wait until tomorrow"

    genStrings = OrderedDict{String, Any}(
        "jsonrpc" => "2.0",
        "params" => Dict{String, Any}(
                "n" => n,
                "length" => length,
                "characters" => characters,
                "replacement" => replace
        ),
        "id" => 22407
    )

    apiType=="basic" ? push!(genStrings, "method" => "generateStrings") : push!(genStrings, "method" => "generateSignedStrings")
    apiType=="basic" ? push!(genStrings["params"], "apiKey" =>   myapikeyzeros) : push!(genStrings["params"], "apiKey" =>  myapikey)

    @show JSON.json(genStrings)

    r = HTTP.request("POST", url,
                ["Content-Type" => "application/json"],
                JSON.json(genStrings))
    out = JSON.parse(String(r.body))

    return out;
end

"""

    generate_gaussians(n = 10::Number, mean = 0.0, stdev = 1.0, digits = 10; check = true, replace = true, apiType = "basic")

Get n numbers from a Gaussian distribution with `mean` and `stdev`.
Returns strings in `dec` decimal places.
Scientific notation only for now.

# Arguments
- `mean`, `stdev` : between [-1e6, 1e6]
- `significantDigits` : decimal places [2,14]
- `check::Bool`: perform a call to `checkQuota` before making request

# Example requesting 500 Gaussians with mean = 5.0, stdev = 2.5, 10 significant digits
```
julia> generate_gaussians(500, 5, 2.5, 10, apiType="basic")
"result" => Dict{String,Any} with 5 entries:
    "bitsLeft"      => 3784434
    "random"        => Dict{String,Any} with 2 entries:
                        "data"           => Any[3.45881, 0.466094, 5.20018, 6.11947, 0.843332, 5.51932, 4.23481, 0.0137821, -1.96319, 3.56091
                                                . . . 6.06737, 7.00533, 12.6734, 1.24946, 4.35043, 5.5967, 5.87513, 4.49188, -0.364909, 4.10236]
                        "completionTime" => "2020-07-31 06:01:50Z"
    "advisoryDelay" => 1990
    "bitsUsed"      => 16610
    "requestsLeft"  => 791750
```
"""
function generate_gaussians(n = 10::Number, mean = 0.0, stdev = 1.0, digits = 10; check = true, apiType = "basic")
    (n < 1 || n > 10000) && return "Requests must be between 1 and 10,000 numbers"

    (mean < -1f+06 || mean > 1f+06) && return "Mean must be between -1e6 and 1e6"

    (stdev < -1f+06 || stdev > 1f+06) && return "Std dev must be between -1e6 and 1e6"

    (digits < 2 || digits > 14) && return "Decimal places must be between 2 and 14"

    (check && !check_usage()) && return "random.org suggests to wait until tomorrow"

    genGaussians = Dict{String, Any}(
                "jsonrpc" => "2.0",
                "params" => Dict{String, Any}(
                        "n" => n,
                        "mean" => mean,
                        "standardDeviation" => stdev,
                        "significantDigits" => digits
                ),
                "id" => 22407
    )

    apiType=="basic" ? push!(genGaussians, "method" => "generateGaussians") : push!(genGaussians, "method" => "generateSignedGaussians")
    apiType=="basic" ? push!(genGaussians["params"], "apiKey" =>   myapikeyzeros) : push!(genGaussians["params"], "apiKey" =>  myapikey)

    r = HTTP.request("POST", url,
                ["Content-Type" => "application/json"],
                JSON.json(genGaussians))
    out = JSON.parse(String(r.body))

    return out
end

"""
    generate_decimal_fractions(n = 10::Number; digits = 5, check = true, replace = true)

Get n decimal fractions on the interval (0,1).
Returns strings in `digits` decimal places.

# Arguments
- `n`: number of fractions to generate
- `digits` : decimal places [1, 14]
- `check::Bool`: perform a call to `checkQuota` before making request
- `replace::Bool`: sample with replacement (duplicates possible)

# Example to obtain 5 decimal fractions on the interval [0, 1)
```
generate_decimal_fractions(5, 10, replace=false)
"result" => Dict{String,Any} with 5 entries:
    "bitsLeft"      => 3722130
    "random"        => Dict{String,Any} with 2 entries:
            "data"           => Any[0.532242, 0.4899, 0.430238, 0.709279, 0.150412]
            "completionTime" => "2020-07-31 07:23:00Z"
    "advisoryDelay" => 2640
    "bitsUsed"      => 166
    "requestsLeft"  => 790030
```
"""
function generate_decimal_fractions(n = 10, digits = 10; check=true, replace = true, apiType = "basic")
    (n < 1 || n > 10000) && return "Requests must be between 1 and 10,000 numbers"

    (digits < 1 || digits > 14) && return "Decimal places must be between 1 and 14"

    (check && !check_usage()) && return "random.org suggests to wait until tomorrow"

    genDecimalFractions = Dict{String, Any}(
        "jsonrpc" => "2.0",
        "params" => Dict{String, Any}(
                "n" => n,
                "decimalPlaces" => digits,
                "replacement" => replace
        ),
        "id" => 22407
    )

    apiType=="basic" ? push!(genDecimalFractions, "method" => "generateDecimalFractions") : push!(genDecimalFractions, "method" => "generateSignedDecimalFractions")
    apiType=="basic" ? push!(genDecimalFractions["params"], "apiKey" =>   myapikeyzeros) : push!(genDecimalFractions["params"], "apiKey" =>  myapikey)

    r = HTTP.request("POST", url,
            ["Content-Type" => "application/json"],
            JSON.json(genDecimalFractions))
    out = JSON.parse(String(r.body))

    return out
end

"""
    generate_uuids(n = 10::Number; check=true, apiType = "basic")

Get n UUIDs.
Returns n UUIDs as strings.

# Arguments
- `n`: number of UUIDs to generate
- `check::Bool`: perform a call to `checkQuota` before making request

# Example to obtain 1 UUIDs
```
generate_uuids(1)
"result" => Dict{String,Any} with 5 entries:
    "bitsLeft"      => 3396647
        "random"        => Dict{String,Any} with 2 entries:
            "data"           => Any["12598e3a-c96a-4b16-bea4-b7d404717080"]
            "completionTime" => "2020-08-01 21:48:39Z"
    "advisoryDelay" => 3070
    "bitsUsed"      => 122
    "requestsLeft"  => 774929
```
"""
function generate_uuids(n = 10::Number; check=true, apiType = "basic")
    (n < 1 || n > 1000) && return "Requests must be between 1 and 1000 numbers"

    (check && !check_usage()) && return "random.org suggests to wait until tomorrow"

    genUUIDs = Dict{String, Any}(
        "jsonrpc" => "2.0",
        "params" => Dict{String, Any}(
                "n" => n,
        ),
        "id" => 22407
)

    apiType=="basic" ? push!(genUUIDs, "method" => "generateUUIDs") : push!(genUUIDs, "method" => "generateSignedUUIDs")
    apiType=="basic" ? push!(genUUIDs["params"], "apiKey" =>   myapikeyzeros) : push!(genUUIDs["params"], "apiKey" =>  myapikey)

    r = HTTP.request("POST", url,
            ["Content-Type" => "application/json"],
            JSON.json(genUUIDs))
    out = JSON.parse(String(r.body))

    return out
end

"""
    generate_blobs(n = 10::Number, size = 80; format = 'base64', check = true, apiType = "basic")

Get n Blobs.
Returns n Blobs as strings in format requested.

# Arguments
- `n`: number of blobs   to generate
- `size` : size in bits (divisible by 8)
- `format` : "base64" [default] or "hex"
- `check::Bool`: perform a call to `checkQuota` before making request

# Example to obtain 1 UUIDs
```
generate_blobs(1)
Dict{String,Any} with 5 entries:
    "bitsLeft"      => 3335873
    "random" => Dict{String,Any} with 2 entries:
        "data"           => Any["342ae094e29e9814ef21"]
        "completionTime" => "2020-08-03 21:14:38Z"
    "advisoryDelay" => 3650
    "bitsUsed"      => 32
    "requestsLeft"  => 774042
```
"""
function generate_blobs(n = 10::Number, size = 100; format = "base64", check=true, apiType = "basic")
    (n < 1 || n > 100) && return "Requests must be between 1 and 100 numbers"

    (size < 1 || size > 2^20) && return "Size must be between 1 and 1048576"

    (mod(size,8) > 0) && return "Size must be evenly divisible by 8"

    (format != "base64" && format != "hex") && return "Format must be 'hex' or 'base64'"

    (check && !check_usage()) && return "random.org suggests to wait until tomorrow"

    genBlobs = Dict{String, Any}(
        "jsonrpc" => "2.0",
        "params" => Dict{String, Any}(
                "n" => n,
                "size" => size,
                "format" => format
        ),
        "id" => 22407
    )

    apiType=="basic" ? push!(genBlobs, "method" => "generateBlobs") : push!(genBlobs, "method" => "generateSignedBlobs")
    apiType=="basic" ? push!(genBlobs["params"], "apiKey" =>   myapikeyzeros) : push!(genBlobs["params"], "apiKey" =>  myapikey)

    r = HTTP.request("POST", url,
            ["Content-Type" => "application/json"],
            JSON.json(genBlobs))
    out = JSON.parse(String(r.body))

    return out
end

"""
    pull_data(RDOresp)

Extract data from RDO response (Dictionary)
Returns numerical data in form of Vector{Real}
Otherwise, returns vector of data

# Arguments
- `RDOresp`: response from RDO in form of dictionary
"""
function pull_data(RDOresp)
    dataval = RDOresp["result"]["random"]["data"]
    if isa(dataval[1], Number)
        return convert(Vector{Real},dataval)
    elseif isa(dataval[1], Vector{Any})
        if isa(dataval[1][1], Number)
            return [[convert(Real, x) for x in y] for y in dataval]
        else
            return dataval
        end
    else
        return dataval
    end
end

end # module


text:
include("plot_domain.jl")
include("plot_field.jl")
include("plot_moments.jl")

# Plot the result across angular frequency for a specific position (x)
@recipe function plot(simres::SimulationResult;
        x = simres.x,
        x_indices = [findmin([norm(z - y) for z in simres.x])[2] for y in x],
        ω_indices = Colon(), field_apply = real)

    for x_ind in x_indices

        fs = field_apply.(field(simres)[x_ind, ω_indices])
        xguide = ((typeof(simres) <: FrequencySimulationResult) ? "ω" : "t")

        @series begin
            label --> "$field_apply x=$(simres.x[x_ind])"
            xguide --> xguide
            (getfield(simres, 3)[ω_indices], fs)
        end
    end
end

"Plot just the particles"
@recipe function plot(sim::FrequencySimulation; bounds = :none)

    # println("Plotting a simulation on its own")

    @series begin
        if bounds != :none
            # bounds = bounding_box(sim.particles)
            xlims --> (bottomleft(bounds)[1], topright(bounds)[1])
            ylims --> (bottomleft(bounds)[2], topright(bounds)[2])
        end

        sim.particles
    end

end

"Plot the field for a particular wavenumber"
@recipe function plot(sim::FrequencySimulation{T,3}, ω::T;
        res=10, xres=res, yres=res,
        y = :auto,
        field_apply=real,
        region_shape = :auto,
        bounds = :auto,
        exclude_region = EmptyShape{T,3}(),
        drawparticles=false) where {T}

    # If user wants us to, generate bounding rectangle around particles
    region_shape = (region_shape != :auto) ? region_shape :
        if isempty(sim.particles)
            if bounds == :auto
                @warn "What region to plot? For example, use keyword bounds = Box([[-1.0,-1.0],[1.0,1.0]])"
                Box([[-one(T),-one(T)],[one(T),one(T)]])
            else bounds
            end
        else
            region_shape = bounding_box(sim.particles)
        end

    bounds = bounding_box(region_shape)

    # If user has not set xlims and ylims, set them to the rectangle
    cs = corners(bounds)
    xlims --> (cs[1][1], cs[end][1])
    ylims --> (cs[1][end], cs[end][end])

    if y == :auto
        y = cs[1][2]
    end

    # Incase the user did set the xlims and ylims, generate a new bounding box with them
    # p_xlims = plotattributes[:xlims]
    # p_ylims = plotattributes[:ylims]
    # bounds = Box([[T(p_xlims[1]),T(p_ylims[1])], [T(p_xlims[2]),T(p_ylims[2])]])
    #
    # region_shape = (bounds ⊆ region_shape) ? bounds : region_shape

    field_sim = run(sim, region_shape, [ω]; y=y, xres=xres, zres=yres, exclude_region=exclude_region)
    xy_mat = reshape(field_sim.x, (xres+1, yres+1))
    x_pixels = [x[1] for x in xy_mat[:,1]]
    y_pixels = [x[end] for x in xy_mat[1,:]]

    @series begin

        # Turn the responses (a big long vector) into a matrix, so that the heatmap will understand us
        response_mat = transpose(reshape(field(field_sim), (xres+1, yres+1)))
        seriestype --> :contour
        fill --> true
        grid --> false
        aspect_ratio := 1.0
        seriescolor --> :balance
        title --> "Field at ω=$ω"

        (x_pixels, y_pixels, field_apply.(response_mat))
    end

    if drawparticles
        @series begin
            sim.particles
        end
    end

end

"Plot the field for a particular wavenumber"
@recipe function plot(sim::FrequencySimulation{T,2}, ω::T;
        res=10, xres=res, yres=res,
        field_apply=real,
        region_shape = :auto,
        bounds = :auto,
        exclude_region = EmptyShape{T,2}(),
        drawparticles=false) where {T}

    # If user wants us to, generate bounding rectangle around particles
    region_shape = (region_shape != :auto) ? region_shape :
        if isempty(sim.particles)
            if bounds == :auto
                @warn "What region to plot? For example, use keyword bounds = Box([[-1.0,-1.0],[1.0,1.0]])"
                Box([[-one(T),-one(T)],[one(T),one(T)]])
            else bounds
            end
        else
            region_shape = bounding_box(sim.particles)
        end

    bounds = bounding_box(region_shape)
    # If user has not set xlims and ylims, set them to the rectangle
    xlims --> (bottomleft(bounds)[1], topright(bounds)[1])
    ylims --> (bottomleft(bounds)[2], topright(bounds)[2])

    # Incase the user did set the xlims and ylims, generate a new bounding
    # rectangle with them
    p_xlims = plotattributes[:xlims]
    p_ylims = plotattributes[:ylims]
    bounds = Box([[T(p_xlims[1]),T(p_ylims[1])], [T(p_xlims[2]),T(p_ylims[2])]])

    region_shape = (bounds ⊆ region_shape) ? bounds : region_shape

    field_sim = run(sim, region_shape, [ω]; xres=xres, yres=yres, zres=yres, exclude_region=exclude_region)
    xy_mat = reshape(field_sim.x, (xres+1, yres+1))
    x_pixels = [x[1] for x in xy_mat[:,1]]
    y_pixels = [x[2] for x in xy_mat[1,:]]

    @series begin

        # Turn the responses (a big long vector) into a matrix, so that the heatmap will understand us
        response_mat = transpose(reshape(field(field_sim), (xres+1, yres+1)))
        seriestype --> :contour
        fill --> true
        grid --> false
        aspect_ratio := 1.0
        seriescolor --> :balance
        title --> "Field at ω=$ω"

        (x_pixels, y_pixels, field_apply.(response_mat))
    end

    if drawparticles
        @series begin
            sim.particles
        end
    end

end

"Plot the source field for a particular wavenumber"
@recipe function plot(s::RegularSource, ω)
    (FrequencySimulation(s), ω)
end


text:
module ICARTTExceptions

using JLLUtils.JLLExceptions; # provides @msgexc

"""
    ICARTTParsingException(msg)

An exception class raised if any issues arise while parsing an ICARTT-formatted
file.
"""
@msgexc ICARTTParsingException

"""
    ICARTTNotImplementedException(msg)

An exception raised if a particular case is recognized, but not implemented, by
the parser.
"""
@msgexc ICARTTNotImplementedException

"""
    ICARTTUnitsException(msg)

An exception raised if there is a problem interpreting a unit in an ICARTT file.
"""
@msgexc ICARTTUnitsException

end

max_stars_repo_name: joshua-laughner/Julia-ICARTT

max_issues_repo_name: joshua-laughner/Julia-ICARTT

max_forks_repo_name: joshua-laughner/Julia-ICARTT


text:
using Term, Documenter

DocMeta.setdocmeta!(Term, :DocTestSetup, :(using Term); recursive = true)

makedocs(;
    modules = [Term],
    authors = "FedeClaudi <federicoclaudi@protonmail.com> and contributors",
    repo = "https://github.com/FedeClaudi/Term.jl/blob/{commit}{path}#{line}",
    sitename = "Term.jl",
    format = Documenter.HTML(;
        prettyurls = get(ENV, "CI", "false") == "true",
        canonical = "https://FedeClaudi.github.io/Term.jl",
        assets = String[],
        collapselevel = 1,
    ),
    strict=false,
    pages = [
        "Home" => "index.md",
        "Basics" => Any[
            "basics/basics.md",
            "basics/styled_text.md",
            "basics/colors.md",
            "basics/renderables.md",
            "basics/content_layout.md",
            "basics/tprint.md"
        ],
        "Renderables" => Any[
            "ren/intro.md",
            "ren/text.md",
            "ren/panel.md",
            "ren/tbox.md",
            "ren/layout_rens.md",
            "ren/tree.md",
        ],
        "Advanced" => Any[
            "adv/adv.md",
            "adv/progressbars.md",
            "adv/logging.md",
            "adv/errors_tracebacks.md",
            "adv/introspection.md",
        ],
        "API" => Any[
            "api/api_term.md",
            "api/api_box.md",
            "api/api_color.md",
            "api/api_console.md",
            "api/api_errors.md",
            "api/api_introspection.md",
            "api/api_layout.md",
            "api/api_logging.md",
            "api/api_markup.md",
            "api/api_measure.md",
            "api/api_panel.md",
            "api/api_renderables.md",
            "api/api_segment.md",
            "api/api_style.md",
            "api/api_tprint.md",
            "api/api_tree.md",
        ],
    ],
)

deploydocs(; repo = "github.com/FedeClaudi/Term.jl",)


text:
# Plot ROC curve
function plot_ROC_curve(y_cond::Vector{Int}, y_prob::Vector{Float64})

    # Fixed counts
    population = length(y_cond)
    cond_pos = sum(y_cond .== 1)
    cond_neg = sum(y_cond .== 0)
    accuracy = round(sum(y_cond .== round(Int, y_prob)) / population, 2)

    # Number of plot points
    points = 1001
    max = points - 1

    # Calculate TPR and FPR for varied discrimination
    TPR, FPR = zeros(points), zeros(points)
    for d = 0:1:max
        TP, FP = 0, 0
        for i=1:population
            if y_prob[i] >= (max - d) / max
                y_cond[i] == 1 ? TP += 1 : FP += 1
            end
        end
        TPR[d + 1] = TP / cond_pos
        FPR[d + 1] = FP / cond_neg
    end

    # Accuracy
    TP, FP = 0, 0
    for i = 1:population
        if y_prob[i] >= 0.5
            y_cond[i] == 1 ? TP += 1 : FP += 1
        end
    end
    TPR1 = TP / cond_pos
    FPR1 = FP / cond_neg
    df1 = DataFrame(x = FPR1, y = TPR1, label = "Accuracy = $accuracy")

    # AUC
    AUC = 0
    for i=2:length(TPR)
        AUC += (TPR[i] + TPR[i - 1]) * 0.5 * (FPR[i] - FPR[i - 1])
    end
    AUC = round(AUC, 2)
    df2 = DataFrame(x = [0.65], y = [0.3], label = "AUC = $AUC")

    Sigma = Char(931)

    # Plot ROC curve
    plot(
         layer(df1, x="x", y="y", label="label", Geom.point, Geom.label,
               Stat.xticks(ticks=[FPR1]), Stat.yticks(ticks=[TPR1]),
               Theme(default_color=colorant"red", default_point_size=3pt)),
         layer(x=FPR, y=TPR, Geom.step),
         layer(x=[0.0; 1.0], y=[0.0; 1.0], Geom.line,
               Stat.xticks(ticks=[0.0, 0.5, 1.0]), Stat.yticks(ticks=[0.0, 0.5, 1.0]),
               Theme(default_color=colorant"grey", line_style=Gadfly.get_stroke_vector(:dash))),
         layer(df2, x="x", y="y", label="label", Geom.label),
         Scale.x_continuous(minvalue=0., maxvalue=1., labels=x -> @sprintf("%0.2f", x)),
         Scale.y_continuous(minvalue=0., maxvalue=1., labels=y -> @sprintf("%0.2f", y)),
         Guide.ylabel("True Positive Rate \n
($Sigma True Positives / $Sigma Condition Positive)"),
         Guide.xlabel("False Positive Rate \n
($Sigma False Positives / $Sigma Condition Negative)")
         )
end


text:
using AlgebraicDynamics
using AlgebraicDynamics.UWDDynam
using AlgebraicDynamics.DWDDynam
using AlgebraicDynamics.CPortGraphDynam
using AlgebraicDynamics.CPortGraphDynam: barbell
using Catlab.WiringDiagrams
using Catlab.Programs # for @relation macro
using Catlab.Graphics # for to_graphviz
using OrdinaryDiffEq
using DataFrames
using Plots
using StatsPlots

δt = 0.1
tmax = 40
tspan = (0.0,tmax)
t = 0:δt:tmax;

u0 = [990.0,10.0];

β, γ = 0.05*10/1000, 0.25; # in other tutorials, βc/N and γ

dots(u, x, p, t) = [-β*u[1]*x[1]]
doti(u, x, p, t) = [β*x[1]*u[1] - γ*u[1]];

susceptible_cm = ContinuousMachine{Float64}(1,1,1, dots, u -> u)
infected_cm    = ContinuousMachine{Float64}(1,1,1, doti, u -> u);

directed_pattern = WiringDiagram([], [])
susceptible_box = add_box!(directed_pattern, Box(:susceptible_cm, [:pop], [:pop]))
infected_box = add_box!(directed_pattern, Box(:infected_cm, [:pop], [:pop]));

add_wires!(directed_pattern, Pair[
    (susceptible_box, 1) => (infected_box, 1),
    (infected_box, 1) => (susceptible_box, 1)
])

to_graphviz(directed_pattern)

directed_system = oapply(directed_pattern, [susceptible_cm, infected_cm]);

directed_prob = ODEProblem(directed_system, u0, tspan)
directed_sol = solve(directed_prob, FRK65(0));

plot(directed_sol)

cpg_pattern = barbell(1)

cpg_system = oapply(cpg_pattern, [susceptible_cm, infected_cm]);

cpg_prob = ODEProblem(cpg_system, u0, tspan)
cpg_sol = solve(cpg_prob, FRK65(0));

plot(cpg_sol)

dotsi(u,p,t) = [-β*u[1]*u[2],β*u[1]*u[2]]
doti(u,p,t) = -γ*u;

si_infection = ContinuousResourceSharer{Float64}(2, dotsi)
i_recovery = ContinuousResourceSharer{Float64}(1, doti);

undirected_pattern = @relation (S, I) begin
    si_infection(S,I)
    i_recovery(I)
end

to_graphviz(undirected_pattern, box_labels = :name, junction_labels = :variable, edge_attrs=Dict(:len => ".75"))

undirected_system = oapply(undirected_pattern, [si_infection, i_recovery]);

undirected_prob = ODEProblem(undirected_system, u0, tspan)
undirected_sol = solve(undirected_prob,FRK65(0));

plot(undirected_sol)

μ = 1.0/10
dotis(u,p,t) = [-μ*u[1],μ*u[1]];

is_birthdeath = ContinuousResourceSharer{Float64}(2, dotis);

undirected_open_pattern = @relation (S, I) begin
    si_infection(S,I)
    i_recovery(I)
    is_birthdeath(I,S)
end

to_graphviz(undirected_open_pattern, box_labels = :name, junction_labels = :variable, edge_attrs=Dict(:len => ".75"))

undirected_open_system = oapply(undirected_open_pattern, [si_infection, i_recovery, is_birthdeath]);

undirected_open_prob = ODEProblem(undirected_open_system, u0, tspan)
undirected_open_sol = solve(undirected_open_prob,FRK65(0));

plot(undirected_open_sol)

nstages = 4
δ = nstages*γ;

sub(i::Int) = i<0 ? error("$i is negative") : join('₀'+d for d in reverse(digits(i)))
sub(x::String,i::Int) = x*sub(i)
istages = [sub("I",i) for i=1:nstages]

dotsii(u,p,t) = [-β*u[1]*u[3],β*u[1]*u[3],0.0]
dotii(u,p,t) = [-δ*u[1],δ*u[1]]
dotilast(u,p,t) = [-δ*u[1]];

sii_infection = ContinuousResourceSharer{Float64}(3, dotsii)
i_transition = ContinuousResourceSharer{Float64}(2, dotii)
ilast_recovery = ContinuousResourceSharer{Float64}(1, dotilast);

undirected_pattern_stages = @relation (S, I₁, I₂, I₃, I₄) begin
    si_infection(S,I₁)
    sii_infection(S,I₁,I₂)
    sii_infection(S,I₁,I₃)
    sii_infection(S,I₁,I₄)
    i_transition(I₁,I₂)
    i_transition(I₂,I₃)
    i_transition(I₃,I₄)
    ilast_recovery(I₄)
end

to_graphviz(undirected_pattern_stages, box_labels = :name, junction_labels = :variable, edge_attrs=Dict(:len => ".75"))

undirected_system_stages = oapply(undirected_pattern_stages, [
    si_infection
    sii_infection
    sii_infection
    sii_infection
    i_transition
    i_transition
    i_transition
    ilast_recovery])

u0stages = [990.0,10.0,0.0,0.0,0.0];

undirected_stages_prob = ODEProblem(undirected_system_stages, u0stages, tspan)
undirected_stages_sol = solve(undirected_stages_prob,FRK65(0));

undirected_stages_df = DataFrame(undirected_stages_sol(t)')
rename!(undirected_stages_df,["S";istages])
undirected_stages_df[!,:I] = undirected_stages_df[!,:I₁] +
                              undirected_stages_df[!,:I₂] +
                              undirected_stages_df[!,:I₃] +
                              undirected_stages_df[!,:I₄]
undirected_stages_df[:t] = t;

plot(undirected_stages_df[!,:t],
     [undirected_stages_df[!,:S],undirected_stages_df[!,:I]])
plot!(undirected_sol)

max_stars_repo_name: mehalter/sir-julia

max_issues_repo_name: mehalter/sir-julia

max_forks_repo_name: mehalter/sir-julia


text:
# following Joe Fowler's drift_correct from mass

    # """Compute a drift correction that minimizes the spectral entropy of
    # <uncorrected> (a filtered pulse height vector) after correction.
    # The <indicator> vector should be of equal length and should have some
    # linear correlation with the pulse gain. Generally it will be the
    # pretrigger mean of the pulses, but you can experiment with other choices.

    # The entropy will be computed on corrected values only in the range
    # [0, <limit>], so <limit> should be set to a characteristic large
    # value of <uncorrected>. If <limit> is None (the default), then
    # in will be compute as 25% larger than the 99%ile point of
    # <uncorrected>

    # The model is that the filtered pulse height PH should be scaled by
    # (1 + a*PTM) where a is an arbitrary parameter computed here, and
    # PTM is the difference between each record's pretrigger mean and the
    # median value of all pretrigger means. (Or replace "pretrigger mean"
    # with whatever quantity you passed in as <indicator>.)
    # """
using KernelDensity, Distributions, Optim

abstract PerPulseTransformation

immutable LinearDriftCorrect <: PerPulseTransformation
    indicator_median::Float64
    param::Float64
end
uvhist(edges, y) = UnivariateKDE(edges, convert(Array{Float64,1},hist(y,edges)[2]))
density(h::UnivariateKDE) = h.density
edges(h::UnivariateKDE) = h.x
Base.midpoints(h::UnivariateKDE) = midpoints(edges(h))
Base.conv(h::Histogram, y) = conv(UnivariateKDE(edges(h), convert(Vector{Float64},counts(h))), y)
function drift_correct(indicator, uncorrected, σ_smooth, limit=NaN)
    median_indicator = median(indicator)
    centered_indicator = indicator-median_indicator
    if isnan(limit)
        limit = 1.25*StatsBase.percentile(uncorrected, 99)
    end
    edges = pow2edges(σ_smooth*0.4, (0, limit))
    smoother = Normal(0,σ_smooth)
    function tooptimize_entropy(param)
        corrected = uncorrected+centered_indicator*param[1]
        uvh = uvhist(edges, corrected)
        entropy(uvh,σ_smooth)
    end
    optimum = optimize(tooptimize_entropy,-1.0,1.0) # optimize(f(x), lower_bound_x, upper_bound_x) minimizes f with respect to x using brent method
    #check for sucess?
    LinearDriftCorrect(median_indicator, optimum.minimum)
end

function pow2edges(stepsize_guess, limits)
    nbins_guess = int(round( 0.5+(last(limits)-first(limits))/stepsize_guess ))
    pow2=1024
    while pow2<nbins_guess
        pow2*=2
    end
    nbins = pow2
    edges = first(limits):(last(limits)-first(limits))/pow2:last(limits)
end
function entropy(h::UnivariateKDE, σ)
    hsmooth=conv(h, Normal(0,σ))
    -sum(StatsBase.xlogx(density(hsmooth)))
end

using Base.Test

@test length(pow2edges(0.1,1:100))==1025
@test isapprox(first(pow2edges(0.1,1:100)),1)
@test isapprox(last(pow2edges(0.1,1:100)),100)

n=10000
indicator = randn(n)*10
uncorrected = 1000+randn(n)+indicator*0.1
dc = drift_correct(indicator, uncorrected, 0.5)
@test abs(dc.indicator_median)<0.5
@test abs(dc.param+0.1)<0.005


text:
using Distances
using Clustering
using StatsPlots

function hclustplot(hc::Hclust, useheight::Bool)
    o = enumerate(hc.order)
    n = [x for x in 1:length(o)]

    pos = treepositions(hc, useheight)

    xs = []
    ys = []
    for i in 1: size(hc.merge, 1)
        x1 = pos[hc.merge[i,1]][1]
        x2 = pos[hc.merge[i,2]][1]
        append!(xs, [x1,x1,x2,x2])

        y1 = pos[hc.merge[i,1]][2]
        y2 = pos[hc.merge[i,2]][2]
        useheight ? h = hc.height[i] : h = 1
        newy = maximum([y1,y2]) + h
        append!(ys, [y1,newy,newy,y2])
    end
    ret = (reshape(xs, 4, size(hc.merge, 1)), reshape(ys, 4, size(hc.merge, 1)))
    ret = convert.(Array{Float64,2}, ret)
    return ret
end

function treepositions(hc::Hclust, useheight::Bool)
    order = enumerate(hc.order)
    positions = Dict{}()
    for (k,v) in order
        positions[-k] = (v, 0)
    end
    for i in 1:size(hc.merge,1)
        xpos = mean([positions[hc.merge[i,1]][1], positions[hc.merge[i,2]][1]])
        if hc.merge[i,1] < 0 && hc.merge[i,2] < 0
            useheight ? ypos = hc.height[i] : ypos = 1
        else
            useheight ? h = hc.height[i] : h = 1
            ypos = maximum([positions[hc.merge[i,1]][2], positions[hc.merge[i,2]][2]]) + h
        end

        positions[i] = (xpos, ypos)
    end
    return positions
end

max_stars_repo_name: SuReLI/RL-STDP

max_issues_repo_name: SuReLI/RL-STDP

max_forks_repo_name: SuReLI/RL-STDP


text:
#= Here I really would like to start a calculation of the IC50,
because in theory Julia it shoul be very good for this? =#
using Pkg
using Queryverse

println("Function for calculation of the linear curve associated with the data points")
#I would like to make a function that is asking the user to save the data
data = readlines("pippo.txt")
data = load("2018-10-06__Paper_table.xlsx", "Jennifer Keiser Author   Full ") |> DataFrame
pipo = data |> @filter(_.year==2018)
save("gigio.csv", delim=';', pipo)

using Queryverse
#I will make a function that is asking for the DataFrame of the effect
calling_data_frame = function()
    #entering the file name and all the basic thigs that I would need from this function
    println("Please enter the file name, complete of extension")
    file_toload = readline(stdin)

    #Now staring playing with the file just loaded
    #Loading the dataset now and globally saving the dataset
    global data_df = load(file_toload) |> DataFrame
end


text:
# Series evaluation of:
#    N_n/(4*b*r)^{n/2} = \int_{-\kappa/2}^{\kappa/2} (k^2-\sin^2{x})^{n/2} sin^2{x} dx
# for the top two values of n = n_max-1 to n_max.
function Nn_series!(t::Transit_Struct{T}) where {T <: Real}
# Use series expansion to compute M_n:
k2 = t.k2; n_max = t.n_max
# Computing leading coefficient (n=0):
if k2 < 1
  tol = eps(k2); term = zero(T)
  fac = t.sqonembmr2^(n_max-1)*t.k*t.k2
# Now, compute higher order terms until desired precision is reached:
  @inbounds for j=1:2
    # Add leading term to M_n:
    t.Nn[n_max-1+j] = t.Nn_coeff[j,1]
    k2n = one(T)   # k^{2n}
    @inbounds for n =1:t.jmax-1
      k2n *= k2
      term = k2n*t.Nn_coeff[j,n+1]
      t.Nn[n_max-1+j] += term
      if abs(term) < tol
        break
      end
    end
    t.Nn[n_max-1+j] *= fac
    fac *= t.sqonembmr2
  end
  return
else # k^2 >= 1
  # No need for this one.
end
end

# Recursive computation of N_n starting at m=0 to 1, and raising to n_max:
function Nn_raise!(t::Transit_Struct{T}) where {T <: Real}
r=t.r; b=t.b; n_max=t.n_max; k2 = t.k2; kc = t.kc; t.k2inv = inv(t.k2); sqarea=t.sqarea; onemr2mb2=t.onemr2mb2
# Computes the integrals:
# N_n(r,b) = (4*b*r)^(n/2) \int_{-\kappa/2}^{\kappa/2} (k^2 - \sin^2{x})^{n/2} \sin^2{x} dx
# where k^2 = (1-(r-b)^2)/(4*b*r), \kappa/2 = \sin^{-1}(k) for k<1; otherwise \kappa=\pi.
# Compute first two integrals:
Nn_two!(t)
@inbounds for m=2:n_max
  t.Nn[m+1]=(t.Mn[m+1]+m*t.onembpr2*t.Nn[m-1])*t.ninv[m+2]
end
return
end

# Compute N_n for m=0 to 1:
function Nn_two!(t::Transit_Struct{T}) where {T <: Real}
r=t.r; b=t.b; n_max=t.n_max; k2 = t.k2; kc = t.kc; t.k2inv = inv(t.k2)
# Computes the integrals:
# N_n(r,b) = (4*b*r)^(m/2) \int_{-\kappa/2}^{\kappa/2} (k^2 - \sin^2{x})^{m/2} \sin^2{x} dx
# where k^2 = (1-(r-b)^2)/(4*b*r), \kappa/2 = \sin^{-1}(k) for k<1; otherwise \kappa=\pi.
# Compute first two integrals:
if k2 <= 1.0
  t.Nn[1] = 0.5*t.kap0-t.k*t.kc
  t.Nn[2] = t.twothird*2*t.sqbr*t.k2*(-t.Eofk+2*t.Em1mKdm)
else
  t.Nn[1] = 0.5*pi
  t.Nn[2] =  t.twothird*2*t.sqbr*t.k*(2*t.Eofk - t.Em1mKdm)
end
return
end

# Recursive computation of N_n starting at m=n_max to n_max-1, and lowering to m=1:
function Nn_lower!(t::Transit_Struct{T}) where {T <: Real}
r=t.r; b=t.b; n_max=t.n_max; k2 = t.k2
# Computes the integrals:
# N_n(r,b) = (4*b*r)^(m/2) \int_{-\kappa/2}^{\kappa/2} (k^2 - \sin^2{x})^{m/2} \sin^2{x} dx
# where k^2 = (1-(r-b)^2)/(4*b*r), \kappa/2 = \sin^{-1}(k) for k<1; otherwise \kappa=\pi.
# Compute series version for top two:
Nn_series!(t)
# Now iterate downwards:
invonembpr2 = inv(t.onembpr2)
@inbounds for m=n_max-2:-1:2
  t.Nn[m+1]=((m+4)*t.Nn[m+3]-t.Mn[m+3])*invonembpr2*t.ninv[m+2]
end
# Now, compute lowest two exactly:
Nn_two!(t)
return
end


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

mutable struct ApplicationGatewayRewriteRuleActionSet <: SwaggerModel
    requestHeaderConfigurations::Any # spec type: Union{ Nothing, Vector{ApplicationGatewayHeaderConfiguration} } # spec name: requestHeaderConfigurations
    responseHeaderConfigurations::Any # spec type: Union{ Nothing, Vector{ApplicationGatewayHeaderConfiguration} } # spec name: responseHeaderConfigurations
    urlConfiguration::Any # spec type: Union{ Nothing, ApplicationGatewayUrlConfiguration } # spec name: urlConfiguration

    function ApplicationGatewayRewriteRuleActionSet(;requestHeaderConfigurations=nothing, responseHeaderConfigurations=nothing, urlConfiguration=nothing)
        o = new()
        validate_property(ApplicationGatewayRewriteRuleActionSet, Symbol("requestHeaderConfigurations"), requestHeaderConfigurations)
        setfield!(o, Symbol("requestHeaderConfigurations"), requestHeaderConfigurations)
        validate_property(ApplicationGatewayRewriteRuleActionSet, Symbol("responseHeaderConfigurations"), responseHeaderConfigurations)
        setfield!(o, Symbol("responseHeaderConfigurations"), responseHeaderConfigurations)
        validate_property(ApplicationGatewayRewriteRuleActionSet, Symbol("urlConfiguration"), urlConfiguration)
        setfield!(o, Symbol("urlConfiguration"), urlConfiguration)
        o
    end
end # type ApplicationGatewayRewriteRuleActionSet

const _property_map_ApplicationGatewayRewriteRuleActionSet = Dict{Symbol,Symbol}(Symbol("requestHeaderConfigurations")=>Symbol("requestHeaderConfigurations"), Symbol("responseHeaderConfigurations")=>Symbol("responseHeaderConfigurations"), Symbol("urlConfiguration")=>Symbol("urlConfiguration"))
const _property_types_ApplicationGatewayRewriteRuleActionSet = Dict{Symbol,String}(Symbol("requestHeaderConfigurations")=>"Vector{ApplicationGatewayHeaderConfiguration}", Symbol("responseHeaderConfigurations")=>"Vector{ApplicationGatewayHeaderConfiguration}", Symbol("urlConfiguration")=>"ApplicationGatewayUrlConfiguration")
Base.propertynames(::Type{ ApplicationGatewayRewriteRuleActionSet }) = collect(keys(_property_map_ApplicationGatewayRewriteRuleActionSet))
Swagger.property_type(::Type{ ApplicationGatewayRewriteRuleActionSet }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_ApplicationGatewayRewriteRuleActionSet[name]))}
Swagger.field_name(::Type{ ApplicationGatewayRewriteRuleActionSet }, property_name::Symbol) =  _property_map_ApplicationGatewayRewriteRuleActionSet[property_name]

function check_required(o::ApplicationGatewayRewriteRuleActionSet)
    true
end

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


text:
"""
$(TYPEDEF)
$(TYPEDFIELDS)

One Petri Net of a PNML model.
"""
struct PnmlNet{PNTD<:PnmlType,D}
    type::PNTD #TODO
    id::Symbol
    pages::Vector{Page}
    declaration::D

    name::Maybe{Name}
    com::ObjectCommon
    xml::XMLNode
end

function PnmlNet(pntd::PNTD, id::Symbol, pages, declare, name, oc::ObjectCommon, xml::XMLNode) where {PNTD<:PnmlType}
    PnmlNet{PNTD, typeof(declare)}(pntd, id, pages, declare, name, oc, xml)
end

pid(net::PnmlNet) = net.id
pages(net::PnmlNet) = net.pages
declarations(net::PnmlNet) = declarations(net.declaration)

has_labels(net::PnmlNet) = has_labels(net.com)
xmlnode(net::PnmlNet) = net.xml

has_name(net::PnmlNet) = !isnothing(net.name)
name(net::PnmlNet) = has_name(net) ? net.name.text : ""

"Usually the only interesting page."
firstpage(net::PnmlNet) = first(net.pages)


text:
# Warning: setupLLCgrid no longer exists
# Example: read monthly average temperature/salinity fields on poseidon

using Revise
using ECCOtour, MeshArrays, MITgcmTools

# other useful packages in Reemergence project
# using Statistics, PyPlot, Distributions, FFTW, LinearAlgebra, StatsBase
path_grid="../inputs/GRID_LLC90/"
γ = setupLLCgrid(path_grid)

path = "/batou/ECCOv4r4/MITgcm/exps/iter129_bulkformula/run/diags/"
fileroot = "state_3d_set1" # contains θ, S, dρ/dz, see accompanying meta file

# get a list of all files in directory.
filelist = searchdir(path,fileroot) # first filter for state_3d_set1

# just save the list of ".data" files, eliminate ".meta"
datafilelist  = filter(x -> occursin("data",x),filelist) # second filter for "data"

# just save the list of ".meta" files, eliminate ".data"
metafilelist  = filter(x -> occursin("meta",x),filelist) # second filter for "meta"
meta = read_meta(path,metafilelist[1]);
pairs([meta[1]])
[meta[i].dimList for i in 1:length(meta)]

# read the first one into a 4D array: lon, lat, depth, variable type
@time state = read_mdsio(path,datafilelist[1]); # 1.78 secs, MeshArrays reading method is faster

# read into gcmarray
@time state = γ.read(path*datafilelist[1],MeshArray(γ,Float32,nz*nc));

# read into a dictionary that handles the meta data
vars = read_state_3d(path,datafile,γ)


text:
"""
checkinputs \\

for checking that the inputs are consistent \\
Arguments:  x = random variable, usually input \\
            y = random variable, usually output/classification \\
            xmax = the maximum number of possible values across the dimensions in x \\
            ymax = the maximum number of possible values across the dimensions in y \\
            xn = number of variables in x \\
            yn = number of variables in y \\
Returns:    nothing
"""

function checkinputs(x::Array{Int}, y::Array{Int}, xmax::Int, ymax::Int, xn::Int, yn::Int)
  maximum(x) >= xmax && error("x values must be between 0 and ", string(xmax), ".")
  maximum(y) >= ymax && error("x values must be between 0 and ", string(ymax), ".")
  size(x)[1] != xn && error("size(x)[1] must equal xn.")
  size(y)[1] != yn && error("size(y)[1] must equal yn.")
  size(x)[2] != size(y)[2] && error("x and y must have the same number of trials.")
  return nothing
end

function checkinputs(x::Array{Int}, xmax::Int, xn::Int)
  maximum(x) >= xmax && error("x values must be between 0 and ", string(xmax), ".")
  size(x)[1] != xn && error("size(x)[1] must equal xn.")
  return nothing
end


text:
using Vectorize

if VERSION >= v"0.5-"
    using Base.Test
else
    using BaseTestNext
    const Test = BaseTestNext
end

println("===== Testing Vectorize Benchmarked Functions =====")
# Length of test array
N = 1000
srand(13)

# Real Arithmetic
# for T in [Float32, Float64]
#     @testset "Vectorize: Basic Arithmetic::$T" begin
#         X = convert(Vector{T}, abs(randn(N)))
#         Y = convert(Vector{T}, abs(randn(N)))
#         @testset "Testing $f::$T" for (f, fb) in [(:add, .+), (:sub, .-), (:mul, .*), (:div, ./),
#                                                   (:hypot, :hypot), (:atan2, :atan2), (:pow, .^)]
#             @eval fbase = $fb
#             try  # need to catch functions that aren't available from any framework
#                 @eval fvec = Vectorize.$f
#             catch err
#                 println("UndefVar on $f")
#                 if !(isa(err, UndefVarError))
#                     error("Testing $(fvec) has failed")
#                 end

#             end

#             @test fbase(X, Y) ≈ fvec(X, Y)
#             println("Testing $f")
#         end
#     end
# end

# # Complex Arithmetic
# for T in [Complex{Float32}, Complex{Float64}]
#     @testset "Vectorize: Complex Arithmetic::$T" begin
#         X = convert(Vector{T}, randn(N)) + convert(Vector{T}, randn(N))*im
#         Y = convert(Vector{T}, randn(N)) + convert(Vector{T}, randn(N))*im
#         @testset "Testing $f::$T" for (f, fb) in [(:add, .+), (:sub, .-), (:div, ./), (:mul, .*),
#                                                   (:pow, .^)]
#             @eval fbase = $fb

#             try # need to catch functions that aren't available from any framework
#                 @eval fvec = Vectorize.$f
#             catch err
#                 if !(isa(err, UndefVarError))
#                     error("Testing $(fvec) has failed")
#                 end
#                 continue
#             end
#             @test fbase(X, Y) ≈ fvec(X, Y)
#         end

#         @testset "Testing $f::$T" for f in [:mulbyconj]

#             try # need to catch functions that aren't available from any framework
#                 @eval fvec = Vectorize.$f
#             catch err
#                 if !(isa(err, UndefVarError))
#                     error("Testing $(fvec) has failed")
#                 end
#                 continue
#             end
#             @test X .* conj(Y) ≈ fvec(X, Y)
#         end
#     end
# end

# # Real Vector Math
# for T in [Float32, Float64]
#     @testset "Vectorize: Real Vector Math::$T" begin
#         X = convert(Vector{T}, randn(N)) # unrestricted
#         Y = convert(Vector{T}, abs(randn(N))) .+ 1 # greater than 1
#         Z = clamp(convert(Vector{T}, randn(N)), -1,  1) # [-1, 1]
#         W = clamp(convert(Vector{T}, randn(N)), 0,  2) # [0, 2]
#         @testset "Testing $f::$T" for f in [:exp, :abs, :ceil, :floor, :round,
#                                             :trunc, :cos, :sin, :tan, :cosh, :sinh, :tanh,
#                                             :cbrt, :erf, :erfc, :lgamma, :gamma]
#             @eval fb = $f
#             @eval fvec = Vectorize.$f
#             # need to catch functions that aren't available from any framework
#             if testexist(fvec, X) == true
#                 @test fb(X) ≈ fvec(X)
#             end
#         end

#         @testset "Testing $f::$T" for f in [:sqrt, :log10, :acosh, :asinh, :log]
#             @eval fb = $f
#             @eval fvec = Vectorize.$f
#             # need to catch functions that aren't available from any framework
#             if testexist(fvec, Y) == true
#                 @test fb(Y) ≈ fvec(Y)
#             end
#         end

#         # @testset "Testing $f::$T" for f in [:erfcinv]
#         #     @eval fb = $f
#         #     @eval fvec = Vectorize.$f
#         #     @test fb(W) ≈ fvec(W)
#         # end

#         @testset "Testing $f::$T" for f in [:acos, :asin, :atan]#, :atanh] #, :erfinv]
#             @eval fb = $f
#             @eval fvec = Vectorize.$f
#             # need to catch functions that aren't available from any framework
#             if testexist(fvec, Z) == true
#                 @test fb(Z) ≈ fvec(Z)
#             end
#         end
#     end
# end

# # Complex Vector Math
# for T in [Complex{Float32}, Complex{Float64}]
#     @testset "Vectorize: Complex Vector Functions::$T" begin
#         X = convert(Vector{T}, complex(randn(N)))
#         @testset "Testing $f::$T" for f in [:exp, :abs, :cos, :sin, :tan, :cosh, :sinh, :tanh,
#                                             :sqrt, :log10, :acosh, :asinh, :log, :acos, :asin,
#                                             :atan, :conj] #, :atanh]
#             @eval fb = $f
#             @eval fvec = Vectorize.$f
#             if testexist(fvec, X) == true
#                 @test fb(X) ≈ fvec(X)
#             end
#         end
#     end
# end

println("===== Vectorize Tests Successful =====\n\n")


text:
using OpenCL

str_kernel_test_01 = "

__kernel void fwsim(__global const float *semus,
                     __global const float *fluxes,
                     __global float *mds_out)
   {
     const int gid = get_group_id(0);
     const int idx = get_local_id(0);
     const int idy = get_local_id(1);

     //int offset_semus   = 0;
     //int offset_fluxes = gid * 17;
     //int offset_out    = gid * 3;

     // Here we store the sizes of the given matrices in a step
     int mAm = 0;
     int mAn = 0;
     int mBm = 0;
     int mBn = 0;
     int mXm = 0;
     int mXn = 0;
     int mYm = 0;
     int mYn = 0;

__local float mA[25];
__local float mA_L[25];
__local float mB[10];
__local float mX[12];
__local float mY[8];
__local float mBY[12];

//Required EMRs: 11  Required Fluxes: 6  Target MDs: 4
int offset_emrs   = 0;
int offset_fluxes = gid * 6;
int offset_target = gid * 4;
//-----Init EMRs-------------------------------------------------------------------------------------
float mr_sA_2_M_0 = semus[offset_emrs + 1];
float mr_sA_2_M_1 = semus[offset_emrs + 2];
float mr_sA_3_M_0 = semus[offset_emrs + 3];
float mr_sA_3_M_1 = semus[offset_emrs + 4];
float mr_sA_2_3_M_0 = semus[offset_emrs + 5];
float mr_sA_2_3_M_1 = semus[offset_emrs + 6];
float mr_sA_2_3_M_2 = semus[offset_emrs + 7];
float mr_sA_1_2_3_M_0 = semus[offset_emrs + 8];
float mr_sA_1_2_3_M_1 = semus[offset_emrs + 9];
float mr_sA_1_2_3_M_2 = semus[offset_emrs + 10];
float mr_sA_1_2_3_M_3 = semus[offset_emrs + 11];
//-----Init Fluxes-----------------------------------------------------------------------------------
float v_R01 = fluxes[ offset_fluxes + 1];
float v_R02 = fluxes[ offset_fluxes + 2];
float v_R03 = fluxes[ offset_fluxes + 3];
float v_R04 = fluxes[ offset_fluxes + 4];
float v_R05 = fluxes[ offset_fluxes + 5];
float v_R06 = fluxes[ offset_fluxes + 6];
//---------------------------------------------------------------------------------------------------
//------Declare All Temporary EMR Float values--------------------------------------------------------------
float mr_sC_1_M_1;
float mr_sD_3_M_0;
float mr_sB_1_2_3_M_1;
float mr_sD_2_M_1;
float mr_sF_1_2_3_M_1;
float mr_sD_2_3_M_0;
float mr_sD_1_2_3_M_0;
float mr_sF_1_2_3_M_2;
float mr_sB_3_M_1;
float mr_sB_2_3_M_0;
float mr_sB_1_2_3_M_2;
float mr_sB_2_M_0;
float mr_sD_1_2_3_M_3;
float mr_sB_2_M_1;
float mr_sD_1_2_3_M_1;
float mr_sD_2_3_M_1;
float mr_sB_3_M_0;
float mr_sB_2_3_M_1;
float mr_sF_1_2_3_M_0;
float mr_sD_2_M_0;
float mr_sB_1_2_3_M_3;
float mr_sB_2_3_M_2;
float mr_sD_2_3_M_2;
float mr_sD_1_2_3_M_2;
float mr_sC_1_M_0;
float mr_sF_1_2_3_M_3;
float mr_sD_3_M_1;
float mr_sB_1_2_3_M_0;
//---------------------------------------------------------------------------------------------------
//----Init A---------------------------------------------------------------------------------
if(idx==0 && idy==0){
mAm = 5;
mAn = 5;
mA[0] = ( (  -1.0f  * v_R01   )  +(  -1.0f  * v_R03   )   ) ;
mA[1] = 0.0f;
mA[2] =  v_R04 ;
mA[3] =  v_R02 ;
mA[4] = 0.0f;
mA[5] = 0.0f;
mA[6] = ( (  -1.0f  * v_R01   )  +(  -1.0f  * v_R03   )   ) ;
mA[7] = 0.0f;
mA[8] =  v_R05 ;
mA[9] =  v_R02 ;
mA[10] = 0.0f;
mA[11] = 0.0f;
mA[12] = ( (  -1.0f  * v_R04   )   ) ;
mA[13] = 0.0f;
mA[14] =  v_R05 ;
mA[15] =  v_R03 ;
mA[16] = 0.0f;
mA[17] = 0.0f;
mA[18] = ( (  -1.0f  * v_R02   )  +(  -1.0f  * v_R05   )   ) ;
mA[19] = 0.0f;
mA[20] = 0.0f;
mA[21] =  v_R03 ;
mA[22] = 0.0f;
mA[23] = 0.0f;
mA[24] = ( (  -1.0f  * v_R05   )  +(  -1.0f  * v_R02   )   ) ;
}
//----Init B---------------------------------------------------------------------------
if(idx==1 && idy==0){
mBm = 5;
mBn = 2;
mB[0] = 0.0f;
mB[1] = (  -1.0f  * v_R01   ) ;
mB[2] = 0.0f;
mB[3] = 0.0f;
mB[4] = 0.0f;
mB[5] = (  -1.0f  * v_R01   ) ;
mB[6] = 0.0f;
mB[7] = 0.0f;
mB[8] = 0.0f;
mB[9] = 0.0f;
}
//----Init Y---------------------------------------------------------------------------
if(idx==2 && idy==0){
mYm = 2;
mYn = 2;
mY[0] = ( (  mr_sA_3_M_0   )   ) ;
mY[1] = ( (  mr_sA_2_M_0   )   ) ;
mY[2] = ( (  mr_sA_3_M_1   )   ) ;
mY[3] = ( (  mr_sA_2_M_1   )   ) ;
}

for(int i=0;i<mAm*mAn;i++){ printf(\" %f \",mA[i]); }
printf(\" \\n \")
for(int i=0;i<mBm*mBn;i++){ printf(\" %f \",mB[i]); }
printf(\" \\n \")
for(int i=0;i<mYm*mYn;i++){ printf(\" %f \",mY[i]); }
printf(\" \\n \")

//------COMPUTE!-----------------------------------------------------------------------
//LU-DECOMPOSITION (in place): mA_L * mA = mA
if( idx<5 && idy<5 ){ mA_L[ idx+5*idy ] = (idx==idy)?1.0f:0.0f; }
barrier(CLK_LOCAL_MEM_FENCE);
//LU-DECOMPOSITION STEP 0 / 4
if( idx==0 && idy > 0 && idy < 5 ){
mA_L[ idy + 0 ] = mA[ idy + 0 ] / mA[ 0 ];
}
barrier(CLK_LOCAL_MEM_FENCE);
if( idx>0 && idx < 5 && idy < 5 ){
mA[ idx + 5*idy ] -= mA_L[ idx + 0 ] * mA[0+ idy*5 ];
}
barrier(CLK_LOCAL_MEM_FENCE);
//LU-DECOMPOSITION STEP 1 / 4
if( idx==0 && idy > 1 && idy < 5 ){
mA_L[ idy + 5 ] = mA[ idy + 5 ] / mA[ 6 ];
}
barrier(CLK_LOCAL_MEM_FENCE);
if( idx>1 && idx < 5 && idy < 5 ){
mA[ idx + 5*idy ] -= mA_L[ idx + 5 ] * mA[1+ idy*5 ];
}
barrier(CLK_LOCAL_MEM_FENCE);
//LU-DECOMPOSITION STEP 2 / 4
if( idx==0 && idy > 2 && idy < 5 ){
mA_L[ idy + 10 ] = mA[ idy + 10 ] / mA[ 12 ];
}
barrier(CLK_LOCAL_MEM_FENCE);
if( idx>2 && idx < 5 && idy < 5 ){
mA[ idx + 5*idy ] -= mA_L[ idx + 10 ] * mA[2+ idy*5 ];
}
barrier(CLK_LOCAL_MEM_FENCE);
//LU-DECOMPOSITION STEP 3 / 4
if( idx==0 && idy > 3 && idy < 5 ){
mA_L[ idy + 15 ] = mA[ idy + 15 ] / mA[ 18 ];
}
barrier(CLK_LOCAL_MEM_FENCE);
if( idx>3 && idx < 5 && idy < 5 ){
mA[ idx + 5*idy ] -= mA_L[ idx + 15 ] * mA[3+ idy*5 ];
}
barrier(CLK_LOCAL_MEM_FENCE);
//LU-DECOMPOSITION STEP 4 / 4
if( idx==0 && idy > 4 && idy < 5 ){
mA_L[ idy + 20 ] = mA[ idy + 20 ] / mA[ 24 ];
}
barrier(CLK_LOCAL_MEM_FENCE);
if( idx>4 && idx < 5 && idy < 5 ){
mA[ idx + 5*idy ] -= mA_L[ idx + 20 ] * mA[4+ idy*5 ];
}
barrier(CLK_LOCAL_MEM_FENCE);
//LU-DECOMPOSITION DONE!
//----------------------------------
//MATRIX MUL: mBY = mB * mY
//Work item (idx,idy) compute mBY[idx,idy]
if( idx < 5 && idy < 2 ) {
float mulsum = 0.0f;
mulsum += mB[ 5*0 + idx ] * mY[ idy * 2 + 0 ];
mulsum += mB[ 5*1 + idx ] * mY[ idy * 2 + 1 ];
mBY[idx+5*idy] = mulsum;
}
//MATRIX MUL COMPLETE!
barrier(CLK_LOCAL_MEM_FENCE);
// Solve L System
// Solve L System - Row 0
if(idx==0 && idy<2 ){
float sum_i = 0.0f;
mBY[ 0 + 5*idy ] = (mBY[ 0 + 5*idy ] - sum_i ) / mA_L[ 0 ] ;
}
// Solve L System - Row 1
if(idx==0 && idy<2 ){
float sum_i = 0.0f;
sum_i += mBY[ 0 + 5*idy ];
mBY[ 1 + 5*idy ] = (mBY[ 1 + 5*idy ] - sum_i ) / mA_L[ 6 ] ;
}
// Solve L System - Row 2
if(idx==0 && idy<2 ){
float sum_i = 0.0f;
sum_i += mBY[ 0 + 5*idy ];
sum_i += mBY[ 1 + 5*idy ];
mBY[ 2 + 5*idy ] = (mBY[ 2 + 5*idy ] - sum_i ) / mA_L[ 12 ] ;
}
// Solve L System - Row 3
if(idx==0 && idy<2 ){
float sum_i = 0.0f;
sum_i += mBY[ 0 + 5*idy ];
sum_i += mBY[ 1 + 5*idy ];
sum_i += mBY[ 2 + 5*idy ];
mBY[ 3 + 5*idy ] = (mBY[ 3 + 5*idy ] - sum_i ) / mA_L[ 18 ] ;
}
// Solve L System - Row 4
if(idx==0 && idy<2 ){
float sum_i = 0.0f;
sum_i += mBY[ 0 + 5*idy ];
sum_i += mBY[ 1 + 5*idy ];
sum_i += mBY[ 2 + 5*idy ];
sum_i += mBY[ 3 + 5*idy ];
mBY[ 4 + 5*idy ] = (mBY[ 4 + 5*idy ] - sum_i ) / mA_L[ 24 ] ;
}
// Solve L System DONE!
barrier(CLK_LOCAL_MEM_FENCE);
// Solve U System
// Solve U System - Row 4
if(idx==0 && idy<2 ){
float sum_i = 0.0f;
sum_i += mBY[ 5 + 5*idy ];
sum_i += mBY[ 4 + 5*idy ];
sum_i += mBY[ 3 + 5*idy ];
sum_i += mBY[ 2 + 5*idy ];
sum_i += mBY[ 1 + 5*idy ];
mBY[ 4 + 5*idy ] = (mBY[ 4 + 5*idy ] - sum_i ) / mA[ 24 ] ;
}
// Solve U System - Row 3
if(idx==0 && idy<2 ){
float sum_i = 0.0f;
sum_i += mBY[ 5 + 5*idy ];
sum_i += mBY[ 4 + 5*idy ];
sum_i += mBY[ 3 + 5*idy ];
sum_i += mBY[ 2 + 5*idy ];
mBY[ 3 + 5*idy ] = (mBY[ 3 + 5*idy ] - sum_i ) / mA[ 18 ] ;
}
// Solve U System - Row 2
if(idx==0 && idy<2 ){
float sum_i = 0.0f;
sum_i += mBY[ 5 + 5*idy ];
sum_i += mBY[ 4 + 5*idy ];
sum_i += mBY[ 3 + 5*idy ];
mBY[ 2 + 5*idy ] = (mBY[ 2 + 5*idy ] - sum_i ) / mA[ 12 ] ;
}
// Solve U System - Row 1
if(idx==0 && idy<2 ){
float sum_i = 0.0f;
sum_i += mBY[ 5 + 5*idy ];
sum_i += mBY[ 4 + 5*idy ];
mBY[ 1 + 5*idy ] = (mBY[ 1 + 5*idy ] - sum_i ) / mA[ 6 ] ;
}
// Solve U System - Row 0
if(idx==0 && idy<2 ){
float sum_i = 0.0f;
sum_i += mBY[ 5 + 5*idy ];
mBY[ 0 + 5*idy ] = (mBY[ 0 + 5*idy ] - sum_i ) / mA[ 0 ] ;
}
// Solve U System DONE!
barrier(CLK_LOCAL_MEM_FENCE);
//-------------------------------------------------------------------------------------
if(idx==0&&idy==0){
mr_sB_2_M_0 = mBY[ 6 ];
mr_sB_3_M_0 = mBY[ 7 ];
mr_sC_1_M_0 = mBY[ 8 ];
mr_sD_2_M_0 = mBY[ 9 ];
mr_sD_3_M_0 = mBY[ 10 ];
mr_sB_2_M_1 = mBY[ 11 ];
mr_sB_3_M_1 = mBY[ 12 ];
mr_sC_1_M_1 = mBY[ 13 ];
mr_sD_2_M_1 = mBY[ 14 ];
mr_sD_3_M_1 = mBY[ 15 ];
}
barrier(CLK_LOCAL_MEM_FENCE);
//-------------------------------------------------------------------------------------
//----Init A---------------------------------------------------------------------------------
if(idx==0 && idy==0){
mAm = 2;
mAn = 2;
mA[0] = ( (  -1.0f  * v_R01   )  +(  -1.0f  * v_R03   )   ) ;
mA[1] =  v_R02 ;
mA[2] =  v_R03 ;
mA[3] = ( (  -1.0f  * v_R02   )  +(  -1.0f  * v_R05   )   ) ;
}
//----Init B---------------------------------------------------------------------------
if(idx==1 && idy==0){
mBm = 2;
mBn = 2;
mB[0] = (  -1.0f  * v_R01   ) ;
mB[1] = 0.0f;
mB[2] = 0.0f;
mB[3] = (  -1.0f  * v_R05   ) ;
}
//----Init Y---------------------------------------------------------------------------
if(idx==2 && idy==0){
mYm = 2;
mYn = 3;
mY[0] = ( (  mr_sA_2_3_M_0   )   ) ;
mY[1] = ( (  mr_sB_3_M_0  * mr_sC_1_M_0   )   ) ;
mY[2] = ( (  mr_sA_2_3_M_1   )   ) ;
mY[3] = ( (  mr_sB_3_M_0  * mr_sC_1_M_1   )  +(  mr_sB_3_M_1  * mr_sC_1_M_0   )   ) ;
mY[4] = ( (  mr_sA_2_3_M_2   )   ) ;
mY[5] = ( (  mr_sB_3_M_1  * mr_sC_1_M_1   )   ) ;
}

//------COMPUTE!-----------------------------------------------------------------------
//LU-DECOMPOSITION (in place): mA_L * mA = mA
if( idx<2 && idy<2 ){ mA_L[ idx+2*idy ] = (idx==idy)?1.0f:0.0f; }
barrier(CLK_LOCAL_MEM_FENCE);
//LU-DECOMPOSITION STEP 0 / 1
if( idx==0 && idy > 0 && idy < 2 ){
mA_L[ idy + 0 ] = mA[ idy + 0 ] / mA[ 0 ];
}
barrier(CLK_LOCAL_MEM_FENCE);
if( idx>0 && idx < 2 && idy < 2 ){
mA[ idx + 2*idy ] -= mA_L[ idx + 0 ] * mA[0+ idy*2 ];
}
barrier(CLK_LOCAL_MEM_FENCE);
//LU-DECOMPOSITION STEP 1 / 1
if( idx==0 && idy > 1 && idy < 2 ){
mA_L[ idy + 2 ] = mA[ idy + 2 ] / mA[ 3 ];
}
barrier(CLK_LOCAL_MEM_FENCE);
if( idx>1 && idx < 2 && idy < 2 ){
mA[ idx + 2*idy ] -= mA_L[ idx + 2 ] * mA[1+ idy*2 ];
}
barrier(CLK_LOCAL_MEM_FENCE);
//LU-DECOMPOSITION DONE!
//----------------------------------
//MATRIX MUL: mBY = mB * mY
//Work item (idx,idy) compute mBY[idx,idy]
if( idx < 2 && idy < 3 ) {
float mulsum = 0.0f;
mulsum += mB[ 2*0 + idx ] * mY[ idy * 2 + 0 ];
mulsum += mB[ 2*1 + idx ] * mY[ idy * 2 + 1 ];
mBY[idx+2*idy] = mulsum;
}
//MATRIX MUL COMPLETE!
barrier(CLK_LOCAL_MEM_FENCE);
// Solve L System
// Solve L System - Row 0
if(idx==0 && idy<3 ){
float sum_i = 0.0f;
mBY[ 0 + 2*idy ] = (mBY[ 0 + 2*idy ] - sum_i ) / mA_L[ 0 ] ;
}
// Solve L System - Row 1
if(idx==0 && idy<3 ){
float sum_i = 0.0f;
sum_i += mBY[ 0 + 2*idy ];
mBY[ 1 + 2*idy ] = (mBY[ 1 + 2*idy ] - sum_i ) / mA_L[ 3 ] ;
}
// Solve L System DONE!
barrier(CLK_LOCAL_MEM_FENCE);
// Solve U System
// Solve U System - Row 1
if(idx==0 && idy<3 ){
float sum_i = 0.0f;
sum_i += mBY[ 2 + 2*idy ];
sum_i += mBY[ 1 + 2*idy ];
mBY[ 1 + 2*idy ] = (mBY[ 1 + 2*idy ] - sum_i ) / mA[ 3 ] ;
}
// Solve U System - Row 0
if(idx==0 && idy<3 ){
float sum_i = 0.0f;
sum_i += mBY[ 2 + 2*idy ];
mBY[ 0 + 2*idy ] = (mBY[ 0 + 2*idy ] - sum_i ) / mA[ 0 ] ;
}
// Solve U System DONE!
barrier(CLK_LOCAL_MEM_FENCE);
//-------------------------------------------------------------------------------------
if(idx==0&&idy==0){
mr_sB_2_3_M_0 = mBY[ 3 ];
mr_sD_2_3_M_0 = mBY[ 4 ];
mr_sB_2_3_M_1 = mBY[ 5 ];
mr_sD_2_3_M_1 = mBY[ 6 ];
mr_sB_2_3_M_2 = mBY[ 7 ];
mr_sD_2_3_M_2 = mBY[ 8 ];
}
barrier(CLK_LOCAL_MEM_FENCE);
//-------------------------------------------------------------------------------------
//----Init A---------------------------------------------------------------------------------
if(idx==0 && idy==0){
mAm = 3;
mAn = 3;
mA[0] = ( (  -1.0f  * v_R03   )  +(  -1.0f  * v_R01   )   ) ;
mA[1] =  v_R02 ;
mA[2] = 0.0f;
mA[3] =  v_R03 ;
mA[4] = ( (  -1.0f  * v_R05   )  +(  -1.0f  * v_R02   )   ) ;
mA[5] =  v_R06 ;
mA[6] = 0.0f;
mA[7] = 0.0f;
mA[8] = ( (  -1.0f  * v_R06   )   ) ;
}
//----Init B---------------------------------------------------------------------------
if(idx==1 && idy==0){
mBm = 3;
mBn = 2;
mB[0] = 0.0f;
mB[1] = (  -1.0f  * v_R05   ) ;
mB[2] = 0.0f;
mB[3] = (  -1.0f  * v_R01   ) ;
mB[4] = 0.0f;
mB[5] = 0.0f;
}
//----Init Y---------------------------------------------------------------------------
if(idx==2 && idy==0){
mYm = 2;
mYn = 4;
mY[0] = ( (  mr_sB_2_3_M_0  * mr_sC_1_M_0   )   ) ;
mY[1] = ( (  mr_sA_1_2_3_M_0   )   ) ;
mY[2] = ( (  mr_sB_2_3_M_0  * mr_sC_1_M_1   )  +(  mr_sB_2_3_M_1  * mr_sC_1_M_0   )   ) ;
mY[3] = ( (  mr_sA_1_2_3_M_1   )   ) ;
mY[4] = ( (  mr_sB_2_3_M_1  * mr_sC_1_M_1   )  +(  mr_sB_2_3_M_2  * mr_sC_1_M_0   )   ) ;
mY[5] = ( (  mr_sA_1_2_3_M_2   )   ) ;
mY[6] = ( (  mr_sB_2_3_M_2  * mr_sC_1_M_1   )   ) ;
mY[7] = ( (  mr_sA_1_2_3_M_3   )   ) ;
}

//------COMPUTE!-----------------------------------------------------------------------
//LU-DECOMPOSITION (in place): mA_L * mA = mA
if( idx<3 && idy<3 ){ mA_L[ idx+3*idy ] = (idx==idy)?1.0f:0.0f; }
barrier(CLK_LOCAL_MEM_FENCE);
//LU-DECOMPOSITION STEP 0 / 2
if( idx==0 && idy > 0 && idy < 3 ){
mA_L[ idy + 0 ] = mA[ idy + 0 ] / mA[ 0 ];
}
barrier(CLK_LOCAL_MEM_FENCE);
if( idx>0 && idx < 3 && idy < 3 ){
mA[ idx + 3*idy ] -= mA_L[ idx + 0 ] * mA[0+ idy*3 ];
}
barrier(CLK_LOCAL_MEM_FENCE);
//LU-DECOMPOSITION STEP 1 / 2
if( idx==0 && idy > 1 && idy < 3 ){
mA_L[ idy + 3 ] = mA[ idy + 3 ] / mA[ 4 ];
}
barrier(CLK_LOCAL_MEM_FENCE);
if( idx>1 && idx < 3 && idy < 3 ){
mA[ idx + 3*idy ] -= mA_L[ idx + 3 ] * mA[1+ idy*3 ];
}
barrier(CLK_LOCAL_MEM_FENCE);
//LU-DECOMPOSITION STEP 2 / 2
if( idx==0 && idy > 2 && idy < 3 ){
mA_L[ idy + 6 ] = mA[ idy + 6 ] / mA[ 8 ];
}
barrier(CLK_LOCAL_MEM_FENCE);
if( idx>2 && idx < 3 && idy < 3 ){
mA[ idx + 3*idy ] -= mA_L[ idx + 6 ] * mA[2+ idy*3 ];
}
barrier(CLK_LOCAL_MEM_FENCE);
//LU-DECOMPOSITION DONE!
//----------------------------------
//MATRIX MUL: mBY = mB * mY
//Work item (idx,idy) compute mBY[idx,idy]
if( idx < 3 && idy < 4 ) {
float mulsum = 0.0f;
mulsum += mB[ 3*0 + idx ] * mY[ idy * 2 + 0 ];
mulsum += mB[ 3*1 + idx ] * mY[ idy * 2 + 1 ];
mBY[idx+3*idy] = mulsum;
}
//MATRIX MUL COMPLETE!
barrier(CLK_LOCAL_MEM_FENCE);
// Solve L System
// Solve L System - Row 0
if(idx==0 && idy<4 ){
float sum_i = 0.0f;
mBY[ 0 + 3*idy ] = (mBY[ 0 + 3*idy ] - sum_i ) / mA_L[ 0 ] ;
}
// Solve L System - Row 1
if(idx==0 && idy<4 ){
float sum_i = 0.0f;
sum_i += mBY[ 0 + 3*idy ];
mBY[ 1 + 3*idy ] = (mBY[ 1 + 3*idy ] - sum_i ) / mA_L[ 4 ] ;
}
// Solve L System - Row 2
if(idx==0 && idy<4 ){
float sum_i = 0.0f;
sum_i += mBY[ 0 + 3*idy ];
sum_i += mBY[ 1 + 3*idy ];
mBY[ 2 + 3*idy ] = (mBY[ 2 + 3*idy ] - sum_i ) / mA_L[ 8 ] ;
}
// Solve L System DONE!
barrier(CLK_LOCAL_MEM_FENCE);
// Solve U System
// Solve U System - Row 2
if(idx==0 && idy<4 ){
float sum_i = 0.0f;
sum_i += mBY[ 3 + 3*idy ];
sum_i += mBY[ 2 + 3*idy ];
sum_i += mBY[ 1 + 3*idy ];
mBY[ 2 + 3*idy ] = (mBY[ 2 + 3*idy ] - sum_i ) / mA[ 8 ] ;
}
// Solve U System - Row 1
if(idx==0 && idy<4 ){
float sum_i = 0.0f;
sum_i += mBY[ 3 + 3*idy ];
sum_i += mBY[ 2 + 3*idy ];
mBY[ 1 + 3*idy ] = (mBY[ 1 + 3*idy ] - sum_i ) / mA[ 4 ] ;
}
// Solve U System - Row 0
if(idx==0 && idy<4 ){
float sum_i = 0.0f;
sum_i += mBY[ 3 + 3*idy ];
mBY[ 0 + 3*idy ] = (mBY[ 0 + 3*idy ] - sum_i ) / mA[ 0 ] ;
}
// Solve U System DONE!
barrier(CLK_LOCAL_MEM_FENCE);
//-------------------------------------------------------------------------------------
if(idx==0&&idy==0){
mr_sB_1_2_3_M_0 = mBY[ 4 ];
mr_sD_1_2_3_M_0 = mBY[ 5 ];
mr_sF_1_2_3_M_0 = mBY[ 6 ];
mr_sB_1_2_3_M_1 = mBY[ 7 ];
mr_sD_1_2_3_M_1 = mBY[ 8 ];
mr_sF_1_2_3_M_1 = mBY[ 9 ];
mr_sB_1_2_3_M_2 = mBY[ 10 ];
mr_sD_1_2_3_M_2 = mBY[ 11 ];
mr_sF_1_2_3_M_2 = mBY[ 12 ];
mr_sB_1_2_3_M_3 = mBY[ 13 ];
mr_sD_1_2_3_M_3 = mBY[ 14 ];
mr_sF_1_2_3_M_3 = mBY[ 15 ];
}
barrier(CLK_LOCAL_MEM_FENCE);
//-------------------------------------------------------------------------------------
//---RESULTS-----------------------------------------------------------------------------------------
if(idx==0&&idy==0){
mds_out[ offset_target + 1 ] = mr_sF_1_2_3_M_0;
mds_out[ offset_target + 2 ] = mr_sF_1_2_3_M_1;
mds_out[ offset_target + 3 ] = mr_sF_1_2_3_M_2;
mds_out[ offset_target + 4 ] = mr_sF_1_2_3_M_3;
}
//---------------------------------------------------------------------------------------------------
}
"

function test_0()
    a = convert( Array{Float32,1} ,[100.; 110 ; 50 ; 20 ; 20 ; 80 ; 0;0;0;0;0;0;0;0;0;0] )
    b = convert( Array{Float32,1} ,[0.00 ; 1.00; 1.00; 0.00; 0.00; 1.0 ;0.0;0.0;1.0;0.0;0.0;0;0;0;0;0] )

    device, ctx, queue = OpenCL.cl.create_compute_context()

    a_buff = cl.Buffer(Float32, ctx, (:r, :copy), hostbuf=a)
    b_buff = cl.Buffer(Float32, ctx, (:r, :copy), hostbuf=b)
    c_buff = cl.Buffer(Float32, ctx, :w, length(a))

    kernel_str = str_kernel_test_01
    p = cl.Program(ctx, source=kernel_str) |> cl.build!

    #p = cl.Program(ctx, source=kernel_02) |> cl.build!

    #cl.CL_device_info(1)

    k = cl.Kernel(p, "fwsim")
    queue(k, [16;16], [16;16], a_buff, b_buff, c_buff)

    #cl.api.clGetKernelWorkGroupInfo(cl.CL_KERNEL_WORK_GROUP_SIZE)
    r = cl.read(queue, c_buff)

    cl.work_group_info( k , cl.CL_KERNEL_WORK_GROUP_SIZE ,device )
    cl.work_group_info( k , cl.CL_KERNEL_PREFERRED_WORK_GROUP_SIZE_MULTIPLE  ,device )
end

test_0()


text:
include(joinpath(@__DIR__, "..", "src", "RBFBase.jl"))
using .RBF
using ForwardDiff: gradient, hessian, jacobian     # to automatically differentiate the model outputs for comparison
using Test
import FiniteDiff;
fd = FiniteDiff;
m = RBF;

# RBFModel with 2 input and 2 outputs

f1(x) = 1 .+ sum(x.^2);
f2(x) = sin(sum(x));
f(x) = [f1(x), f2(x)];

sites = vec( [[x;y] for x = range( -2, 2, length = 5),
                    y = range( -2, 2, length = 5) ] )

sites = [rand(2) for i = 1 : 5]

vals_1 = f1.(sites);
vals_2 = f2.(sites);
vals = f.(sites);

shape_parameters = [ 1 ]
kernels = [:multiquadric, :exp, :cubic, :thin_plate_spline]

println("Looping through different RBF settings.")
for poly_deg ∈ [ -1, 0, 1 ]
    for kernel ∈ kernels
        for sp ∈ shape_parameters

            A = 1e-3    # somewhat loose tolerance because of possibly bad conditioning
            println("Kernel $(kernel), sp = $sp, poly_deg = $poly_deg")
            rbf = m.RBFModel(
                training_sites = sites,
                training_values = vals,
                kernel = kernel,
                shape_parameter = sp,
                polynomial_degree = poly_deg,
            );
            m.train!(rbf);

            # Test interpolation properties
            rbf1 = x -> m.output( rbf, 1, x );
            rbf2 = x -> m.output( rbf, 2, x );

            evals_1 = rbf1.( sites )
            @test all( isapprox.( evals_1 , vals_1,  atol=A ) )

            evals_2 = rbf2.( sites )
            @test all( isapprox.( evals_2 , vals_2,  atol=A ) )

            rbf_output = x -> m.output( rbf, x );
            @test rbf_output( sites[1] ) == rbf.function_handle( sites[1] )

            # Test gradient calculations
            d_rbf1 = x -> m.grad( rbf, 1, x);
            d_rbf2 = x -> m.grad( rbf, 2, x);

            fwd_rbf1 = χ -> gradient( x -> m.output( rbf, 1, x ), χ);
            fwd_rbf2 = χ -> gradient(x -> m.output( rbf, 2, x ), χ);

            dvals_1 = d_rbf1.(sites)
            dvals_2 = d_rbf2.(sites)
            fwdvals_1 = fwd_rbf1.(sites)
            fwdvals_2 = fwd_rbf2.(sites)

            @test all( isapprox.( dvals_1, fwdvals_1, atol = A) )
            @test all( isapprox.( dvals_2, fwdvals_2, atol = A) )

            jacobian_handle = x -> m.jac( rbf, x )
            @test jacobian_handle(sites[1]) ≈ [d_rbf1(sites[1])'; d_rbf2(sites[1])']

            # test hessian correctness
            h_rbf1 = x -> m.hess( rbf, 1, x )

            H_fwd1 = χ -> fd.finite_difference_hessian( x -> m.output( rbf, 1, x), χ )

            X = rand(2)
            @test h_rbf1( X ) ≈ H_fwd1( X ) atol = A

        end
    end
end


text:
using JupyterSite
using Base.test

@test 1==1


text:
using EventSimulation
using Distributions
using Printf
using Random

# Objectives of the example:
# * show how SimQueue object can be used
# * handling actions by controling state of objects using enums

@enum CUSTOMER_STATE WAITING CASHING RENEGED FINISHED

mutable struct Customer
    name::String
    state::CUSTOMER_STATE
    arrival_time::Float64
end

mutable struct BankState <: AbstractState
    all_cust::Int # how many customers we want to process
    n::Int # how many customers were processed
    reneged_count::Int # number of customers that reneged
    arrival::Exponential{Float64}
    patience::Uniform{Float64}
    cashing::Exponential{Float64}
    counter::SimQueue{Customer}
    report::Bool
end

function source(s::Scheduler)
    s.state.n += 1
    c = Customer(string("Customer", s.state.n), WAITING, s.now)
    s.state.report && @printf("%6.4f %s: Here I am\n", s.now, c.name)
    provide!(s, s.state.counter, c)
    register!(s, x -> renege(x, c), rand(s.state.patience))
    if s.state.n < s.state.all_cust
        register!(s, source, rand(s.state.arrival))
    end
end

function renege(s, c)
    # we do not remove actions from event_queue so we have to check state
    if c.state == WAITING
        s.state.reneged_count += 1
        c.state = RENEGED
        s.state.report && @printf("%6.4f %s: RENEGED after %6.4f\n",
                           s.now, c.name, s.now-c.arrival_time)
    end
end

function serve(s, c)
    # we do not remove actions from event_queue so we have to check state
    if c.state == WAITING
        s.state.report && @printf("%6.4f %s: waited %6.4f\n",
                                  s.now, c.name, s.now-c.arrival_time)
        c.state = CASHING
        register!(s, x -> finish(x, c), rand(s.state.cashing))
    else
        # actually customer reneged so take the next from queue
        request!(s, s.state.counter, serve)
    end
end

function finish(s, c)
    s.state.report && @printf("%6.4f %s: Finished\n", s.now, c.name)
    c.state = FINISHED
    request!(s, s.state.counter, serve)
end

function run(ccount, report)
    bs = BankState(ccount, 0, 0, Exponential(10.0),
                   Uniform(1.0, 3.0), Exponential(12.0),
                   SimQueue{Customer}(), report)
    report && println("Bank renege")
    s = Scheduler(bs)
    request!(s, bs.counter, serve)
    source(s)
    go!(s, Inf)
    println("Renege probability: ", bs.reneged_count / bs.all_cust)
end

Random.seed!(1)
@time run(10, true)
@time run(1_000_000, false)


text:
using Observables
using Test

@testset "ambiguities" begin
    if VERSION < v"1.2"  # Julia 1.0 is bad at detecting ambiguities
    elseif VERSION < v"1.6"
        @test isempty(detect_ambiguities(Base, Observables))
    else
        @test isempty(detect_ambiguities(Observables))
    end
end

@testset "construct and show" begin
    obs = Observable(5)
    @test string(obs) == "Observable{$Int} with 0 listeners. Value:\n5"
    on(identity, obs)
    @test string(obs) == "Observable{$Int} with 1 listeners. Value:\n5"
    on(x->nothing, obs)
    @test string(obs) == "Observable{$Int} with 2 listeners. Value:\n5"
    obs[] = 7
    @test string(obs) == "Observable{$Int} with 2 listeners. Value:\n7"
    obs = Observable{Any}()
    @test string(obs) == "Observable{Any} with 0 listeners. Value:\nnot assigned yet!"
end

@testset "listeners" begin
    r = Observable(0)
    @test r[] == 0
    r[] = 1
    @test r[] == 1

    f = on(r) do x
        r[] == 2
        @test x == 2
    end
    r[] = 2

    @test off(r, f)
    @test !(f in r.listeners)
    @test off(r, f) == false
    r[] = 3 # shouldn't call test
end

@testset "onany and map" begin
    r = Observable(0)
    tested = Ref(false)
    onany(1, r) do x, y
        @test x == 1
        @test y == 2
        tested[] = true
    end
    r[] = 2
    @test tested[]

    r1 = Observable(0)
    r2 = Observable(0)
    map!(x->x+1, r2, r1)
    @test r2[] == 1
    r2 = Observable(0)
    map!(x->x+1, r2, r1; update=false)
    @test r2[] == 0
    r1[] = 1
    @test r2[] == 2

    r1 = Observable(2)
    r2 = @inferred(map(+, r1, 1))
    @test r2[] === 3
    @test eltype(r2) === Int
    r1[] = 3
    @test r2[] == 4

    r3 = @inferred(map!(+, Observable{Float32}(), r1, 1))
    @test r3[] === 4.0f0
    @test eltype(r3) === Float32
    r1[] = 4
    @test r3[] === 5.0f0

    # Make sure `precompile` doesn't error
    precompile(r1)
end

@testset "disconnect observerfuncs" begin

    x = Observable(1)
    y = Observable(2)
    z = Observable(3)

    of1 = on(x) do x
        println(x)
    end

    of2_3 = onany(y, z) do y, z
        println(y, z)
    end

    off(of1)
    off.(of2_3)
    for obs in (x, y, z)
        @test isempty(obs.listeners)
    end
end

# this struct is just supposed to show that a value in memory was released
mutable struct ToFinalize
    val

    function ToFinalize(val, finalized_flag::Ref)
        tf = new(val)
        finalizer(tf) do tf
            finalized_flag[] = true
        end
    end
end

@testset "weak connections" begin
    a = Observable(1)

    finalized_flag = Ref(false)

    function create_a_dangling_listener()
        t = ToFinalize(1, finalized_flag)

        obsfunc = on(a; weak = true) do a
            t.val += 1
        end

        a[] = 2
        @test t.val == 2

        # obsfunc falls out of scope here and should deregister the closure
        # when it gets garbage collected, which should in turn free t
        nothing
    end

    GC.enable(false)
    create_a_dangling_listener()
    @test length(Observables.listeners(a)) == 1
    @test finalized_flag[] == false

    GC.enable(true)
    GC.gc()
    # somehow this needs a double sweep, maybe first obsfunc, then ToFinalize?
    GC.gc()
    @test isempty(Observables.listeners(a))
    @test finalized_flag[] == true
end

@testset "world age" begin
    # issue #50
    obs = Observable(0)
    t = @task for n=1:10
        obs[] = n
        sleep(0.1)
    end
    schedule(t)
    sleep(0.1)
    map(x->2x, obs)
    sleep(0.2)
    @test t.state !== :failed # !istaskfailed(t) is only available on Julia 1.3+
end

@testset "macros" begin
    a = Observable(2)
    b = Observable(3)
    c = Observables.@map &a + &b
    @test c isa Observable
    @test c[] == 5
    a[] = 100
    @test c[] == 103

    a = Observable(2)
    b = Observable(3)
    c = Observable(10)
    Observables.@map! c &a + &b
    @test c[] == 5
    a[] = 100
    @test c[] == 103

    a = Observable(2)
    b = Observable(3)
    c = Observable(10)
    Observables.@on c[] = &a + &b
    @test c[] == 10
    a[] = 100
    @test c[] == 103
end

@testset "observe_changes" begin
    o = Observable(0)
    changes = observe_changes(o)
    changes_approx = observe_changes(o, (x,y) -> abs(x-y) < 2)
    @test eltype(changes) == eltype(changes_approx) == Int
    @test changes[] == changes_approx[] == 0

    values = Int[]
    values_approx = Int[]
    on(changes) do v
        push!(values, v)
    end
    on(changes_approx) do v
        push!(values_approx, v)
    end

    for i in 1:100
        o[] = floor(Int, i/10)
    end
    @test values == 1:10
    @test values_approx == 2:2:10
end

@testset "async_latest" begin
    o = Observable(0)
    cnt = Ref(0)
    function compute_something(x)
        for i=1:10^7; rand() end
        cnt[] = cnt[] + 1
    end
    o_latest = async_latest(o, 1)

    on(compute_something, o_latest) # compute something on the latest update
    for i=1:5
        o[] = i
    end
    sleep(1)

    @test o_latest[] == 5
    @test cnt[] == 1

    for i=5:-1:1
        @async o[] = i
    end
    sleep(1)

    @test o_latest[] == 1
    @test cnt[] == 2 # only one more

    o = Observable(0)
    cnt[] = 0
    function compute_something(x)
        for i=1:10^8; rand() end
        cnt[] = cnt[] + 1
    end
    sleep(1)
    o_latest = async_latest(o, 3)

    on(compute_something, o_latest) # compute something on the latest update
    for i=1:5
        o[] = i
    end
    sleep(.1)

    @test o_latest[] == 5
    @test cnt[] == 3
end

@testset "connect!" begin
    a, b = Observable(""), Observable("")
    connect!(a, b)
    b[] = "Hi!"
    @test a[] == "Hi!"
end

@testset "throttle" begin
    obs = Observable(1)
    throttled = throttle(2, obs)
    @test throttled[] == 1
    obs[] = 2
    sleep(0.1)
    @test throttled[] == 2
    obs[] = 3
    sleep(1)
    @test throttled[] == 2
    sleep(2)
    @test throttled[] == 3
end

struct A<:Observables.AbstractObservable{Int}; end

struct B{T}<:Observables.AbstractObservable{T}
    output::Observable{T}
end

Observables.observe(b::B) = b.output

@testset "interface" begin
    a = A()
    @test_throws ErrorException a[]
    @test_throws ErrorException a[] = 2
    @test_throws ErrorException Observables.obsid(a)
    @test_throws ErrorException Observables.listeners(a)
    @test eltype(a) == Int

    b = B(Observable("test"))
    @test b[] == "test"
    b[] = "test2"
    sleep(0.1)
    @test b[] == "test2"
    @test isempty(Observables.listeners(b))
    @test eltype(b) == String
end

@testset "pair" begin
    v = Observables.ObservablePair(Observable(1.0), f = exp, g = log)
    @test length(v) == 2
    @test firstindex(v) == 1
    @test lastindex(v) == 2
    firstobs, lastobs = v
    @test firstobs[] == 1.0
    @test lastobs[] ≈ ℯ
    @test v.second[] ≈ ℯ
    @test first(v)[] ≈ 1.0
    @test last(v)[] ≈ ℯ
    @test v[2][] ≈ ℯ

    v.first[] = 0
    @test v.second[] ≈ 1
    v.second[] = 2
    @test v.first[] ≈ log(2)

    obs = Observable(Observable(2))
    o2 = Observables.flatten(obs)

    o2[] = 12
    sleep(0.1)
    @test obs[][] == 12
    obs[][] = 22
    sleep(0.1)
    @test o2[] == 22
    obs[] = Observable(11)
    sleep(0.1)
    @test o2[] == 11

    obs = Observable(Observable(Observable(10)))
    o2 = Observables.flatten(obs)
    @test o2[] == 10

    obs[] = Observable(Observable(13))
    @test o2[] == 13
end

@testset "async" begin
    y = Observable(@async begin
        sleep(0.5)
        return 1 + 1
    end)
    @test !isdefined(y, :val)
    while !isdefined(y, :val)
        sleep(0.001)
    end
    @test y[] == 2

    x = Observable(1)
    y = map(x) do val
        @async begin
            sleep(0.5)
            return val + 1
        end
    end
    @test !isdefined(y, :val)
    while !isdefined(y, :val)
        sleep(0.001)
    end
    @test y[] == 2

    x = Observable(1)
    y = map(x) do val
        Channel() do channel
            for i in 1:10
                put!(channel, i + val)
            end
        end
    end
    @test !isdefined(y, :val)
    vals = Int[]
    on(y) do val
        push!(vals, val)
    end
    while length(vals) != 10
        sleep(0.001)
    end
    @test vals == 2:11
    y = Observable(Channel() do channel
        for i in 1:10
            put!(channel, i + 1)
        end
    end)
    @test !isdefined(y, :val)
    vals = Int[]
    on(y) do val
        push!(vals, val)
    end
    while length(vals) != 10
        sleep(0.001)
    end
    @test vals == 2:11
end

@testset "copy" begin
    obs = Observable{Any}(:hey)
    obs_copy = copy(obs)
    obs[] = 1
    @test obs_copy[] == 1
end

@testset "methodlist" begin
    _only(list) = (@assert length(list) == 1; return list[1])  # `only` is avail in Julia 1.4+
    obs = Observable(1)
    obsf = on(obs) do x
        x + 1
    end
    line1 = -2 + @__LINE__
    obs2 = map(obs) do x
        x + 1
    end
    line2 = -2 + @__LINE__
    obsflist = onany(obs, 5) do x, y
        x + y
    end
    line3 = -2 + @__LINE__
    obsf2 = _only(obsflist)
    m = _only(Observables.methodlist(obsf).ms)
    @test occursin("runtests.jl:$line1", string(m))
    m = _only(Observables.methodlist(obsf.f).ms)
    @test occursin("runtests.jl:$line1", string(m))
    m = _only(Observables.methodlist(obs.listeners[2]).ms)
    @test occursin("runtests.jl:$line2", string(m))
    m = _only(Observables.methodlist(obsf2).ms)
    @test occursin("runtests.jl:$line3", string(m))
    obsf3 = on(sqrt, obs)
    @test Observables.methodlist(obsf3).mt.name === :sqrt
end

@testset "deprecations" begin
    a = Observable(1)
    b = @test_deprecated map(sqrt, a; init=13.0)
    @test b[] == 13.0
    a[] = 4
    @test b[] == 2.0
end


text:
using SpecialFunctions

export hann_window, rect_window, kaiser_window, tri_window, lanczos_window

"""
    rect_window(width = 0.5)

Rectangular window which is 1 from 0 to `width` and 0 otherwise.

Can be used to calculate a moving average, for example.
"""
rect_window(width = 0.5) = WindowFunction(x -> 1.0, width)

"""
    hann_window(width = 1)

Hann window (raised cosine) that is 1 at 0 and reaches 0 at `width`.
"""
hann_window(width = 1.0) = WindowFunction(
    x -> cos(pi * x / 2 / width)^2,
    width
)

"""
    kaiser_window(width = 1.2; alpha = 3.5)

Kaiser window cut off at `width` with the shape parameter `alpha`.
"""
kaiser_window(width = 1.2; alpha = 3.5) = WindowFunction(
    x -> besseli(0, pi * alpha * sqrt(1 - (x / width)^2)) / besseli(0, pi * alpha),
    width
)

"""
    tri_window(width = 1)

Triangular window which is 1 at 0 and reaches 0 at `width`.

Can be used to perform linear interpolation when used as upsampling_function
with `width=1`.
"""
tri_window(width=1.) = WindowFunction(
    x -> 1 - x/width,
    width
)

"""
    lanczos_window(lobes = 3; width = 1)

Lanczos window which is 1 at 0 and 0 at 1*`width`, 2*`width`, 3*`width`, ...

The higher the number of lobes, the higher the accuracy of sinc approximation.
"""
function lanczos_window(lobes::Int = 3; width=1.)
    @assert lobes > 0 "Lanczos window is only defined for at least one lobe"

    WindowFunction(
    x -> begin
        x == 0 && return 1.
        x = x/width
        lobes*sin(pi*x)*sin(pi*x/lobes)/(pi^2*x^2)
    end,
    lobes*width
    )
end


text:
# cd social-simulation/Santos2008
# julia src/exec.jl
using LightGraphs: SimpleGraph, ne, nv, barabasi_albert
using DataFrames: DataFrame
using Statistics: mean
using Agents: run!, dummystep
using Dates

include("model.jl")

# η を r に変換する
calc_r(G::SimpleGraph, η::Float64)::Float64 = η * (2 * ne(G) / nv(G) + 1) # 握手の補題: 次数の総和 は 総枝数の2倍に等しい

function run_simulation(network_type::Symbol, cost_model::Symbol, η::Float64; N_gen::Int = 10^5, N_sim::Int = 10)::Float64
    # Log
    (abspath(PROGRAM_FILE) == @__FILE__) && println("network_type = $network_type, cost_model = $cost_model, η = $η")

    f_c = 0.0

    Threads.@threads for _ in 1:N_sim
        G = barabasi_albert(10^3, 2)
        for _ in 1:N_sim
            model = Model.build_model(;G, cost_model, r = calc_r(G, η))
            temp_df, _ = run!(
                model,
                dummystep,
                Model.model_step!,
                N_gen;
                adata = [(:strategy, mean)],
                when = (N_gen - 1999):(N_gen + 1)
            )
            f_c += mean(temp_df.mean_strategy)
        end
    end

    return f_c / N_sim^2
end

if abspath(PROGRAM_FILE) == @__FILE__
    # 10分 x 10 x 7 = 700分 (12時間)
    date_time = Dates.format(now(), "yyyymmdd_HHMMSS")
    file_name = "data/output_$(date_time).csv"

    for cost_model in [:fixed_cost, :variable_cost]
        for η in 0.5:0.05:0.8 # 0.2:0.05:0.8
            @time c_rate = run_simulation(:scale_free, cost_model, η; N_gen = 10^5, N_sim = 7)
            open(file_name, "a") do f
                println(f, "$(:scale_free),$(cost_model),$(η),$(c_rate)")
            end
        end
    end
end